[gcalctool/gcalctool-newui2] ...
- From: Robert Ancell <rancell src gnome org>
- To: svn-commits-list gnome org
- Subject: [gcalctool/gcalctool-newui2] ...
- Date: Sun, 12 Jul 2009 04:57:31 +0000 (UTC)
commit cf06af78e6d44555982281a7390aa91af38f685b
Author: Robert Ancell <robert ancell gmail com>
Date: Sun Jul 12 12:56:59 2009 +0800
...
src/mp-convert.c | 48 ++++----
src/mp-internal.h | 18 ++-
src/mp-trigonometric.c | 18 ++--
src/mp.c | 289 ++++++++++++++++++++++--------------------------
4 files changed, 176 insertions(+), 197 deletions(-)
---
diff --git a/src/mp-convert.c b/src/mp-convert.c
index 9285fd3..eccaf4e 100644
--- a/src/mp-convert.c
+++ b/src/mp-convert.c
@@ -48,7 +48,7 @@ mp_set_from_float(float rx, MPNumber *z)
int i, k, i2, ib, ie, tp;
float rb, rj;
- i2 = MP.t + 4;
+ i2 = MP_T + 4;
memset(z, 0, sizeof(MPNumber));
@@ -82,7 +82,7 @@ mp_set_from_float(float rx, MPNumber *z)
* SET EXPONENT TO 0
*/
z->exponent = 0;
- rb = (float) MP.b;
+ rb = (float) MP_BASE;
/* CONVERSION LOOP (ASSUME SINGLE-PRECISION OPS. EXACT) */
for (i = 0; i < i2; i++) {
@@ -95,7 +95,7 @@ mp_set_from_float(float rx, MPNumber *z)
mp_normalize(z, 0);
/* Computing MAX */
- ib = max(MP.b * 7 * MP.b, 32767) / 16;
+ ib = max(MP_BASE * 7 * MP_BASE, 32767) / 16;
tp = 1;
/* NOW MULTIPLY BY 16**IE */
@@ -103,7 +103,7 @@ mp_set_from_float(float rx, MPNumber *z)
k = -ie;
for (i = 1; i <= k; ++i) {
tp <<= 4;
- if (tp <= ib && tp != MP.b && i < k)
+ if (tp <= ib && tp != MP_BASE && i < k)
continue;
mp_divide_integer(z, tp, z);
tp = 1;
@@ -111,7 +111,7 @@ mp_set_from_float(float rx, MPNumber *z)
} else if (ie > 0) {
for (i = 1; i <= ie; ++i) {
tp <<= 4;
- if (tp <= ib && tp != MP.b && i < ie)
+ if (tp <= ib && tp != MP_BASE && i < ie)
continue;
mp_multiply_integer(z, tp, z);
tp = 1;
@@ -140,7 +140,7 @@ mp_set_from_double(double dx, MPNumber *z)
int i, k, i2, ib, ie, tp;
double db, dj;
- i2 = MP.t + 4;
+ i2 = MP_T + 4;
memset(z, 0, sizeof(MPNumber));
@@ -168,7 +168,7 @@ mp_set_from_double(double dx, MPNumber *z)
*/
z->exponent = 0;
- db = (double) MP.b;
+ db = (double) MP_BASE;
/* CONVERSION LOOP (ASSUME DOUBLE-PRECISION OPS. EXACT) */
for (i = 0; i < i2; i++) {
@@ -181,7 +181,7 @@ mp_set_from_double(double dx, MPNumber *z)
mp_normalize(z, 0);
/* Computing MAX */
- ib = max(MP.b * 7 * MP.b, 32767) / 16;
+ ib = max(MP_BASE * 7 * MP_BASE, 32767) / 16;
tp = 1;
/* NOW MULTIPLY BY 16**IE */
@@ -189,7 +189,7 @@ mp_set_from_double(double dx, MPNumber *z)
k = -ie;
for (i = 1; i <= k; ++i) {
tp <<= 4;
- if (tp <= ib && tp != MP.b && i < k)
+ if (tp <= ib && tp != MP_BASE && i < k)
continue;
mp_divide_integer(z, tp, z);
tp = 1;
@@ -197,7 +197,7 @@ mp_set_from_double(double dx, MPNumber *z)
} else if (ie > 0) {
for (i = 1; i <= ie; ++i) {
tp <<= 4;
- if (tp <= ib && tp != MP.b && i < ie)
+ if (tp <= ib && tp != MP_BASE && i < ie)
continue;
mp_multiply_integer(z, tp, z);
tp = 1;
@@ -227,13 +227,13 @@ mp_set_from_integer(int ix, MPNumber *z)
z->sign = 1;
/* SET EXPONENT TO T */
- z->exponent = MP.t;
+ z->exponent = MP_T;
/* CLEAR FRACTION */
- memset(z->fraction, 0, (MP.t-1)*sizeof(int));
+ memset(z->fraction, 0, (MP_T-1)*sizeof(int));
/* INSERT IX */
- z->fraction[MP.t - 1] = ix;
+ z->fraction[MP_T - 1] = ix;
/* NORMALIZE BY CALLING MPMUL2 */
mpmul2(z, 1, z, 1);
@@ -285,8 +285,8 @@ mp_cast_to_int(const MPNumber *x)
for (i = 0; i < x2; i++) {
int izs;
izs = ret_val;
- ret_val = MP.b * ret_val;
- if (i < MP.t)
+ ret_val = MP_BASE * ret_val;
+ if (i < MP_T)
ret_val += x->fraction[i];
/* CHECK FOR SIGNS OF INTEGER OVERFLOW */
@@ -301,11 +301,11 @@ mp_cast_to_int(const MPNumber *x)
for (i = x2 - 1; i >= 0; i--) {
int j1, kx;
- j1 = j / MP.b;
+ j1 = j / MP_BASE;
kx = 0;
- if (i < MP.t)
+ if (i < MP_T)
kx = x->fraction[i];
- if (kx != j - MP.b * j1)
+ if (kx != j - MP_BASE * j1)
return 0;
j = j1;
}
@@ -363,8 +363,8 @@ mp_cast_to_float(const MPNumber *x)
if (x->sign == 0)
return 0.0;
- rb = (float) MP.b;
- for (i = 0; i < MP.t; i++) {
+ rb = (float) MP_BASE;
+ for (i = 0; i < MP_T; i++) {
rz = rb * rz + (float)x->fraction[i];
/* CHECK IF FULL SINGLE-PRECISION ACCURACY ATTAINED */
@@ -378,7 +378,7 @@ mp_cast_to_float(const MPNumber *x)
/* CHECK REASONABLENESS OF RESULT */
/* LHS SHOULD BE <= 0.5, BUT ALLOW FOR SOME ERROR IN ALOG */
if (rz <= (float)0. ||
- fabs((float) x->exponent - (log(rz) / log((float) MP.b) + (float).5)) > (float).6) {
+ fabs((float) x->exponent - (log(rz) / log((float) MP_BASE) + (float).5)) > (float).6) {
/* FOLLOWING MESSAGE INDICATES THAT X IS TOO LARGE OR SMALL -
* TRY USING MPCMRE INSTEAD.
*/
@@ -429,8 +429,8 @@ mp_cast_to_double(const MPNumber *x)
if (x->sign == 0)
return 0.0;
- db = (double) MP.b;
- for (i = 0; i < MP.t; i++) {
+ db = (double) MP_BASE;
+ for (i = 0; i < MP_T; i++) {
ret_val = db * ret_val + (double) x->fraction[i];
tm = i;
@@ -451,7 +451,7 @@ mp_cast_to_double(const MPNumber *x)
/* LHS SHOULD BE .LE. 0.5 BUT ALLOW FOR SOME ERROR IN DLOG */
if (ret_val <= 0. ||
((d__1 = (double) ((float) x->exponent) - (log(ret_val) / log((double)
- ((float) MP.b)) + .5), abs(d__1)) > .6)) {
+ ((float) MP_BASE)) + .5), abs(d__1)) > .6)) {
/* FOLLOWING MESSAGE INDICATES THAT X IS TOO LARGE OR SMALL -
* TRY USING MPCMDE INSTEAD.
*/
diff --git a/src/mp-internal.h b/src/mp-internal.h
index 80f8de4..f6ba4ea 100644
--- a/src/mp-internal.h
+++ b/src/mp-internal.h
@@ -34,17 +34,21 @@
/* Evil global variables that must be removed */
struct {
- /* Base */
- int b;
-
- /* Number of digits */
- // This number is temporarily changed across calls to add/subtract/multiply/divide - it should be passed to those calls
- int t;
-
/* Min/max exponent value */
+ // FIXME: MP.m modified inside functions, needs to be fixed to be thread safe
int m;
} MP;
+#define MP_BASE 10000
+
+//2E0 BELOW ENSURES AT LEAST ONE GUARD DIGIT
+//MP.t = (int) ((float) (accuracy) * log((float)10.) / log((float) MP.b) + (float) 2.0);
+//if (MP.t > MP_SIZE) {
+// mperr("MP_SIZE TOO SMALL IN CALL TO MPSET, INCREASE MP_SIZE AND DIMENSIONS OF MP ARRAYS TO AT LEAST %d ***", MP.t);
+// MP.t = MP_SIZE;
+//}
+#define MP_T 55
+
void mperr(const char *format, ...) __attribute__((format(printf, 1, 2)));
void mpgcd(int *, int *);
void mpmul2(const MPNumber *, int, MPNumber *, int);
diff --git a/src/mp-trigonometric.c b/src/mp-trigonometric.c
index a439ba5..e809c82 100644
--- a/src/mp-trigonometric.c
+++ b/src/mp-trigonometric.c
@@ -70,7 +70,7 @@ mpsin1(const MPNumber *x, MPNumber *z, int do_sin)
return;
}
- b2 = max(MP.b, 64) << 1;
+ b2 = max(MP_BASE, 64) << 1;
mp_multiply(x, x, &t2);
if (mp_compare_mp_to_int(&t2, 1) > 0) {
mperr("*** ABS(X) > 1 IN CALL TO MPSIN1 ***");
@@ -91,7 +91,7 @@ mpsin1(const MPNumber *x, MPNumber *z, int do_sin)
/* POWER SERIES LOOP. REDUCE T IF POSSIBLE */
for (; ; i += 2) {
- if (MP.t + t1.exponent <= 0)
+ if (MP_T + t1.exponent <= 0)
break;
/* IF I*(I+1) IS NOT REPRESENTABLE AS AN INTEGER, THE FOLLOWING
@@ -139,7 +139,7 @@ mp_atan1N(int n, MPNumber *z)
mp_set_from_mp(z, &t2);
/* ASSUME AT LEAST 16-BIT WORD. */
- b2 = max(MP.b, 64);
+ b2 = max(MP_BASE, 64);
if (n < b2)
id = b2 * 7 * b2 / (n * n);
else
@@ -147,7 +147,7 @@ mp_atan1N(int n, MPNumber *z)
/* MAIN LOOP. FIRST REDUCE T IF POSSIBLE */
for (i = 1; ; i += 2) {
- if (MP.t + 2 + t2.exponent - z->exponent <= 1)
+ if (MP_T + 2 + t2.exponent - z->exponent <= 1)
break;
/* IF (I+2)*N**2 IS NOT REPRESENTABLE AS AN INTEGER THE DIVISION
@@ -514,7 +514,7 @@ mp_atan(const MPNumber *x, MPAngleUnit unit, MPNumber *z)
q = 1;
while (t2.exponent >= 0)
{
- if (t2.exponent == 0 && (t2.fraction[0] + 1) << 1 <= MP.b)
+ if (t2.exponent == 0 && (t2.fraction[0] + 1) << 1 <= MP_BASE)
break;
q <<= 1;
@@ -531,7 +531,7 @@ mp_atan(const MPNumber *x, MPAngleUnit unit, MPNumber *z)
/* SERIES LOOP. REDUCE T IF POSSIBLE. */
for (i = 1; ; i += 2) {
- if (MP.t + 2 + t2.exponent <= 1)
+ if (MP_T + 2 + t2.exponent <= 1)
break;
mp_multiply(&t2, &t1, &t2);
@@ -657,7 +657,7 @@ mp_cosh(const MPNumber *x, MPNumber *z)
/* RESTORE M. IF RESULT OVERFLOWS OR UNDERFLOWS, mp_divide_integer WILL
* ACT ACCORDINGLY.
*/
- MP.m += -2;
+ MP.m -= 2;
mp_divide_integer(z, 2, z);
}
@@ -701,7 +701,7 @@ mp_sinh(const MPNumber *x, MPNumber *z)
/* RESTORE M. IF RESULT OVERFLOWS OR UNDERFLOWS, mp_divide_integer AT
* STATEMENT 30 WILL ACT ACCORDINGLY.
*/
- MP.m += -2;
+ MP.m -= 2;
}
/* DIVIDE BY TWO AND RESTORE SIGN */
@@ -730,7 +730,7 @@ mp_tanh(const MPNumber *x, MPNumber *z)
mp_abs(x, &t);
/* SEE IF ABS(X) SO LARGE THAT RESULT IS +-1 */
- r__1 = (float) MP.t * 0.5 * log((float) MP.b);
+ r__1 = (float) MP_T * 0.5 * log((float) MP_BASE);
mp_set_from_float(r__1, z);
if (mp_compare_mp_to_mp(&t, z) > 0) {
/* HERE ABS(X) IS VERY LARGE */
diff --git a/src/mp.c b/src/mp.c
index f2247c2..45e9a53 100644
--- a/src/mp.c
+++ b/src/mp.c
@@ -28,8 +28,6 @@
#include "mp-internal.h"
#include "calctool.h" // FIXME: Required for doerr() and MAXLINE
-// FIXME: MP.t and MP.m modified inside functions, needs to be fixed to be thread safe
-
/* SETS X TO THE LARGEST POSSIBLE POSITIVE MP NUMBER
* CHECK LEGALITY OF B, T, M AND MXR
*/
@@ -38,10 +36,10 @@ mpmaxr(MPNumber *x)
{
int i, it;
- it = MP.b - 1;
+ it = MP_BASE - 1;
/* SET FRACTION DIGITS TO B-1 */
- for (i = 0; i < MP.t; i++)
+ for (i = 0; i < MP_T; i++)
x->fraction[i] = it;
/* SET SIGN AND EXPONENT */
@@ -89,27 +87,6 @@ mpunfl(MPNumber *z)
}
-static int
-pow_ii(int x, int n)
-{
- int p = 1;
-
- if (n <= 0)
- return 1;
-
- for (;;) {
- if (n & 01)
- p *= x;
- if (n >>= 1)
- x *= x;
- else
- break;
- }
-
- return p;
-}
-
-
/* ROUTINE CALLED BY MP_DIVIDE AND MP_SQRT TO ENSURE THAT
* RESULTS ARE REPRESENTED EXACTLY IN T-2 DIGITS IF THEY
* CAN BE. X IS AN MP NUMBER, I AND J ARE INTEGERS.
@@ -119,26 +96,26 @@ mpext(int i, int j, MPNumber *x)
{
int q, s;
- if (x->sign == 0 || MP.t <= 2 || i == 0)
+ if (x->sign == 0 || MP_T <= 2 || i == 0)
return;
/* COMPUTE MAXIMUM POSSIBLE ERROR IN THE LAST PLACE */
q = (j + 1) / i + 1;
- s = MP.b * x->fraction[MP.t - 2] + x->fraction[MP.t - 1];
+ s = MP_BASE * x->fraction[MP_T - 2] + x->fraction[MP_T - 1];
/* SET LAST TWO DIGITS TO ZERO */
if (s <= q) {
- x->fraction[MP.t - 2] = 0;
- x->fraction[MP.t - 1] = 0;
+ x->fraction[MP_T - 2] = 0;
+ x->fraction[MP_T - 1] = 0;
return;
}
- if (s + q < MP.b * MP.b)
+ if (s + q < MP_BASE * MP_BASE)
return;
/* ROUND UP HERE */
- x->fraction[MP.t - 2] = MP.b - 1;
- x->fraction[MP.t - 1] = MP.b;
+ x->fraction[MP_T - 2] = MP_BASE - 1;
+ x->fraction[MP_T - 1] = MP_BASE;
/* NORMALIZE X (LAST DIGIT B IS OK IN MP_MULTIPLY_INTEGER) */
mp_multiply_integer(x, 1, x);
@@ -170,7 +147,8 @@ mpext(int i, int j, MPNumber *x)
void
mp_init(int accuracy)
{
- int i, k, w, b;
+ int i, k, w;
+ //int b, n;
/* DETERMINE LARGE REPRESENTABLE INTEGER W OF FORM 2**K - 1 */
/* ON CYBER 76 HAVE TO FIND K <= 47, SO ONLY LOOP
@@ -199,29 +177,26 @@ mp_init(int accuracy)
/* SET MAXIMUM EXPONENT TO (W-1)/4 */
MP.m = w / 4;
- if (accuracy <= 0) {
- mperr("*** ACCURACY <= 0 IN CALL TO MPSET ***");
- return;
- }
/* B IS THE LARGEST POWER OF 2 SUCH THAT (8*B*B-1) <= W */
- MP.b = pow_ii(2, (k - 3) / 2);
+ /*MP.b = 1;
+ n = (k - 3) / 2;
+ if (n > 0) {
+ for (;;) {
+ if (n & 01)
+ MP.b *= x;
+ if (n >>= 1)
+ x *= x;
+ else
+ break;
+ }
+ }*/
/* Make a multiple of 10 so fractions can be represented exactly */
- b = 1;
+ /*b = 1;
while (MP.b % (10 * b) != MP.b)
b *= 10;
- MP.b = b;
-
- /* 2E0 BELOW ENSURES AT LEAST ONE GUARD DIGIT */
- MP.t = (int) ((float) (accuracy) * log((float)10.) / log((float) MP.b) +
- (float) 2.0);
-
- /* SEE IF T TOO LARGE FOR DIMENSION STATEMENTS */
- if (MP.t > MP_SIZE) {
- mperr("MP_SIZE TOO SMALL IN CALL TO MPSET, INCREASE MP_SIZE AND DIMENSIONS OF MP ARRAYS TO AT LEAST %d ***", MP.t);
- MP.t = MP_SIZE;
- }
+ MP.b = b;*/
}
@@ -253,24 +228,24 @@ mp_add3(const MPNumber *x, const MPNumber *y, int *r, int s, int med)
/* CLEAR GUARD DIGITS TO RIGHT OF X DIGITS */
for(i = 3; i >= med; i--)
- r[MP.t + i] = 0;
+ r[MP_T + i] = 0;
if (s >= 0) {
/* HERE DO ADDITION, EXPONENT(Y) >= EXPONENT(X) */
- for (i = MP.t + 3; i >= MP.t; i--)
+ for (i = MP_T + 3; i >= MP_T; i--)
r[i] = x->fraction[i - med];
c = 0;
for (; i >= med; i--) {
c = y->fraction[i] + x->fraction[i - med] + c;
- if (c < MP.b) {
+ if (c < MP_BASE) {
/* NO CARRY GENERATED HERE */
r[i] = c;
c = 0;
} else {
/* CARRY GENERATED HERE */
- r[i] = c - MP.b;
+ r[i] = c - MP_BASE;
c = 1;
}
}
@@ -278,7 +253,7 @@ mp_add3(const MPNumber *x, const MPNumber *y, int *r, int s, int med)
for (; i >= 0; i--)
{
c = y->fraction[i] + c;
- if (c < MP.b) {
+ if (c < MP_BASE) {
r[i] = c;
i--;
@@ -295,7 +270,7 @@ mp_add3(const MPNumber *x, const MPNumber *y, int *r, int s, int med)
/* MUST SHIFT RIGHT HERE AS CARRY OFF END */
if (c != 0) {
- for (i = MP.t + 3; i > 0; i--)
+ for (i = MP_T + 3; i > 0; i--)
r[i] = r[i - 1];
r[0] = 1;
return 1;
@@ -305,7 +280,7 @@ mp_add3(const MPNumber *x, const MPNumber *y, int *r, int s, int med)
}
c = 0;
- for (i = MP.t + med - 1; i >= MP.t; i--) {
+ for (i = MP_T + med - 1; i >= MP_T; i--) {
/* HERE DO SUBTRACTION, ABS(Y) > ABS(X) */
r[i] = c - x->fraction[i - med];
c = 0;
@@ -313,7 +288,7 @@ mp_add3(const MPNumber *x, const MPNumber *y, int *r, int s, int med)
/* BORROW GENERATED HERE */
if (r[i] < 0) {
c = -1;
- r[i] += MP.b;
+ r[i] += MP_BASE;
}
}
@@ -325,7 +300,7 @@ mp_add3(const MPNumber *x, const MPNumber *y, int *r, int s, int med)
c = 0;
} else {
/* BORROW GENERATED HERE */
- r[i] = c + MP.b;
+ r[i] = c + MP_BASE;
c = -1;
}
}
@@ -344,7 +319,7 @@ mp_add3(const MPNumber *x, const MPNumber *y, int *r, int s, int med)
return 0;
}
- r[i] = c + MP.b;
+ r[i] = c + MP_BASE;
c = -1;
}
@@ -389,7 +364,7 @@ mp_add2(const MPNumber *x, const MPNumber *y, int y_sign, MPNumber *z)
med = abs(exp_diff);
if (exp_diff < 0) {
/* HERE EXPONENT(Y) > EXPONENT(X) */
- if (med > MP.t) {
+ if (med > MP_T) {
/* 'y' so much larger than 'x' that 'x+-y = y'. Warning: still true with rounding?? */
mp_set_from_mp(y, z);
z->sign = y_sign;
@@ -398,7 +373,7 @@ mp_add2(const MPNumber *x, const MPNumber *y, int y_sign, MPNumber *z)
x_largest = 0;
} else if (exp_diff > 0) {
/* HERE EXPONENT(X) > EXPONENT(Y) */
- if (med > MP.t) {
+ if (med > MP_T) {
/* 'x' so much larger than 'y' that 'x+-y = x'. Warning: still true with rounding?? */
mp_set_from_mp(x, z);
return;
@@ -409,7 +384,7 @@ mp_add2(const MPNumber *x, const MPNumber *y, int y_sign, MPNumber *z)
if (sign_prod < 0) {
/* Signs are not equal. find out which mantissa is larger. */
int j;
- for (j = 0; j < MP.t; j++) {
+ for (j = 0; j < MP_T; j++) {
int i = x->fraction[j] - y->fraction[j];
if (i == 0)
continue;
@@ -421,7 +396,7 @@ mp_add2(const MPNumber *x, const MPNumber *y, int y_sign, MPNumber *z)
}
/* Both mantissas equal, so result is zero. */
- if (j >= MP.t) {
+ if (j >= MP_T) {
mp_set_from_integer(0, z);
return;
}
@@ -489,7 +464,7 @@ mp_fractional_component(const MPNumber *x, MPNumber *z)
{
/* RETURN 0 IF X = 0
OR IF EXPONENT SO LARGE THAT NO FRACTIONAL PART */
- if (x->sign == 0 || x->exponent >= MP.t) {
+ if (x->sign == 0 || x->exponent >= MP_T) {
mp_set_from_integer(0, z);
return;
}
@@ -506,8 +481,8 @@ mp_fractional_component(const MPNumber *x, MPNumber *z)
memset(z->fraction, 0, z->exponent*sizeof(int));
if (x != z)
memcpy(z->fraction + z->exponent, x->fraction + x->exponent,
- (MP.t - x->exponent)*sizeof(int));
- memset(z->fraction + MP.t, 0, 4*sizeof(int));
+ (MP_T - x->exponent)*sizeof(int));
+ memset(z->fraction + MP_T, 0, 4*sizeof(int));
/* NORMALIZE RESULT AND RETURN */
mp_normalize(z, 1);
@@ -528,7 +503,7 @@ mp_integer_component(const MPNumber *x, MPNumber *z)
return;
/* IF EXPONENT LARGE ENOUGH RETURN Y = X */
- if (z->exponent >= MP.t) {
+ if (z->exponent >= MP_T) {
return;
}
@@ -539,7 +514,7 @@ mp_integer_component(const MPNumber *x, MPNumber *z)
}
/* SET FRACTION TO ZERO */
- for (i = z->exponent; i < MP.t; i++) {
+ for (i = z->exponent; i < MP_T; i++) {
z->fraction[i] = 0;
}
}
@@ -590,7 +565,7 @@ mp_compare_mp_to_mp(const MPNumber *x, const MPNumber *y)
/* ABS(X) > ABS(Y) */
return x->sign;
}
- for (i = 0; i < MP.t; i++) {
+ for (i = 0; i < MP_T; i++) {
i2 = x->fraction[i] - y->fraction[i];
if (i2 < 0) {
/* ABS(X) < ABS(Y) */
@@ -647,7 +622,7 @@ mp_divide(const MPNumber *x, const MPNumber *y, MPNumber *z)
mpext(i, iz3, z);
/* RESTORE M, CORRECT EXPONENT AND RETURN */
- MP.m += -2;
+ MP.m -= 2;
z->exponent += ie;
/* Check for overflow or underflow */
@@ -684,7 +659,7 @@ mp_divide_integer(const MPNumber *x, int iy, MPNumber *z)
z->exponent = x->exponent;
/* CHECK FOR DIVISION BY B */
- if (iy == MP.b) {
+ if (iy == MP_BASE) {
if (x->exponent <= -MP.m)
{
mpunfl(z);
@@ -704,7 +679,7 @@ mp_divide_integer(const MPNumber *x, int iy, MPNumber *z)
}
c = 0;
- i2 = MP.t + 4;
+ i2 = MP_T + 4;
i = 0;
/* IF IY*B NOT REPRESENTABLE AS AN INTEGER HAVE TO SIMULATE
@@ -712,12 +687,12 @@ mp_divide_integer(const MPNumber *x, int iy, MPNumber *z)
*/
/* Computing MAX */
- b2 = max(MP.b << 3, 32767 / MP.b);
+ b2 = max(MP_BASE << 3, 32767 / MP_BASE);
if (iy < b2) {
/* LOOK FOR FIRST NONZERO DIGIT IN QUOTIENT */
do {
- c = MP.b * c;
- if (i < MP.t)
+ c = MP_BASE * c;
+ if (i < MP_T)
c += x->fraction[i];
i++;
r1 = c / iy;
@@ -728,14 +703,14 @@ mp_divide_integer(const MPNumber *x, int iy, MPNumber *z)
/* ADJUST EXPONENT AND GET T+4 DIGITS IN QUOTIENT */
z->exponent += 1 - i;
z->fraction[0] = r1;
- c = MP.b * (c - iy * r1);
+ c = MP_BASE * (c - iy * r1);
kh = 1;
- if (i < MP.t) {
- kh = MP.t + 1 - i;
+ if (i < MP_T) {
+ kh = MP_T + 1 - i;
for (k = 1; k < kh; k++) {
c += x->fraction[i];
z->fraction[k] = c / iy;
- c = MP.b * (c - iy * z->fraction[k]);
+ c = MP_BASE * (c - iy * z->fraction[k]);
i++;
}
if (c < 0)
@@ -744,7 +719,7 @@ mp_divide_integer(const MPNumber *x, int iy, MPNumber *z)
for (k = kh; k < i2; k++) {
z->fraction[k] = c / iy;
- c = MP.b * (c - iy * z->fraction[k]);
+ c = MP_BASE * (c - iy * z->fraction[k]);
}
if (c < 0)
goto L210;
@@ -756,15 +731,15 @@ mp_divide_integer(const MPNumber *x, int iy, MPNumber *z)
}
/* HERE NEED SIMULATED DOUBLE-PRECISION DIVISION */
- j1 = iy / MP.b;
+ j1 = iy / MP_BASE;
/* LOOK FOR FIRST NONZERO DIGIT */
c2 = 0;
- j2 = iy - j1 * MP.b;
+ j2 = iy - j1 * MP_BASE;
do {
- c = MP.b * c + c2;
+ c = MP_BASE * c + c2;
i__1 = c - j1;
- c2 = i < MP.t ? x->fraction[i] : 0;
+ c2 = i < MP_T ? x->fraction[i] : 0;
i++;
} while (i__1 < 0 || (i__1 == 0 && c2 < j2));
@@ -787,14 +762,14 @@ mp_divide_integer(const MPNumber *x, int iy, MPNumber *z)
iq -= j1;
}
- iq = iq * MP.b - ir * j2;
+ iq = iq * MP_BASE - ir * j2;
if (iq < 0) {
/* HERE IQ NEGATIVE SO IR WAS TOO LARGE */
ir--;
iq += iy;
}
- if (i < MP.t)
+ if (i < MP_T)
iq += x->fraction[i];
i++;
iqj = iq / iy;
@@ -901,7 +876,7 @@ mp_epowy(const MPNumber *x, MPNumber *z)
/* SEE IF ABS(X) SO LARGE THAT EXP(X) WILL CERTAINLY OVERFLOW
* OR UNDERFLOW. 1.01 IS TO ALLOW FOR ERRORS IN ALOG.
*/
- rlb = log((float) MP.b) * (float)1.01;
+ rlb = log((float) MP_BASE) * (float)1.01;
if (mp_compare_mp_to_float(x, -(double)((float) (MP.m + 1)) * rlb) < 0) {
/* UNDERFLOW SO CALL MPUNFL AND RETURN */
mpunfl(z);
@@ -940,10 +915,10 @@ mp_epowy(const MPNumber *x, MPNumber *z)
/* COMPUTE E-2 OR 1/E USING TWO EXTRA DIGITS IN CASE ABS(X) LARGE
* (BUT ONLY ONE EXTRA DIGIT IF T < 4)
*/
- if (MP.t < 4)
- tss = MP.t + 1;
+ if (MP_T < 4)
+ tss = MP_T + 1;
else
- tss = MP.t + 2;
+ tss = MP_T + 2;
/* LOOP FOR E COMPUTATION. DECREASE T IF POSSIBLE. */
/* Computing MIN */
@@ -1038,19 +1013,19 @@ mpexp1(const MPNumber *x, MPNumber *z)
}
mp_set_from_mp(x, &t1);
- rlb = log((float) MP.b);
+ rlb = log((float) MP_BASE);
/* COMPUTE APPROXIMATELY OPTIMAL Q (AND DIVIDE X BY 2**Q) */
- q = (int) (sqrt((float) MP.t * (float).48 * rlb) + (float) x->exponent *
+ q = (int) (sqrt((float) MP_T * (float).48 * rlb) + (float) x->exponent *
(float)1.44 * rlb);
/* HALVE Q TIMES */
if (q > 0) {
- ib = MP.b << 2;
+ ib = MP_BASE << 2;
ic = 1;
for (i = 1; i <= q; ++i) {
ic <<= 1;
- if (ic < ib && ic != MP.b && i < q)
+ if (ic < ib && ic != MP_BASE && i < q)
continue;
mp_divide_integer(&t1, ic, &t1);
ic = 1;
@@ -1066,7 +1041,7 @@ mpexp1(const MPNumber *x, MPNumber *z)
/* SUM SERIES, REDUCING T WHERE POSSIBLE */
for (i = 2; ; i++) {
- if (MP.t + t2.exponent - z->exponent <= 0)
+ if (MP_T + t2.exponent - z->exponent <= 0)
break;
mp_multiply(&t1, &t2, &t2);
@@ -1203,8 +1178,8 @@ mplns(const MPNumber *x, MPNumber *z)
mp_multiply(x, &t1, z);
/* START NEWTON ITERATION USING SMALL T, LATER INCREASE */
- t = max(5, 13 - (MP.b << 1));
- if (t <= MP.t)
+ t = max(5, 13 - (MP_BASE << 1));
+ if (t <= MP_T)
{
it0 = (t + 5) / 2;
@@ -1217,7 +1192,7 @@ mplns(const MPNumber *x, MPNumber *z)
mp_add(&t3, &t1, &t3);
mp_add(&t2, &t3, &t3);
mp_subtract(z, &t3, z);
- if (t >= MP.t)
+ if (t >= MP_T)
break;
/* FOLLOWING LOOP COMPUTES NEXT VALUE OF T TO USE.
@@ -1225,7 +1200,7 @@ mplns(const MPNumber *x, MPNumber *z)
* WE CAN ALMOST DOUBLE T EACH TIME.
*/
ts3 = t;
- t = MP.t;
+ t = MP_T;
do {
ts2 = t;
t = (t + it0) / 2;
@@ -1234,7 +1209,7 @@ mplns(const MPNumber *x, MPNumber *z)
}
/* CHECK THAT NEWTON ITERATION WAS CONVERGING AS EXPECTED */
- if (t3.sign != 0 && t3.exponent << 1 > it0 - MP.t) {
+ if (t3.sign != 0 && t3.exponent << 1 > it0 - MP_T) {
mperr("*** ERROR OCCURRED IN MPLNS, NEWTON ITERATION NOT CONVERGING PROPERLY ***");
}
}
@@ -1293,7 +1268,7 @@ mp_ln(const MPNumber *x, MPNumber *z)
/* RESTORE EXPONENT AND COMPUTE SINGLE-PRECISION LOG */
t1.exponent = e;
- rlx = log(rx) + (float) e * log((float) MP.b);
+ rlx = log(rx) + (float) e * log((float) MP_BASE);
mp_set_from_float(-(double)rlx, &t2);
/* UPDATE Z AND COMPUTE ACCURATE EXP OF APPROXIMATE LOG */
@@ -1356,7 +1331,7 @@ mp_multiply(const MPNumber *x, const MPNumber *y, MPNumber *z)
int c, i, j, i2, xi;
MPNumber r;
- i2 = MP.t + 4;
+ i2 = MP_T + 4;
memset(&r, 0, sizeof(MPNumber));
@@ -1370,7 +1345,7 @@ mp_multiply(const MPNumber *x, const MPNumber *y, MPNumber *z)
/* PERFORM MULTIPLICATION */
c = 8;
- for (i = 0; i < MP.t; i++) {
+ for (i = 0; i < MP_T; i++) {
xi = x->fraction[i];
/* FOR SPEED, PUT THE NUMBER WITH MANY ZEROS FIRST */
@@ -1378,14 +1353,14 @@ mp_multiply(const MPNumber *x, const MPNumber *y, MPNumber *z)
continue;
/* Computing MIN */
- for (j = 0; j < min(MP.t, i2 - i - 1); j++)
+ for (j = 0; j < min(MP_T, i2 - i - 1); j++)
r.fraction[i+j+1] += xi * y->fraction[j];
c--;
if (c > 0)
continue;
/* CHECK FOR LEGAL BASE B DIGIT */
- if (xi < 0 || xi >= MP.b) {
+ if (xi < 0 || xi >= MP_BASE) {
mperr("*** ILLEGAL BASE B DIGIT IN CALL TO MP_MULTIPLY, POSSIBLE OVERWRITING PROBLEM ***");
mp_set_from_integer(0, z);
return;
@@ -1401,8 +1376,8 @@ mp_multiply(const MPNumber *x, const MPNumber *y, MPNumber *z)
mp_set_from_integer(0, z);
return;
}
- c = ri / MP.b;
- r.fraction[j] = ri - MP.b * c;
+ c = ri / MP_BASE;
+ r.fraction[j] = ri - MP_BASE * c;
}
if (c != 0) {
mperr("*** ILLEGAL BASE B DIGIT IN CALL TO MP_MULTIPLY, POSSIBLE OVERWRITING PROBLEM ***");
@@ -1413,7 +1388,7 @@ mp_multiply(const MPNumber *x, const MPNumber *y, MPNumber *z)
}
if (c != 8) {
- if (xi < 0 || xi >= MP.b) {
+ if (xi < 0 || xi >= MP_BASE) {
mperr("*** ILLEGAL BASE B DIGIT IN CALL TO MP_MULTIPLY, POSSIBLE OVERWRITING PROBLEM ***");
mp_set_from_integer(0, z);
return;
@@ -1427,8 +1402,8 @@ mp_multiply(const MPNumber *x, const MPNumber *y, MPNumber *z)
mp_set_from_integer(0, z);
return;
}
- c = ri / MP.b;
- r.fraction[j] = ri - MP.b * c;
+ c = ri / MP_BASE;
+ r.fraction[j] = ri - MP_BASE * c;
}
if (c != 0) {
@@ -1468,7 +1443,7 @@ mpmul2(const MPNumber *x, int iy, MPNumber *z, int trunc)
z->sign = -z->sign;
/* CHECK FOR MULTIPLICATION BY B */
- if (iy == MP.b) {
+ if (iy == MP_BASE) {
if (x->exponent < MP.m) {
mp_set_from_mp(x, z);
z->sign = z->sign;
@@ -1486,41 +1461,41 @@ mpmul2(const MPNumber *x, int iy, MPNumber *z, int trunc)
/* FORM PRODUCT IN ACCUMULATOR */
c = 0;
- t1 = MP.t + 1;
- t3 = MP.t + 3;
- t4 = MP.t + 4;
+ t1 = MP_T + 1;
+ t3 = MP_T + 3;
+ t4 = MP_T + 4;
/* IF IY*B NOT REPRESENTABLE AS AN INTEGER WE HAVE TO SIMULATE
* DOUBLE-PRECISION MULTIPLICATION.
*/
/* Computing MAX */
- if (iy >= max(MP.b << 3, 32767 / MP.b)) {
+ if (iy >= max(MP_BASE << 3, 32767 / MP_BASE)) {
/* HERE J IS TOO LARGE FOR SINGLE-PRECISION MULTIPLICATION */
- j1 = iy / MP.b;
- j2 = iy - j1 * MP.b;
+ j1 = iy / MP_BASE;
+ j2 = iy - j1 * MP_BASE;
/* FORM PRODUCT */
for (ij = 1; ij <= t4; ++ij) {
- c1 = c / MP.b;
- c2 = c - MP.b * c1;
+ c1 = c / MP_BASE;
+ c2 = c - MP_BASE * c1;
i = t1 - ij;
ix = 0;
if (i > 0)
ix = x->fraction[i - 1];
ri = j2 * ix + c2;
- is = ri / MP.b;
+ is = ri / MP_BASE;
c = j1 * ix + c1 + is;
- z->fraction[i + 3] = ri - MP.b * is;
+ z->fraction[i + 3] = ri - MP_BASE * is;
}
}
else
{
- for (ij = 1; ij <= MP.t; ++ij) {
+ for (ij = 1; ij <= MP_T; ++ij) {
i = t1 - ij;
ri = iy * x->fraction[i - 1] + c;
- c = ri / MP.b;
- z->fraction[i + 3] = ri - MP.b * c;
+ c = ri / MP_BASE;
+ z->fraction[i + 3] = ri - MP_BASE * c;
}
/* CHECK FOR INTEGER OVERFLOW */
@@ -1534,8 +1509,8 @@ mpmul2(const MPNumber *x, int iy, MPNumber *z, int trunc)
for (ij = 1; ij <= 4; ++ij) {
i = 5 - ij;
ri = c;
- c = ri / MP.b;
- z->fraction[i - 1] = ri - MP.b * c;
+ c = ri / MP_BASE;
+ z->fraction[i - 1] = ri - MP_BASE * c;
}
}
@@ -1559,8 +1534,8 @@ mpmul2(const MPNumber *x, int iy, MPNumber *z, int trunc)
z->fraction[i] = z->fraction[i - 1];
}
ri = c;
- c = ri / MP.b;
- z->fraction[0] = ri - MP.b * c;
+ c = ri / MP_BASE;
+ z->fraction[0] = ri - MP_BASE * c;
z->exponent++;
}
}
@@ -1621,7 +1596,7 @@ mp_invert_sign(const MPNumber *x, MPNumber *y)
* AND RETURNS MP RESULT IN Z. INTEGER ARGUMENTS REG_EXP IS
* NOT PRESERVED. R*-ROUNDING IS USED IF TRUNC == 0
*/
-// FIXME: Is r->fraction large enough? It seems to be in practise but it may be MP.t+4 instead of MP.t
+// FIXME: Is r->fraction large enough? It seems to be in practise but it may be MP_T+4 instead of MP_T
// FIXME: There is some sort of stack corruption/use of unitialised variables here. Some functions are
// using stack variables as x otherwise there are corruption errors. e.g. "Cos(45) - 1/Sqrt(2) = -0"
// (try in scientific mode)
@@ -1641,7 +1616,7 @@ mp_normalize(MPNumber *x, int trunc)
return;
}
- i2 = MP.t + 4;
+ i2 = MP_T + 4;
/* Normalize by shifting the fraction to the left */
if (x->fraction[0] == 0) {
@@ -1665,12 +1640,12 @@ mp_normalize(MPNumber *x, int trunc)
/* SEE IF ROUNDING NECESSARY
* TREAT EVEN AND ODD BASES DIFFERENTLY
*/
- b2 = MP.b / 2;
- if (b2 << 1 != MP.b) {
+ b2 = MP_BASE / 2;
+ if (b2 << 1 != MP_BASE) {
round = 0;
/* ODD BASE, ROUND IF R(T+1)... > 1/2 */
for (i = 0; i < 4; i++) {
- i__1 = x->fraction[MP.t + i] - b2;
+ i__1 = x->fraction[MP_T + i] - b2;
if (i__1 < 0)
break;
else if (i__1 == 0)
@@ -1686,12 +1661,12 @@ mp_normalize(MPNumber *x, int trunc)
* AFTER R(T+2).
*/
round = 1;
- i__1 = x->fraction[MP.t] - b2;
+ i__1 = x->fraction[MP_T] - b2;
if (i__1 < 0)
round = 0;
else if (i__1 == 0) {
- if (x->fraction[MP.t - 1] % 2 != 0) {
- if (x->fraction[MP.t + 1] + x->fraction[MP.t + 2] + x->fraction[MP.t + 3] == 0) {
+ if (x->fraction[MP_T - 1] % 2 != 0) {
+ if (x->fraction[MP_T + 1] + x->fraction[MP_T + 2] + x->fraction[MP_T + 3] == 0) {
round = 0;
}
}
@@ -1700,9 +1675,9 @@ mp_normalize(MPNumber *x, int trunc)
/* ROUND */
if (round) {
- for (j = MP.t - 1; j >= 0; j--) {
+ for (j = MP_T - 1; j >= 0; j--) {
++x->fraction[j];
- if (x->fraction[j] < MP.b)
+ if (x->fraction[j] < MP_BASE)
break;
x->fraction[j] = 0;
}
@@ -1858,14 +1833,14 @@ mp_reciprocal(const MPNumber *x, MPNumber *z)
t1.exponent -= ex;
/* START WITH SMALL T TO SAVE TIME. ENSURE THAT B**(T-1) >= 100 */
- if (MP.b < 10)
- t = it[MP.b - 1];
+ if (MP_BASE < 10)
+ t = it[MP_BASE - 1];
else
t = 3;
it0 = (t + 4) / 2;
/* MAIN ITERATION LOOP */
- if (t <= MP.t) {
+ if (t <= MP_T) {
while(1) {
int ts2, ts3;
@@ -1873,23 +1848,23 @@ mp_reciprocal(const MPNumber *x, MPNumber *z)
mp_add_integer(&t2, -1, &t2);
mp_multiply(&t1, &t2, &t2);
mp_subtract(&t1, &t2, &t1);
- if (t >= MP.t)
+ if (t >= MP_T)
break;
/* FOLLOWING LOOP ALMOST DOUBLES T (POSSIBLE
* BECAUSE NEWTONS METHOD HAS 2ND ORDER CONVERGENCE).
*/
ts3 = t;
- t = MP.t;
+ t = MP_T;
do {
ts2 = t;
t = (t + it0) / 2;
} while (t > ts3);
- t = min(ts2, MP.t);
+ t = min(ts2, MP_T);
}
/* RETURN IF NEWTON ITERATION WAS CONVERGING */
- if (t2.sign != 0 && (t1.exponent - t2.exponent) << 1 < MP.t - it0) {
+ if (t2.sign != 0 && (t1.exponent - t2.exponent) << 1 < MP_T - it0) {
/* THE FOLLOWING MESSAGE MAY INDICATE THAT B**(T-1) IS TOO SMALL,
* OR THAT THE STARTING APPROXIMATION IS NOT ACCURATE ENOUGH.
*/
@@ -1901,7 +1876,7 @@ mp_reciprocal(const MPNumber *x, MPNumber *z)
mp_set_from_mp(&t1, z);
/* RESTORE M AND CHECK FOR OVERFLOW (UNDERFLOW IMPOSSIBLE) */
- MP.m += -2;
+ MP.m -= 2;
if (z->exponent > MP.m)
mpovfl(z, "*** OVERFLOW OCCURRED IN MP_RECIPROCAL ***");
}
@@ -1938,7 +1913,7 @@ mp_root(const MPNumber *x, int n, MPNumber *z)
np = abs(n);
/* LOSS OF ACCURACY IF NP LARGE, SO ONLY ALLOW NP <= MAX (B, 64) */
- if (np > max(MP.b, 64)) {
+ if (np > max(MP_BASE, 64)) {
mperr("*** ABS(N) TOO LARGE IN CALL TO MP_ROOT ***");
mp_set_from_integer(0, z);
return;
@@ -1971,7 +1946,7 @@ mp_root(const MPNumber *x, int n, MPNumber *z)
}
/* USE SINGLE-PRECISION ROOT FOR FIRST APPROXIMATION */
- r__1 = exp(((float) (np * ex - x->exponent) * log((float) MP.b) -
+ r__1 = exp(((float) (np * ex - x->exponent) * log((float) MP_BASE) -
log((fabs(rx)))) / (float) np);
mp_set_from_float(r__1, &t1);
@@ -1983,12 +1958,12 @@ mp_root(const MPNumber *x, int n, MPNumber *z)
/* START WITH SMALL T TO SAVE TIME */
/* ENSURE THAT B**(T-1) >= 100 */
- if (MP.b < 10)
- t = it[MP.b - 1];
+ if (MP_BASE < 10)
+ t = it[MP_BASE - 1];
else
t = 3;
- if (t <= MP.t) {
+ if (t <= MP_T) {
/* IT0 IS A NECESSARY SAFETY FACTOR */
it0 = (t + 4) / 2;
@@ -2006,22 +1981,22 @@ mp_root(const MPNumber *x, int n, MPNumber *z)
/* FOLLOWING LOOP ALMOST DOUBLES T (POSSIBLE BECAUSE
* NEWTONS METHOD HAS 2ND ORDER CONVERGENCE).
*/
- if (t >= MP.t)
+ if (t >= MP_T)
break;
ts3 = t;
- t = MP.t;
+ t = MP_T;
do {
ts2 = t;
t = (t + it0) / 2;
} while (t > ts3);
- t = min(ts2, MP.t);
+ t = min(ts2, MP_T);
}
/* NOW R(I2) IS X**(-1/NP)
* CHECK THAT NEWTON ITERATION WAS CONVERGING
*/
- if (t2.sign != 0 && (t1.exponent - t2.exponent) << 1 < MP.t - it0) {
+ if (t2.sign != 0 && (t1.exponent - t2.exponent) << 1 < MP_T - it0) {
/* THE FOLLOWING MESSAGE MAY INDICATE THAT B**(T-1) IS TOO SMALL,
* OR THAT THE INITIAL APPROXIMATION OBTAINED USING ALOG AND EXP
* IS NOT ACCURATE ENOUGH.
@@ -2052,7 +2027,7 @@ mp_sqrt(const MPNumber *x, MPNumber *z)
MPNumber t;
/* MP_ROOT NEEDS 4T+10 WORDS, BUT CAN OVERLAP SLIGHTLY. */
- i2 = MP.t * 3 + 9;
+ i2 = MP_T * 3 + 9;
if (x->sign < 0) {
mperr("*** X NEGATIVE IN CALL TO SUBROUTINE MP_SQRT ***");
} else if (x->sign == 0) {
[
Date Prev][
Date Next] [
Thread Prev][
Thread Next]
[
Thread Index]
[
Date Index]
[
Author Index]
