[gcalctool] Simplify modifications of MP.t
- From: Robert Ancell <rancell src gnome org>
- To: svn-commits-list gnome org
- Subject: [gcalctool] Simplify modifications of MP.t
- Date: Sun, 10 May 2009 04:52:51 -0400 (EDT)
commit 56f43ec1177032fdc298e3a1a56d5c27b56a8c04
Author: Robert Ancell <robert ancell gmail com>
Date: Sun May 10 18:36:31 2009 +1000
Simplify modifications of MP.t
---
src/mp-internal.h | 1 -
src/mp-trigonometric.c | 115 +++++++++++++++++++++++++++++++++++++++++++-----
src/mp.c | 95 ---------------------------------------
3 files changed, 104 insertions(+), 107 deletions(-)
diff --git a/src/mp-internal.h b/src/mp-internal.h
index 876d0b3..803072e 100644
--- a/src/mp-internal.h
+++ b/src/mp-internal.h
@@ -53,6 +53,5 @@ void mp_normalize(MPNumber *, int trunc);
void mpexp1(const MPNumber *, MPNumber *);
void mpmulq(const MPNumber *, int, int, MPNumber *);
void mp_reciprocal(const MPNumber *, MPNumber *);
-void mp_atan1N(int n, MPNumber *z);
#endif /* MP_INTERNAL_H */
diff --git a/src/mp-trigonometric.c b/src/mp-trigonometric.c
index 3f97882..05a9d9a 100644
--- a/src/mp-trigonometric.c
+++ b/src/mp-trigonometric.c
@@ -46,7 +46,6 @@ mp_compare_mp_to_int(const MPNumber *x, int i)
return mp_compare_mp_to_mp(x, &t);
}
-
/* COMPUTES Z = SIN(X) IF DO_SIN != 0, Z = COS(X) IF DO_SIN == 0,
* USING TAYLOR SERIES. ASSUMES ABS(X) <= 1.
* X AND Y ARE MP NUMBERS, IS AN INTEGER.
@@ -124,6 +123,96 @@ mpsin1(const MPNumber *x, MPNumber *z, int do_sin)
mp_add_integer(z, 1, z);
}
+/* COMPUTES MP Z = ARCTAN(1/N), ASSUMING INTEGER N > 1.
+ * USES SERIES ARCTAN(X) = X - X**3/3 + X**5/5 - ...
+ * DIMENSION OF R IN CALLING PROGRAM MUST BE
+ * AT LEAST 2T+6
+ * CHECK LEGALITY OF B, T, M AND MXR
+ */
+static void
+mp_atan1N(int n, MPNumber *z)
+{
+ int i, b2, id;
+ MPNumber t2;
+
+ mpchk();
+ if (n <= 1) {
+ mperr("*** N <= 1 IN CALL TO MP_ATAN1N ***");
+ z->sign = 0;
+ return;
+ }
+
+ /* SET SUM TO X = 1/N */
+ mp_set_from_fraction(1, n, z);
+
+ /* SET ADDITIVE TERM TO X */
+ mp_set_from_mp(z, &t2);
+
+ /* ASSUME AT LEAST 16-BIT WORD. */
+ b2 = max(MP.b, 64);
+ if (n < b2)
+ id = b2 * 7 * b2 / (n * n);
+ else
+ id = 0;
+
+ /* MAIN LOOP. FIRST REDUCE T IF POSSIBLE */
+ for (i = 1; ; i += 2) {
+ int t, ts;
+
+ t = MP.t + 2 + t2.exponent - z->exponent;
+ if (t <= 1)
+ break;
+ t = min(t, MP.t);
+
+ /* IF (I+2)*N**2 IS NOT REPRESENTABLE AS AN INTEGER THE DIVISION
+ * FOLLOWING HAS TO BE PERFORMED IN SEVERAL STEPS.
+ */
+ ts = MP.t;
+ MP.t = t;
+ if (i >= id) {
+ mpmulq(&t2, -i, i + 2, &t2);
+ mp_divide_integer(&t2, n, &t2);
+ mp_divide_integer(&t2, n, &t2);
+ }
+ else {
+ mpmulq(&t2, -i, (i + 2)*n*n, &t2);
+ }
+ MP.t = ts;
+
+ /* ADD TO SUM */
+ mp_add(&t2, z, z);
+ if (t2.sign == 0)
+ break;
+ }
+}
+
+/* SETS MP Z = PI TO THE AVAILABLE PRECISION.
+ * USES PI/4 = 4.ARCTAN(1/5) - ARCTAN(1/239).
+ * TIME IS O(T**2).
+ * DIMENSION OF R MUST BE AT LEAST 3T+8
+ * CHECK LEGALITY OF B, T, M AND MXR
+ */
+void
+mp_get_pi(MPNumber *z)
+{
+ float prec;
+ MPNumber t;
+
+ mpchk();
+
+ mp_atan1N(5, &t);
+ mp_multiply_integer(&t, 4, &t);
+ mp_atan1N(239, z);
+ mp_subtract(&t, z, z);
+ mp_multiply_integer(z, 4, z);
+
+ /* RETURN IF ERROR IS LESS THAN 0.01 */
+ prec = fabs(mp_cast_to_float(z) - 3.1416);
+ if (prec < 0.01) return;
+
+ /* FOLLOWING MESSAGE MAY INDICATE THAT B**(T-1) IS TOO SMALL */
+ mperr("*** ERROR OCCURRED IN MP_GET_PI, RESULT INCORRECT ***");
+}
/* MP precision arc cosine.
*
@@ -271,7 +360,7 @@ mp_asinh(const MPNumber *x, MPNumber *z)
void
mp_atan(const MPNumber *x, MPNumber *z)
{
- int i, q, ts;
+ int i, q;
float rx = 0.0, ry;
MPNumber t1, t2;
@@ -285,9 +374,8 @@ mp_atan(const MPNumber *x, MPNumber *z)
if (abs(x->exponent) <= 2)
rx = mp_cast_to_float(x);
- q = 1;
-
/* REDUCE ARGUMENT IF NECESSARY BEFORE USING SERIES */
+ q = 1;
while (t2.exponent >= 0)
{
if (t2.exponent == 0 && (t2.fraction[0] + 1) << 1 <= MP.b)
@@ -304,23 +392,28 @@ mp_atan(const MPNumber *x, MPNumber *z)
/* USE POWER SERIES NOW ARGUMENT IN (-0.5, 0.5) */
mp_set_from_mp(&t2, z);
mp_multiply(&t2, &t2, &t1);
- i = 1;
- ts = MP.t;
/* SERIES LOOP. REDUCE T IF POSSIBLE. */
- while ((MP.t = ts + 2 + t2.exponent) > 1) {
- MP.t = min(MP.t,ts);
+ for (i = 1; ; i += 2) {
+ int t, ts;
+
+ t = MP.t + 2 + t2.exponent;
+ if (t <= 1)
+ break;
+ t = min(t, MP.t);
+
+ ts = MP.t;
+ MP.t = t;
mp_multiply(&t2, &t1, &t2);
mpmulq(&t2, -i, i + 2, &t2);
- i += 2;
MP.t = ts;
+
mp_add(z, &t2, z);
if (t2.sign == 0)
break;
}
- /* RESTORE T, CORRECT FOR ARGUMENT REDUCTION, AND EXIT */
- MP.t = ts;
+ /* CORRECT FOR ARGUMENT REDUCTION, AND EXIT */
mp_multiply_integer(z, q, z);
/* CHECK THAT RELATIVE ERROR LESS THAN 0.01 UNLESS EXPONENT
diff --git a/src/mp.c b/src/mp.c
index 2461558..f85ca80 100644
--- a/src/mp.c
+++ b/src/mp.c
@@ -372,71 +372,6 @@ mp_add_fraction(const MPNumber *x, int i, int j, MPNumber *y)
}
-/* COMPUTES MP Z = ARCTAN(1/N), ASSUMING INTEGER N > 1.
- * USES SERIES ARCTAN(X) = X - X**3/3 + X**5/5 - ...
- * DIMENSION OF R IN CALLING PROGRAM MUST BE
- * AT LEAST 2T+6
- * CHECK LEGALITY OF B, T, M AND MXR
- */
-void
-mp_atan1N(int n, MPNumber *z)
-{
- int i, b2, id, ts;
- MPNumber t;
-
- mpchk();
- if (n <= 1) {
- mperr("*** N <= 1 IN CALL TO MP_ATAN1N ***");
- z->sign = 0;
- return;
- }
-
- ts = MP.t;
-
- /* SET SUM TO X = 1/N */
- mp_set_from_fraction(1, n, z);
-
- /* SET ADDITIVE TERM TO X */
- mp_set_from_mp(z, &t);
- i = 1;
- id = 0;
-
- /* ASSUME AT LEAST 16-BIT WORD. */
- b2 = max(MP.b, 64);
- if (n < b2)
- id = b2 * 7 * b2 / (n * n);
-
- /* MAIN LOOP. FIRST REDUCE T IF POSSIBLE */
- while ((MP.t = ts + 2 + t.exponent - z->exponent) > 1) {
-
- MP.t = min(MP.t,ts);
-
- /* IF (I+2)*N**2 IS NOT REPRESENTABLE AS AN INTEGER THE DIVISION
- * FOLLOWING HAS TO BE PERFORMED IN SEVERAL STEPS.
- */
- if (i >= id) {
- mpmulq(&t, -i, i + 2, &t);
- mp_divide_integer(&t, n, &t);
- mp_divide_integer(&t, n, &t);
- }
- else {
- mpmulq(&t, -i, (i + 2)*n*n, &t);
- }
-
- i += 2;
-
- /* RESTORE T */
- MP.t = ts;
-
- /* ADD TO SUM */
- mp_add(&t, z, z);
- if (t.sign == 0)
- break;
- }
- MP.t = ts;
-}
-
-
/* CHECKS LEGALITY OF B, T, M AND MXR.
* THE CONDITION ON MXR (THE DIMENSION OF R IN COMMON) IS THAT
* MXR >= (I*T + J)
@@ -1809,36 +1744,6 @@ mpovfl(MPNumber *x, const char *error)
mperr("*** CALL TO MPOVFL, MP OVERFLOW OCCURRED ***");
}
-
-/* SETS MP Z = PI TO THE AVAILABLE PRECISION.
- * USES PI/4 = 4.ARCTAN(1/5) - ARCTAN(1/239).
- * TIME IS O(T**2).
- * DIMENSION OF R MUST BE AT LEAST 3T+8
- * CHECK LEGALITY OF B, T, M AND MXR
- */
-void
-mp_get_pi(MPNumber *z)
-{
- float prec;
- MPNumber t;
-
- mpchk();
-
- mp_atan1N(5, &t);
- mp_multiply_integer(&t, 4, &t);
- mp_atan1N(239, z);
- mp_subtract(&t, z, z);
- mp_multiply_integer(z, 4, z);
-
- /* RETURN IF ERROR IS LESS THAN 0.01 */
- prec = fabs(mp_cast_to_float(z) - 3.1416);
- if (prec < 0.01) return;
-
- /* FOLLOWING MESSAGE MAY INDICATE THAT B**(T-1) IS TOO SMALL */
- mperr("*** ERROR OCCURRED IN MP_GET_PI, RESULT INCORRECT ***");
-}
-
-
/* RETURNS Y = X**N, FOR MP X AND Y, INTEGER N, WITH 0**0 = 1.
* R MUST BE DIMENSIONED AT LEAST 4T+10 IN CALLING PROGRAM
* (2T+6 IS ENOUGH IF N NONNEGATIVE).
[
Date Prev][
Date Next] [
Thread Prev][
Thread Next]
[
Thread Index]
[
Date Index]
[
Author Index]