[gcalctool] General tidyups
- From: Robert Ancell <rancell src gnome org>
- To: svn-commits-list gnome org
- Subject: [gcalctool] General tidyups
- Date: Tue, 19 May 2009 08:53:48 -0400 (EDT)
commit 41b31f5bdc615160b60483408a04d038cd225fb7
Author: Robert Ancell <robert ancell gmail com>
Date: Tue May 19 14:53:40 2009 +0200
General tidyups
---
src/mp.c | 636 ++++++++++++++++++++++++++++++-------------------------------
1 files changed, 313 insertions(+), 323 deletions(-)
diff --git a/src/mp.c b/src/mp.c
index afe5d50..8ea8330 100644
--- a/src/mp.c
+++ b/src/mp.c
@@ -31,15 +31,121 @@
// FIXME: MP.t and MP.m modified inside functions, needs to be fixed to be thread safe
static int mp_compare_mp_to_float(const MPNumber *, float);
-static int pow_ii(int, int);
-static void mpadd2(const MPNumber *, const MPNumber *, MPNumber *, int, int);
-static int mpadd3(const MPNumber *, const MPNumber *, int *, int, int);
-static void mpext(int, int, MPNumber *);
-static void mplns(const MPNumber *, MPNumber *);
-static void mpmaxr(MPNumber *);
-static void mpovfl(MPNumber *, const char *);
-static void mpunfl(MPNumber *);
+/* SETS X TO THE LARGEST POSSIBLE POSITIVE MP NUMBER
+ * CHECK LEGALITY OF B, T, M AND MXR
+ */
+static void
+mpmaxr(MPNumber *x)
+{
+ int i, it;
+
+ mpchk();
+ it = MP.b - 1;
+
+ /* SET FRACTION DIGITS TO B-1 */
+ for (i = 0; i < MP.t; i++)
+ x->fraction[i] = it;
+
+ /* SET SIGN AND EXPONENT */
+ x->sign = 1;
+ x->exponent = MP.m;
+}
+
+
+/* CALLED ON MULTIPLE-PRECISION OVERFLOW, IE WHEN THE
+ * EXPONENT OF MP NUMBER X WOULD EXCEED M.
+ * AT PRESENT EXECUTION IS TERMINATED WITH AN ERROR MESSAGE
+ * AFTER CALLING MPMAXR(X), BUT IT WOULD BE POSSIBLE TO RETURN,
+ * POSSIBLY UPDATING A COUNTER AND TERMINATING EXECUTION AFTER
+ * A PRESET NUMBER OF OVERFLOWS. ACTION COULD EASILY BE DETERMINED
+ * BY A FLAG IN LABELLED COMMON.
+ * M MAY HAVE BEEN OVERWRITTEN, SO CHECK B, T, M ETC.
+ */
+static void
+mpovfl(MPNumber *x, const char *error)
+{
+ fprintf(stderr, "%s\n", error);
+
+ mpchk();
+
+ /* SET X TO LARGEST POSSIBLE POSITIVE NUMBER */
+ mpmaxr(x);
+
+ /* TERMINATE EXECUTION BY CALLING MPERR */
+ mperr("*** CALL TO MPOVFL, MP OVERFLOW OCCURRED ***");
+}
+
+
+/* CALLED ON MULTIPLE-PRECISION UNDERFLOW, IE WHEN THE
+ * EXPONENT OF MP NUMBER X WOULD BE LESS THAN -M.
+ * SINCE M MAY HAVE BEEN OVERWRITTEN, CHECK B, T, M ETC.
+ */
+static void
+mpunfl(MPNumber *x)
+{
+ mpchk();
+
+ /* THE UNDERFLOWING NUMBER IS SET TO ZERO
+ * AN ALTERNATIVE WOULD BE TO CALL MPMINR (X) AND RETURN,
+ * POSSIBLY UPDATING A COUNTER AND TERMINATING EXECUTION
+ * AFTER A PRESET NUMBER OF UNDERFLOWS. ACTION COULD EASILY
+ * BE DETERMINED BY A FLAG IN LABELLED COMMON.
+ */
+ x->sign = 0;
+}
+
+
+static int
+pow_ii(int x, int n)
+{
+ int pow = 1;
+
+ if (n > 0) {
+ for (;;) {
+ if (n & 01) pow *= x;
+ if (n >>= 1) x *= x;
+ else break;
+ }
+ }
+
+ return(pow);
+}
+
+
+/* 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.
+ */
+static void
+mpext(int i, int j, MPNumber *x)
+{
+ int q, s;
+
+ 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];
+
+ /* SET LAST TWO DIGITS TO ZERO */
+ if (s <= q) {
+ x->fraction[MP.t - 2] = 0;
+ x->fraction[MP.t - 1] = 0;
+ return;
+ }
+
+ if (s + q < MP.b * MP.b)
+ return;
+
+ /* ROUND UP HERE */
+ x->fraction[MP.t - 2] = MP.b - 1;
+ x->fraction[MP.t - 1] = MP.b;
+
+ /* NORMALIZE X (LAST DIGIT B IS OK IN MP_MULTIPLY_INTEGER) */
+ mp_multiply_integer(x, 1, x);
+}
/* SETS BASE (B) AND NUMBER OF DIGITS (T) TO GIVE THE
@@ -119,6 +225,22 @@ mp_init(int accuracy)
}
+/* CHECKS LEGALITY OF B, T, M AND MXR.
+ * THE CONDITION ON MXR (THE DIMENSION OF R IN COMMON) IS THAT
+ * MXR >= (I*T + J)
+ */
+// FIXME: MP.t is changed around calls to add/subtract/multiply/divide - it should not be global
+void
+mpchk()
+{
+ /* CHECK LEGALITY OF T AND M */
+ if (MP.t <= 1)
+ mperr("*** T = %d ILLEGAL IN CALL TO MPCHK, PERHAPS NOT SET BEFORE CALL TO AN MP ROUTINE ***", MP.t);
+ if (MP.m <= MP.t)
+ mperr("*** M <= T IN CALL TO MPCHK, PERHAPS NOT SET BEFORE CALL TO AN MP ROUTINE ***");
+}
+
+
/* SETS Y = ABS(X) FOR MP NUMBERS X AND Y */
void
mp_abs(const MPNumber *x, MPNumber *z)
@@ -128,114 +250,11 @@ mp_abs(const MPNumber *x, MPNumber *z)
z->sign = -z->sign;
}
-/* ADDS X AND Y, FORMING RESULT IN Z, WHERE X, Y AND Z ARE MP
- * NUMBERS. FOUR GUARD DIGITS ARE USED, AND THEN R*-ROUNDING.
- */
-void
-mp_add(const MPNumber *x, const MPNumber *y, MPNumber *z)
-{
- mpadd2(x, y, z, y->sign, 0);
-}
-
-/* CALLED BY MP_ADD, MP_SUBTRACT ETC.
- * X, Y AND Z ARE MP NUMBERS, Y_SIGN AND TRUNC ARE INTEGERS.
- * SETS Z = X + Y_SIGN*ABS(Y), WHERE Y_SIGN = +- y->sign.
- * IF TRUNC == 0 R*-ROUNDING IS USED, OTHERWISE TRUNCATION.
- * R*-ROUNDING IS DEFINED IN KUKI AND CODI, COMM. ACM
- * 16(1973), 223. (SEE ALSO BRENT, IEEE TC-22(1973), 601.)
- * CHECK FOR X OR Y ZERO
- */
-static void
-mpadd2(const MPNumber *x, const MPNumber *y, MPNumber *z, int y_sign, int trunc)
-{
- int sign_prod;
- int exp_diff, med;
- int x_largest = 0;
- MPNumber zt; // Use stack variable because of mp_normalize brokeness
-
- /* X = 0 OR NEGLIGIBLE, SO RESULT = +-Y */
- if (x->sign == 0) {
- mp_set_from_mp(y, z);
- z->sign = y_sign;
- return;
- }
-
- /* Y = 0 OR NEGLIGIBLE, SO RESULT = X */
- if (y_sign == 0 || y->sign == 0) {
- mp_set_from_mp(x, z);
- return;
- }
-
- /* COMPARE SIGNS */
- sign_prod = y_sign * x->sign;
- if (abs(sign_prod) > 1) {
- mpchk();
- mperr("*** SIGN NOT 0, +1 OR -1 IN MPADD2 CALL, POSSIBLE OVERWRITING PROBLEM ***");
- z->sign = 0;
- return;
- }
-
- /* COMPARE EXPONENTS */
- exp_diff = x->exponent - y->exponent;
- med = abs(exp_diff);
- if (exp_diff < 0) {
- /* HERE EXPONENT(Y) > EXPONENT(X) */
- 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;
- return;
- }
- x_largest = 0;
- } else if (exp_diff > 0) {
- /* HERE EXPONENT(X) > EXPONENT(Y) */
- 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;
- }
- x_largest = 1;
- } else {
- /* EXPONENTS EQUAL SO COMPARE SIGNS, THEN FRACTIONS IF NEC. */
- if (sign_prod < 0) {
- /* Signs are not equal. find out which mantissa is larger. */
- int j;
- for (j = 0; j < MP.t; j++) {
- int i = x->fraction[j] - y->fraction[j];
- if (i == 0)
- continue;
- if (i < 0)
- x_largest = 0;
- else if (i > 0)
- x_largest = 1;
- break;
- }
-
- /* Both mantissas equal, so result is zero. */
- if (j >= MP.t) {
- z->sign = 0;
- return;
- }
- }
- }
-
- /* NORMALIZE, ROUND OR TRUNCATE, AND RETURN */
- if (x_largest) {
- zt.sign = x->sign;
- zt.exponent = x->exponent + mpadd3(y, x, zt.fraction, sign_prod, med);
- } else {
- zt.sign = y_sign;
- zt.exponent = y->exponent + mpadd3(x, y, zt.fraction, sign_prod, med);
- }
- mp_normalize(&zt, trunc);
- mp_set_from_mp(&zt, z);
-}
-
/* CALLED BY MPADD2, DOES INNER LOOPS OF ADDITION */
/* return value is amount by which exponent needs to be increased. */
static int
-mpadd3(const MPNumber *x, const MPNumber *y, int *r, int s, int med)
+mp_add3(const MPNumber *x, const MPNumber *y, int *r, int s, int med)
{
int i, c;
@@ -340,6 +359,104 @@ mpadd3(const MPNumber *x, const MPNumber *y, int *r, int s, int med)
}
+/* z = x + y_sign * abs(y) */
+static void
+mp_add2(const MPNumber *x, const MPNumber *y, int y_sign, MPNumber *z)
+{
+ int sign_prod;
+ int exp_diff, med;
+ int x_largest = 0;
+ MPNumber zt; // Use stack variable because of mp_normalize brokeness
+
+ /* X = 0 OR NEGLIGIBLE, SO RESULT = +-Y */
+ if (x->sign == 0) {
+ mp_set_from_mp(y, z);
+ z->sign = y_sign;
+ return;
+ }
+
+ /* Y = 0 OR NEGLIGIBLE, SO RESULT = X */
+ if (y->sign == 0 || y->sign == 0) {
+ mp_set_from_mp(x, z);
+ return;
+ }
+
+ /* COMPARE SIGNS */
+ sign_prod = y_sign * x->sign;
+ if (abs(sign_prod) > 1) {
+ mpchk();
+ mperr("*** SIGN NOT 0, +1 OR -1 IN MP_ADD2 CALL, POSSIBLE OVERWRITING PROBLEM ***");
+ z->sign = 0;
+ return;
+ }
+
+ /* COMPARE EXPONENTS */
+ exp_diff = x->exponent - y->exponent;
+ med = abs(exp_diff);
+ if (exp_diff < 0) {
+ /* HERE EXPONENT(Y) > EXPONENT(X) */
+ 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;
+ return;
+ }
+ x_largest = 0;
+ } else if (exp_diff > 0) {
+ /* HERE EXPONENT(X) > EXPONENT(Y) */
+ 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;
+ }
+ x_largest = 1;
+ } else {
+ /* EXPONENTS EQUAL SO COMPARE SIGNS, THEN FRACTIONS IF NEC. */
+ if (sign_prod < 0) {
+ /* Signs are not equal. find out which mantissa is larger. */
+ int j;
+ for (j = 0; j < MP.t; j++) {
+ int i = x->fraction[j] - y->fraction[j];
+ if (i == 0)
+ continue;
+ if (i < 0)
+ x_largest = 0;
+ else if (i > 0)
+ x_largest = 1;
+ break;
+ }
+
+ /* Both mantissas equal, so result is zero. */
+ if (j >= MP.t) {
+ z->sign = 0;
+ return;
+ }
+ }
+ }
+
+ /* NORMALIZE, ROUND OR TRUNCATE, AND RETURN */
+ if (x_largest) {
+ zt.sign = x->sign;
+ zt.exponent = x->exponent + mp_add3(y, x, zt.fraction, sign_prod, med);
+ } else {
+ zt.sign = y_sign;
+ zt.exponent = y->exponent + mp_add3(x, y, zt.fraction, sign_prod, med);
+ }
+ mp_normalize(&zt, 0);
+ mp_set_from_mp(&zt, z);
+}
+
+
+/* ADDS X AND Y, FORMING RESULT IN Z, WHERE X, Y AND Z ARE MP
+ * NUMBERS. FOUR GUARD DIGITS ARE USED, AND THEN R*-ROUNDING.
+ */
+void
+mp_add(const MPNumber *x, const MPNumber *y, MPNumber *z)
+{
+ mp_add2(x, y, y->sign, z);
+}
+
+
/* ADDS MULTIPLE-PRECISION X TO INTEGER Y
* GIVING MULTIPLE-PRECISION Z.
* DIMENSION OF R IN CALLING PROGRAM MUST BE
@@ -372,21 +489,6 @@ mp_add_fraction(const MPNumber *x, int i, int j, MPNumber *y)
}
-/* CHECKS LEGALITY OF B, T, M AND MXR.
- * THE CONDITION ON MXR (THE DIMENSION OF R IN COMMON) IS THAT
- * MXR >= (I*T + J)
- */
-// FIXME: MP.t is changed around calls to add/subtract/multiply/divide - it should not be global
-void
-mpchk()
-{
- /* CHECK LEGALITY OF T AND M */
- if (MP.t <= 1)
- mperr("*** T = %d ILLEGAL IN CALL TO MPCHK, PERHAPS NOT SET BEFORE CALL TO AN MP ROUTINE ***", MP.t);
- if (MP.m <= MP.t)
- mperr("*** M <= T IN CALL TO MPCHK, PERHAPS NOT SET BEFORE CALL TO AN MP ROUTINE ***");
-}
-
/* FOR MP X AND Y, RETURNS FRACTIONAL PART OF X IN Y,
* I.E., Y = X - INTEGER PART OF X (TRUNCATED TOWARDS 0).
*/
@@ -1026,41 +1128,6 @@ mpexp1(const MPNumber *x, MPNumber *y)
}
-/* 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.
- */
-static void
-mpext(int i, int j, MPNumber *x)
-{
- int q, s;
-
- 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];
-
- /* SET LAST TWO DIGITS TO ZERO */
- if (s <= q) {
- x->fraction[MP.t - 2] = 0;
- x->fraction[MP.t - 1] = 0;
- return;
- }
-
- if (s + q < MP.b * MP.b)
- return;
-
- /* ROUND UP HERE */
- x->fraction[MP.t - 2] = MP.b - 1;
- x->fraction[MP.t - 1] = MP.b;
-
- /* NORMALIZE X (LAST DIGIT B IS OK IN MP_MULTIPLY_INTEGER) */
- mp_multiply_integer(x, 1, x);
-}
-
-
/* RETURNS K = K/GCD AND L = L/GCD, WHERE GCD IS THE
* GREATEST COMMON DIVISOR OF K AND L.
* SAVE INPUT PARAMETERS IN LOCAL VARIABLES
@@ -1138,93 +1205,6 @@ mp_is_less_equal(const MPNumber *x, const MPNumber *y)
}
-/* RETURNS Y = LN(X), FOR MP X AND Y, USING MPLNS.
- * RESTRICTION - INTEGER PART OF LN(X) MUST BE REPRESENTABLE
- * AS A SINGLE-PRECISION INTEGER. TIME IS O(SQRT(T).M(T)).
- * FOR SMALL INTEGER X, MPLNI IS FASTER.
- * ASYMPTOTICALLY FASTER METHODS EXIST (EG THE GAUSS-SALAMIN
- * METHOD, SEE MPLNGS), BUT ARE NOT USEFUL UNLESS T IS LARGE.
- * SEE COMMENTS TO MP_ATAN, MPEXP1 AND MPPIGL.
- * DIMENSION OF R IN CALLING PROGRAM MUST BE AT LEAST 6T+14.
- * CHECK LEGALITY OF B, T, M AND MXR
- */
-void
-mp_ln(const MPNumber *x, MPNumber *z)
-{
- int e, k;
- float rx, rlx;
- MPNumber t1, t2;
-
- mpchk();
-
- /* CHECK THAT X IS POSITIVE */
- if (x->sign <= 0) {
- mperr("*** X NONPOSITIVE IN CALL TO MP_LN ***");
- z->sign = 0;
- return;
- }
-
- /* MOVE X TO LOCAL STORAGE */
- mp_set_from_mp(x, &t1);
- z->sign = 0;
- k = 0;
-
- /* LOOP TO GET APPROXIMATE LN(X) USING SINGLE-PRECISION */
- while(1)
- {
- mp_add_integer(&t1, -1, &t2);
-
- /* IF POSSIBLE GO TO CALL MPLNS */
- if (t2.sign == 0 || t2.exponent + 1 <= 0) {
- /* COMPUTE FINAL CORRECTION ACCURATELY USING MPLNS */
- mplns(&t2, &t2);
- mp_add(z, &t2, z);
- return;
- }
-
- /* REMOVE EXPONENT TO AVOID FLOATING-POINT OVERFLOW */
- e = t1.exponent;
- t1.exponent = 0;
- rx = mp_cast_to_float(&t1);
-
- /* RESTORE EXPONENT AND COMPUTE SINGLE-PRECISION LOG */
- t1.exponent = e;
- rlx = log(rx) + (float) e * log((float) MP.b);
- mp_set_from_float(-(double)rlx, &t2);
-
- /* UPDATE Z AND COMPUTE ACCURATE EXP OF APPROXIMATE LOG */
- mp_subtract(z, &t2, z);
- mp_epowy(&t2, &t2);
-
- /* COMPUTE RESIDUAL WHOSE LOG IS STILL TO BE FOUND */
- mp_multiply(&t1, &t2, &t1);
-
- /* MAKE SURE NOT LOOPING INDEFINITELY */
- ++k;
- if (k >= 10) {
- mperr("*** ERROR IN MP_LN, ITERATION NOT CONVERGING ***");
- return;
- }
- }
-}
-
-
-/* MP precision common log.
- *
- * 1. log10(x) = log10(e) * log(x)
- */
-void
-mp_logarithm(int n, MPNumber *x, MPNumber *z)
-{
- MPNumber t1, t2;
-
- mp_set_from_integer(n, &t1);
- mp_ln(&t1, &t1);
- mp_ln(x, &t2);
- mp_divide(&t2, &t1, z);
-}
-
-
/* RETURNS MP Y = LN(1+X) IF X IS AN MP NUMBER SATISFYING THE
* CONDITION ABS(X) < 1/B, ERROR OTHERWISE.
* USES NEWTONS METHOD TO SOLVE THE EQUATION
@@ -1309,32 +1289,98 @@ mplns(const MPNumber *x, MPNumber *y)
}
-/* RETURNS LOGICAL VALUE OF (X < Y) FOR MP X AND Y. */
-int
-mp_is_less_than(const MPNumber *x, const MPNumber *y)
+/* RETURNS Y = LN(X), FOR MP X AND Y, USING MPLNS.
+ * RESTRICTION - INTEGER PART OF LN(X) MUST BE REPRESENTABLE
+ * AS A SINGLE-PRECISION INTEGER. TIME IS O(SQRT(T).M(T)).
+ * FOR SMALL INTEGER X, MPLNI IS FASTER.
+ * ASYMPTOTICALLY FASTER METHODS EXIST (EG THE GAUSS-SALAMIN
+ * METHOD, SEE MPLNGS), BUT ARE NOT USEFUL UNLESS T IS LARGE.
+ * SEE COMMENTS TO MP_ATAN, MPEXP1 AND MPPIGL.
+ * DIMENSION OF R IN CALLING PROGRAM MUST BE AT LEAST 6T+14.
+ * CHECK LEGALITY OF B, T, M AND MXR
+ */
+void
+mp_ln(const MPNumber *x, MPNumber *z)
{
- return (mp_compare_mp_to_mp(x, y) < 0);
+ int e, k;
+ float rx, rlx;
+ MPNumber t1, t2;
+
+ mpchk();
+
+ /* CHECK THAT X IS POSITIVE */
+ if (x->sign <= 0) {
+ mperr("*** X NONPOSITIVE IN CALL TO MP_LN ***");
+ z->sign = 0;
+ return;
+ }
+
+ /* MOVE X TO LOCAL STORAGE */
+ mp_set_from_mp(x, &t1);
+ z->sign = 0;
+ k = 0;
+
+ /* LOOP TO GET APPROXIMATE LN(X) USING SINGLE-PRECISION */
+ while(1)
+ {
+ mp_add_integer(&t1, -1, &t2);
+
+ /* IF POSSIBLE GO TO CALL MPLNS */
+ if (t2.sign == 0 || t2.exponent + 1 <= 0) {
+ /* COMPUTE FINAL CORRECTION ACCURATELY USING MPLNS */
+ mplns(&t2, &t2);
+ mp_add(z, &t2, z);
+ return;
+ }
+
+ /* REMOVE EXPONENT TO AVOID FLOATING-POINT OVERFLOW */
+ e = t1.exponent;
+ t1.exponent = 0;
+ rx = mp_cast_to_float(&t1);
+
+ /* RESTORE EXPONENT AND COMPUTE SINGLE-PRECISION LOG */
+ t1.exponent = e;
+ rlx = log(rx) + (float) e * log((float) MP.b);
+ mp_set_from_float(-(double)rlx, &t2);
+
+ /* UPDATE Z AND COMPUTE ACCURATE EXP OF APPROXIMATE LOG */
+ mp_subtract(z, &t2, z);
+ mp_epowy(&t2, &t2);
+
+ /* COMPUTE RESIDUAL WHOSE LOG IS STILL TO BE FOUND */
+ mp_multiply(&t1, &t2, &t1);
+
+ /* MAKE SURE NOT LOOPING INDEFINITELY */
+ ++k;
+ if (k >= 10) {
+ mperr("*** ERROR IN MP_LN, ITERATION NOT CONVERGING ***");
+ return;
+ }
+ }
}
-/* SETS X TO THE LARGEST POSSIBLE POSITIVE MP NUMBER
- * CHECK LEGALITY OF B, T, M AND MXR
+/* MP precision common log.
+ *
+ * 1. log10(x) = log10(e) * log(x)
*/
-static void
-mpmaxr(MPNumber *x)
+void
+mp_logarithm(int n, MPNumber *x, MPNumber *z)
{
- int i, it;
+ MPNumber t1, t2;
- mpchk();
- it = MP.b - 1;
+ mp_set_from_integer(n, &t1);
+ mp_ln(&t1, &t1);
+ mp_ln(x, &t2);
+ mp_divide(&t2, &t1, z);
+}
- /* SET FRACTION DIGITS TO B-1 */
- for (i = 0; i < MP.t; i++)
- x->fraction[i] = it;
- /* SET SIGN AND EXPONENT */
- x->sign = 1;
- x->exponent = MP.m;
+/* RETURNS LOGICAL VALUE OF (X < Y) FOR MP X AND Y. */
+int
+mp_is_less_than(const MPNumber *x, const MPNumber *y)
+{
+ return (mp_compare_mp_to_mp(x, y) < 0);
}
@@ -1735,29 +1781,6 @@ mp_normalize(MPNumber *x, int trunc)
}
-/* CALLED ON MULTIPLE-PRECISION OVERFLOW, IE WHEN THE
- * EXPONENT OF MP NUMBER X WOULD EXCEED M.
- * AT PRESENT EXECUTION IS TERMINATED WITH AN ERROR MESSAGE
- * AFTER CALLING MPMAXR(X), BUT IT WOULD BE POSSIBLE TO RETURN,
- * POSSIBLY UPDATING A COUNTER AND TERMINATING EXECUTION AFTER
- * A PRESET NUMBER OF OVERFLOWS. ACTION COULD EASILY BE DETERMINED
- * BY A FLAG IN LABELLED COMMON.
- * M MAY HAVE BEEN OVERWRITTEN, SO CHECK B, T, M ETC.
- */
-static void
-mpovfl(MPNumber *x, const char *error)
-{
- fprintf(stderr, "%s\n", error);
-
- mpchk();
-
- /* SET X TO LARGEST POSSIBLE POSITIVE NUMBER */
- mpmaxr(x);
-
- /* TERMINATE EXECUTION BY CALLING MPERR */
- mperr("*** CALL TO MPOVFL, MP OVERFLOW OCCURRED ***");
-}
-
/* 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).
@@ -2146,42 +2169,9 @@ mp_sqrt(const MPNumber *x, MPNumber *z)
void
mp_subtract(const MPNumber *x, const MPNumber *y, MPNumber *z)
{
- mpadd2(x, y, z, -y->sign, 0);
+ mp_add2(x, y, -y->sign, z);
}
-/* CALLED ON MULTIPLE-PRECISION UNDERFLOW, IE WHEN THE
- * EXPONENT OF MP NUMBER X WOULD BE LESS THAN -M.
- * SINCE M MAY HAVE BEEN OVERWRITTEN, CHECK B, T, M ETC.
- */
-static void
-mpunfl(MPNumber *x)
-{
- mpchk();
-
- /* THE UNDERFLOWING NUMBER IS SET TO ZERO
- * AN ALTERNATIVE WOULD BE TO CALL MPMINR (X) AND RETURN,
- * POSSIBLY UPDATING A COUNTER AND TERMINATING EXECUTION
- * AFTER A PRESET NUMBER OF UNDERFLOWS. ACTION COULD EASILY
- * BE DETERMINED BY A FLAG IN LABELLED COMMON.
- */
- x->sign = 0;
-}
-
-static int
-pow_ii(int x, int n)
-{
- int pow = 1;
-
- if (n > 0) {
- for (;;) {
- if (n & 01) pow *= x;
- if (n >>= 1) x *= x;
- else break;
- }
- }
-
- return(pow);
-}
void
mp_factorial(const MPNumber *x, MPNumber *z)
[
Date Prev][
Date Next] [
Thread Prev][
Thread Next]
[
Thread Index]
[
Date Index]
[
Author Index]
