[gcalctool/gcalctool-newui2] Tridy up MP
- From: Robert Ancell <rancell src gnome org>
- To: svn-commits-list gnome org
- Subject: [gcalctool/gcalctool-newui2] Tridy up MP
- Date: Sun, 12 Jul 2009 04:57:11 +0000 (UTC)
commit 927a9f9b0b5603818a77cd60699e6b5ba26744b5
Author: Robert Ancell <robert ancell gmail com>
Date: Sat Jul 11 01:25:23 2009 +1000
Tridy up MP
src/mp-convert.c | 48 ++++------
src/mp-trigonometric.c | 50 +++-------
src/mp.c | 237 +++++++++++++++++-------------------------------
3 files changed, 120 insertions(+), 215 deletions(-)
---
diff --git a/src/mp-convert.c b/src/mp-convert.c
index 0dea37f..2dedc1f 100644
--- a/src/mp-convert.c
+++ b/src/mp-convert.c
@@ -31,17 +31,10 @@
* SEE IF X AND Y HAVE THE SAME ADDRESS (THEY OFTEN DO)
*/
void
-mp_set_from_mp(const MPNumber *x, MPNumber *y)
+mp_set_from_mp(const MPNumber *x, MPNumber *z)
{
- /* HERE X AND Y MUST HAVE THE SAME ADDRESS */
- if (x == y)
- return;
-
- /* NO NEED TO COPY X[1],X[2],... IF X[0] == 0 */
- if (x->sign == 0)
- y->sign = 0;
- else
- memcpy (y, x, (MP.t + 2) * sizeof(int));
+ if (z != x)
+ memcpy(z, x, sizeof(MPNumber));
}
/* CONVERTS SINGLE-PRECISION NUMBER RX TO MULTIPLE-PRECISION Z.
@@ -55,8 +48,9 @@ mp_set_from_float(float rx, MPNumber *z)
int i, k, i2, ib, ie, tp;
float rb, rj;
- mpchk();
i2 = MP.t + 4;
+
+ memset(z, 0, sizeof(MPNumber));
/* CHECK SIGN */
if (rx < (float) 0.0) {
@@ -67,7 +61,7 @@ mp_set_from_float(float rx, MPNumber *z)
rj = rx;
} else {
/* IF RX = 0E0 RETURN 0 */
- z->sign = 0;
+ mp_set_from_integer(0, z);
return;
}
@@ -146,9 +140,10 @@ mp_set_from_double(double dx, MPNumber *z)
int i, k, i2, ib, ie, tp;
double db, dj;
- mpchk();
i2 = MP.t + 4;
+ memset(z, 0, sizeof(MPNumber));
+
/* CHECK SIGN */
if (dx < 0.) {
z->sign = -1;
@@ -157,7 +152,7 @@ mp_set_from_double(double dx, MPNumber *z)
z->sign = 1;
dj = dx;
} else {
- z->sign = 0;
+ mp_set_from_integer(0, z);
return;
}
@@ -219,12 +214,10 @@ mp_set_from_double(double dx, MPNumber *z)
void
mp_set_from_integer(int ix, MPNumber *z)
{
- mpchk();
+ memset(z, 0, sizeof(MPNumber));
- if (ix == 0) {
- z->sign = 0;
+ if (ix == 0)
return;
- }
if (ix < 0) {
ix = -ix;
@@ -248,23 +241,24 @@ mp_set_from_integer(int ix, MPNumber *z)
/* CONVERTS THE RATIONAL NUMBER I/J TO MULTIPLE PRECISION Q. */
void
-mp_set_from_fraction(int i, int j, MPNumber *q)
+mp_set_from_fraction(int i, int j, MPNumber *z)
{
mpgcd(&i, &j);
if (j == 0) {
- mperr("*** J == 0 IN CALL TO MP_SET_FROM_FRACTION ***\n");
- q->sign = 0;
- return;
+ mperr("*** J == 0 IN CALL TO MP_SET_FROM_FRACTION ***\n");
+ mp_set_from_integer(0, z);
+ return;
}
if (j < 0) {
- i = -i;
- j = -j;
+ i = -i;
+ j = -j;
}
- mp_set_from_integer(i, q);
- if (j != 1) mp_divide_integer(q, j, q);
+ mp_set_from_integer(i, z);
+ if (j != 1)
+ mp_divide_integer(z, j, z);
}
/* CONVERTS MULTIPLE-PRECISION X TO INTEGER, AND
@@ -366,7 +360,6 @@ mp_cast_to_float(const MPNumber *x)
int i;
float rb, rz = 0.0;
- mpchk();
if (x->sign == 0)
return 0.0;
@@ -433,7 +426,6 @@ mp_cast_to_double(const MPNumber *x)
int i, tm = 0;
double d__1, db, dz2, ret_val = 0.0;
- mpchk();
if (x->sign == 0)
return 0.0;
diff --git a/src/mp-trigonometric.c b/src/mp-trigonometric.c
index 698d60f..73c2ba2 100644
--- a/src/mp-trigonometric.c
+++ b/src/mp-trigonometric.c
@@ -39,8 +39,6 @@ mp_compare_mp_to_int(const MPNumber *x, int i)
{
MPNumber t;
- mpchk();
-
/* CONVERT I TO MULTIPLE-PRECISION AND COMPARE */
mp_set_from_integer(i, &t);
return mp_compare_mp_to_mp(x, &t);
@@ -63,13 +61,12 @@ mpsin1(const MPNumber *x, MPNumber *z, int do_sin)
int i, b2;
MPNumber t1, t2;
- mpchk();
-
/* SIN(0) = 0, COS(0) = 1 */
if (x->sign == 0) {
- z->sign = 0;
if (do_sin == 0)
mp_set_from_integer(1, z);
+ else
+ mp_set_from_integer(0, z);
return;
}
@@ -84,12 +81,13 @@ mpsin1(const MPNumber *x, MPNumber *z, int do_sin)
else
mp_set_from_mp(x, &t1);
- z->sign = 0;
i = 1;
if (do_sin != 0) {
mp_set_from_mp(&t1, z);
i = 2;
}
+ else
+ mp_set_from_integer(0, z);
/* POWER SERIES LOOP. REDUCE T IF POSSIBLE */
for (; ; i+= 2) {
@@ -135,10 +133,9 @@ 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;
+ mp_set_from_integer(0, z);
return;
}
@@ -250,8 +247,6 @@ mp_get_pi(MPNumber *z)
float prec;
MPNumber t;
- mpchk();
-
mp_atan1N(5, &t);
mp_multiply_integer(&t, 4, &t);
mp_atan1N(239, z);
@@ -278,12 +273,10 @@ mp_sin(const MPNumber *x, MPAngleUnit unit, MPNumber *z)
float rx = 0.0;
MPNumber t1, t2;
- mpchk();
-
convert_to_radians(x, unit, &t1);
if (t1.sign == 0) {
- z->sign = 0;
+ mp_set_from_integer(0, z);
return;
}
@@ -316,7 +309,7 @@ mp_sin(const MPNumber *x, MPAngleUnit unit, MPNumber *z)
mp_add_fraction(&t1, -1, 2, &t1);
xs = -xs * t1.sign;
if (xs == 0) {
- z->sign = 0;
+ mp_set_from_integer(0, z);
return;
}
@@ -328,7 +321,7 @@ mp_sin(const MPNumber *x, MPAngleUnit unit, MPNumber *z)
mp_add_integer(&t1, -2, &t1);
if (t1.sign == 0) {
- z->sign = 0;
+ mp_set_from_integer(0, z);
return;
}
@@ -376,9 +369,6 @@ mp_cos(const MPNumber *xi, MPAngleUnit unit, MPNumber *z)
return;
}
- /* CHECK LEGALITY OF B, T, M AND MXR */
- mpchk();
-
convert_to_radians(xi, unit, z);
/* SEE IF ABS(X) <= 1 */
@@ -427,9 +417,8 @@ mp_asin(const MPNumber *x, MPAngleUnit unit, MPNumber *z)
{
MPNumber t1, t2;
- mpchk();
if (x->sign == 0) {
- z->sign = 0;
+ mp_set_from_integer(0, z);
return;
}
@@ -485,11 +474,11 @@ mp_acos(const MPNumber *x, MPAngleUnit unit, MPNumber *z)
if (mp_is_greater_than(x, &MP1) || mp_is_less_than(x, &MPn1)) {
mperr("Error");
- z->sign = 0;
+ mp_set_from_integer(0, z);
} else if (x->sign == 0) {
mp_divide_integer(&MPpi, 2, z);
} else if (mp_is_equal(x, &MP1)) {
- z->sign = 0;
+ mp_set_from_integer(0, z);
} else if (mp_is_equal(x, &MPn1)) {
mp_set_from_mp(&MPpi, z);
} else {
@@ -526,9 +515,8 @@ mp_atan(const MPNumber *x, MPAngleUnit unit, MPNumber *z)
float rx = 0.0;
MPNumber t1, t2;
- mpchk();
if (x->sign == 0) {
- z->sign = 0;
+ mp_set_from_integer(0, z);
return;
}
@@ -651,7 +639,7 @@ mp_atanh(const MPNumber *x, MPNumber *z)
if (mp_is_greater_equal(x, &MP1) || mp_is_less_equal(x, &MPn1)) {
mperr("Error");
- z->sign = 0;
+ mp_set_from_integer(0, z);
} else {
mp_add(&MP1, x, &MP2);
mp_subtract(&MP1, x, &MP3);
@@ -677,8 +665,6 @@ mp_cosh(const MPNumber *x, MPNumber *z)
return;
}
- /* CHECK LEGALITY OF B, T, M AND MXR */
- mpchk();
mp_abs(x, &t);
/* IF ABS(X) TOO LARGE MP_EPOWY WILL PRINT ERROR MESSAGE
@@ -709,13 +695,10 @@ mp_sinh(const MPNumber *x, MPNumber *z)
xs = x->sign;
if (xs == 0) {
- z->sign = 0;
+ mp_set_from_integer(0, z);
return;
}
- /* CHECK LEGALITY OF B, T, M AND MXR */
- mpchk();
-
/* WORK WITH ABS(X) */
mp_abs(x, &t2);
@@ -759,13 +742,10 @@ mp_tanh(const MPNumber *x, MPNumber *z)
/* TANH(0) = 0 */
if (x->sign == 0) {
- z->sign = 0;
+ mp_set_from_integer(0, z);
return;
}
- /* CHECK LEGALITY OF B, T, M AND MXR */
- mpchk();
-
/* SAVE SIGN AND WORK WITH ABS(X) */
xs = x->sign;
mp_abs(x, &t);
diff --git a/src/mp.c b/src/mp.c
index c5e57a7..1755fe1 100644
--- a/src/mp.c
+++ b/src/mp.c
@@ -38,7 +38,6 @@ mpmaxr(MPNumber *x)
{
int i, it;
- mpchk();
it = MP.b - 1;
/* SET FRACTION DIGITS TO B-1 */
@@ -65,8 +64,6 @@ mpovfl(MPNumber *x, const char *error)
{
fprintf(stderr, "%s\n", error);
- mpchk();
-
/* SET X TO LARGEST POSSIBLE POSITIVE NUMBER */
mpmaxr(x);
@@ -80,17 +77,15 @@ mpovfl(MPNumber *x, const char *error)
* SINCE M MAY HAVE BEEN OVERWRITTEN, CHECK B, T, M ETC.
*/
static void
-mpunfl(MPNumber *x)
+mpunfl(MPNumber *z)
{
- 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;
+ mp_set_from_integer(0, z);
}
@@ -217,25 +212,6 @@ mp_init(int accuracy)
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;
}
-
- /* CHECK LEGALITY OF B, T, M AND MXR (AT LEAST T+4) */
- mpchk();
-}
-
-
-/* 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 ***");
}
@@ -374,6 +350,8 @@ mp_add2(const MPNumber *x, const MPNumber *y, int y_sign, MPNumber *z)
int exp_diff, med;
int x_largest = 0;
MPNumber zt; // Use stack variable because of mp_normalize brokeness
+
+ memset(&zt, 0, sizeof(MPNumber));
/* X = 0 OR NEGLIGIBLE, SO RESULT = +-Y */
if (x->sign == 0) {
@@ -391,9 +369,8 @@ mp_add2(const MPNumber *x, const MPNumber *y, int y_sign, MPNumber *z)
/* 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;
+ mp_set_from_integer(0, z);
return;
}
@@ -435,7 +412,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) {
- z->sign = 0;
+ mp_set_from_integer(0, z);
return;
}
}
@@ -450,7 +427,7 @@ mp_add2(const MPNumber *x, const MPNumber *y, int y_sign, MPNumber *z)
zt.exponent = y->exponent + mp_add3(x, y, zt.fraction, sign_prod, med);
}
mp_normalize(&zt, 0);
- mp_set_from_mp(&zt, z);
+ mp_set_from_mp(&zt, z);
}
@@ -476,7 +453,6 @@ void
mp_add_integer(const MPNumber *x, int y, MPNumber *z)
{
MPNumber t;
- mpchk();
mp_set_from_integer(y, &t);
mp_add(x, &t, z);
}
@@ -490,7 +466,6 @@ void
mp_add_fraction(const MPNumber *x, int i, int j, MPNumber *y)
{
MPNumber t;
- mpchk();
mp_set_from_fraction(i, j, &t);
mp_add(x, &t, y);
}
@@ -500,32 +475,32 @@ mp_add_fraction(const MPNumber *x, int i, int j, MPNumber *y)
* I.E., Y = X - INTEGER PART OF X (TRUNCATED TOWARDS 0).
*/
void
-mp_fractional_component(const MPNumber *x, MPNumber *y)
+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) {
- y->sign = 0;
+ mp_set_from_integer(0, z);
return;
}
/* IF EXPONENT NOT POSITIVE CAN RETURN X */
if (x->exponent <= 0) {
- mp_set_from_mp(x, y);
+ mp_set_from_mp(x, z);
return;
}
/* MOVE FRACTIONAL PART OF X TO ACCUMULATOR */
- y->sign = x->sign;
- y->exponent = x->exponent;
- memset(y->fraction, 0, y->exponent*sizeof(int));
- if (x != y)
- memcpy(y->fraction + y->exponent, x->fraction + x->exponent,
+ z->sign = x->sign;
+ z->exponent = x->exponent;
+ 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(y->fraction + MP.t, 0, 4*sizeof(int));
+ memset(z->fraction + MP.t, 0, 4*sizeof(int));
/* NORMALIZE RESULT AND RETURN */
- mp_normalize(y, 1);
+ mp_normalize(z, 1);
}
/* RETURNS Y = INTEGER PART OF X (TRUNCATED TOWARDS 0), FOR MP X AND Y.
@@ -534,29 +509,28 @@ mp_fractional_component(const MPNumber *x, MPNumber *y)
* CHECK LEGALITY OF B, T, M AND MXR
*/
void
-mp_integer_component(const MPNumber *x, MPNumber *y)
+mp_integer_component(const MPNumber *x, MPNumber *z)
{
int i;
- mpchk();
- mp_set_from_mp(x, y);
- if (y->sign == 0)
+ mp_set_from_mp(x, z);
+ if (z->sign == 0)
return;
/* IF EXPONENT LARGE ENOUGH RETURN Y = X */
- if (y->exponent >= MP.t) {
+ if (z->exponent >= MP.t) {
return;
}
/* IF EXPONENT SMALL ENOUGH RETURN ZERO */
- if (y->exponent <= 0) {
- y->sign = 0;
+ if (z->exponent <= 0) {
+ mp_set_from_integer(0, z);
return;
}
/* SET FRACTION TO ZERO */
- for (i = y->exponent; i < MP.t; i++) {
- y->fraction[i] = 0;
+ for (i = z->exponent; i < MP.t; i++) {
+ z->fraction[i] = 0;
}
}
@@ -571,7 +545,6 @@ static int
mp_compare_mp_to_float(const MPNumber *x, float ri)
{
MPNumber t;
- mpchk();
/* CONVERT RI TO MULTIPLE-PRECISION AND COMPARE */
mp_set_from_float(ri, &t);
@@ -634,18 +607,16 @@ mp_divide(const MPNumber *x, const MPNumber *y, MPNumber *z)
int i, ie, iz3;
MPNumber t;
- mpchk();
-
/* CHECK FOR DIVISION BY ZERO */
if (y->sign == 0) {
mperr("*** ATTEMPTED DIVISION BY ZERO IN CALL TO MP_DIVIDE ***");
- z->sign = 0;
+ mp_set_from_integer(0, z);
return;
}
/* CHECK FOR X = 0 */
if (x->sign == 0) {
- z->sign = 0;
+ mp_set_from_integer(0, z);
return;
}
@@ -689,7 +660,7 @@ mp_divide_integer(const MPNumber *x, int iy, MPNumber *z)
if (iy == 0) {
mperr("*** ATTEMPTED DIVISION BY ZERO IN CALL TO MP_DIVIDE_INTEGER ***");
- z->sign = 0;
+ mp_set_from_integer(0, z);
return;
}
@@ -834,9 +805,8 @@ mp_divide_integer(const MPNumber *x, int iy, MPNumber *z)
L210:
/* CARRY NEGATIVE SO OVERFLOW MUST HAVE OCCURRED */
- mpchk();
mperr("*** INTEGER OVERFLOW IN MP_DIVIDE_INTEGER, B TOO LARGE ***");
- z->sign = 0;
+ mp_set_from_integer(0, z);
}
@@ -904,8 +874,6 @@ mp_epowy(const MPNumber *x, MPNumber *z)
float rx, rz, rlb;
MPNumber t1, t2;
- mpchk();
-
/* CHECK FOR X == 0 */
if (x->sign == 0) {
mp_set_from_integer(1, z);
@@ -969,10 +937,7 @@ mp_epowy(const MPNumber *x, MPNumber *z)
/* LOOP FOR E COMPUTATION. DECREASE T IF POSSIBLE. */
/* Computing MIN */
- ts = MP.t;
- MP.t = tss;
mp_set_from_integer(xs, &t1);
- MP.t = ts;
t2.sign = 0;
for (i = 2 ; ; i++) {
@@ -1056,24 +1021,22 @@ mp_epowy(const MPNumber *x, MPNumber *z)
* CHECK LEGALITY OF B, T, M AND MXR
*/
void
-mpexp1(const MPNumber *x, MPNumber *y)
+mpexp1(const MPNumber *x, MPNumber *z)
{
int i, q, ib, ic;
float rlb;
MPNumber t1, t2;
- mpchk();
-
/* CHECK FOR X = 0 */
if (x->sign == 0) {
- y->sign = 0;
+ mp_set_from_integer(0, z);
return;
}
/* CHECK THAT ABS(X) < 1 */
if (x->exponent > 0) {
mperr("*** ABS(X) NOT LESS THAN 1 IN CALL TO MPEXP1 ***");
- y->sign = 0;
+ mp_set_from_integer(0, z);
return;
}
@@ -1098,17 +1061,17 @@ mpexp1(const MPNumber *x, MPNumber *y)
}
if (t1.sign == 0) {
- y->sign = 0;
+ mp_set_from_integer(0, z);
return;
}
- mp_set_from_mp(&t1, y);
+ mp_set_from_mp(&t1, z);
mp_set_from_mp(&t1, &t2);
/* SUM SERIES, REDUCING T WHERE POSSIBLE */
for (i = 2; ; i++) {
int t, ts;
- t = MP.t + 2 + t2.exponent - y->exponent;
+ t = MP.t + 2 + t2.exponent - z->exponent;
if (t <= 2)
break;
t = min(t, MP.t);
@@ -1119,7 +1082,7 @@ mpexp1(const MPNumber *x, MPNumber *y)
mp_divide_integer(&t2, i, &t2);
MP.t = ts;
- mp_add(&t2, y, y);
+ mp_add(&t2, z, z);
if (t2.sign == 0)
break;
}
@@ -1129,8 +1092,8 @@ mpexp1(const MPNumber *x, MPNumber *y)
/* APPLY (X+1)**2 - 1 = X(2 + X) FOR Q ITERATIONS */
for (i = 1; i <= q; ++i) {
- mp_add_integer(y, 2, &t1);
- mp_multiply(&t1, y, y);
+ mp_add_integer(z, 2, &t1);
+ mp_multiply(&t1, z, z);
}
}
@@ -1221,23 +1184,21 @@ mp_is_less_equal(const MPNumber *x, const MPNumber *y)
* CHECK LEGALITY OF B, T, M AND MXR
*/
static void
-mplns(const MPNumber *x, MPNumber *y)
+mplns(const MPNumber *x, MPNumber *z)
{
int t, it0;
MPNumber t1, t2, t3;
- mpchk();
-
/* CHECK FOR X = 0 EXACTLY */
if (x->sign == 0) {
- y->sign = 0;
+ mp_set_from_integer(0, z);
return;
}
/* CHECK THAT ABS(X) < 1/B */
if (x->exponent + 1 > 0) {
mperr("*** ABS(X) >= 1/B IN CALL TO MPLNS ***");
- y->sign = 0;
+ mp_set_from_integer(0, z);
return;
}
@@ -1249,7 +1210,7 @@ mplns(const MPNumber *x, MPNumber *y)
mp_add_fraction(&t1, 1, 2, &t1);
mp_multiply(x, &t1, &t1);
mp_add_integer(&t1, -1, &t1);
- mp_multiply(x, &t1, y);
+ mp_multiply(x, &t1, z);
/* START NEWTON ITERATION USING SMALL T, LATER INCREASE */
t = max(5, 13 - (MP.b << 1));
@@ -1263,11 +1224,11 @@ mplns(const MPNumber *x, MPNumber *y)
ts = MP.t;
MP.t = t;
- mpexp1(y, &t3);
+ mpexp1(z, &t3);
mp_multiply(&t2, &t3, &t1);
mp_add(&t3, &t1, &t3);
mp_add(&t2, &t3, &t3);
- mp_subtract(y, &t3, y);
+ mp_subtract(z, &t3, z);
MP.t = ts;
if (t >= MP.t)
break;
@@ -1292,7 +1253,7 @@ mplns(const MPNumber *x, MPNumber *y)
}
/* REVERSE SIGN OF Y AND RETURN */
- y->sign = -y->sign;
+ z->sign = -z->sign;
}
@@ -1313,18 +1274,16 @@ mp_ln(const MPNumber *x, MPNumber *z)
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;
+ mp_set_from_integer(0, z);
return;
}
/* MOVE X TO LOCAL STORAGE */
mp_set_from_mp(x, &t1);
- z->sign = 0;
+ mp_set_from_integer(0, z);
k = 0;
/* LOOP TO GET APPROXIMATE LN(X) USING SINGLE-PRECISION */
@@ -1410,9 +1369,10 @@ mp_multiply(const MPNumber *x, const MPNumber *y, MPNumber *z)
int c, i, j, i2, xi;
MPNumber r;
- mpchk();
i2 = MP.t + 4;
-
+
+ memset(&r, 0, sizeof(MPNumber));
+
/* FORM SIGN OF PRODUCT */
z->sign = x->sign * y->sign;
if (z->sign == 0)
@@ -1421,10 +1381,6 @@ mp_multiply(const MPNumber *x, const MPNumber *y, MPNumber *z)
/* FORM EXPONENT OF PRODUCT */
z->exponent = x->exponent + y->exponent;
- /* CLEAR ACCUMULATOR */
- for (i = 0; i < i2; i++)
- r.fraction[i] = 0;
-
/* PERFORM MULTIPLICATION */
c = 8;
for (i = 0; i < MP.t; i++) {
@@ -1444,7 +1400,7 @@ mp_multiply(const MPNumber *x, const MPNumber *y, MPNumber *z)
/* CHECK FOR LEGAL BASE B DIGIT */
if (xi < 0 || xi >= MP.b) {
mperr("*** ILLEGAL BASE B DIGIT IN CALL TO MP_MULTIPLY, POSSIBLE OVERWRITING PROBLEM ***");
- z->sign = 0;
+ mp_set_from_integer(0, z);
return;
}
@@ -1455,7 +1411,7 @@ mp_multiply(const MPNumber *x, const MPNumber *y, MPNumber *z)
int ri = r.fraction[j] + c;
if (ri < 0) {
mperr("*** INTEGER OVERFLOW IN MP_MULTIPLY, B TOO LARGE ***");
- z->sign = 0;
+ mp_set_from_integer(0, z);
return;
}
c = ri / MP.b;
@@ -1463,7 +1419,7 @@ mp_multiply(const MPNumber *x, const MPNumber *y, MPNumber *z)
}
if (c != 0) {
mperr("*** ILLEGAL BASE B DIGIT IN CALL TO MP_MULTIPLY, POSSIBLE OVERWRITING PROBLEM ***");
- z->sign = 0;
+ mp_set_from_integer(0, z);
return;
}
c = 8;
@@ -1472,7 +1428,7 @@ mp_multiply(const MPNumber *x, const MPNumber *y, MPNumber *z)
if (c != 8) {
if (xi < 0 || xi >= MP.b) {
mperr("*** ILLEGAL BASE B DIGIT IN CALL TO MP_MULTIPLY, POSSIBLE OVERWRITING PROBLEM ***");
- z->sign = 0;
+ mp_set_from_integer(0, z);
return;
}
@@ -1481,7 +1437,7 @@ mp_multiply(const MPNumber *x, const MPNumber *y, MPNumber *z)
int ri = r.fraction[j] + c;
if (ri < 0) {
mperr("*** INTEGER OVERFLOW IN MP_MULTIPLY, B TOO LARGE ***");
- z->sign = 0;
+ mp_set_from_integer(0, z);
return;
}
c = ri / MP.b;
@@ -1490,14 +1446,14 @@ mp_multiply(const MPNumber *x, const MPNumber *y, MPNumber *z)
if (c != 0) {
mperr("*** ILLEGAL BASE B DIGIT IN CALL TO MP_MULTIPLY, POSSIBLE OVERWRITING PROBLEM ***");
- z->sign = 0;
+ mp_set_from_integer(0, z);
return;
}
}
/* NORMALIZE AND ROUND RESULT */
// FIXME: Use stack variable because of mp_normalize brokeness
- for (i = 0; i < i2; i++)
+ for (i = 0; i < MP_SIZE; i++)
z->fraction[i] = r.fraction[i];
mp_normalize(z, 0);
}
@@ -1516,7 +1472,7 @@ mpmul2(const MPNumber *x, int iy, MPNumber *z, int trunc)
z->sign = x->sign;
if (z->sign == 0 || iy == 0) {
- z->sign = 0;
+ mp_set_from_integer(0, z);
return;
}
@@ -1532,7 +1488,6 @@ mpmul2(const MPNumber *x, int iy, MPNumber *z, int trunc)
z->exponent = x->exponent + 1;
}
else {
- mpchk();
mpovfl(z, "*** OVERFLOW OCCURRED IN MPMUL2 ***");
}
return;
@@ -1583,9 +1538,8 @@ mpmul2(const MPNumber *x, int iy, MPNumber *z, int trunc)
/* CHECK FOR INTEGER OVERFLOW */
if (ri < 0) {
- mpchk();
mperr("*** INTEGER OVERFLOW IN MPMUL2, B TOO LARGE ***");
- z->sign = 0;
+ mp_set_from_integer(0, z);
return;
}
@@ -1608,9 +1562,8 @@ mpmul2(const MPNumber *x, int iy, MPNumber *z, int trunc)
}
if (c < 0) {
- mpchk();
mperr("*** INTEGER OVERFLOW IN MPMUL2, B TOO LARGE ***");
- z->sign = 0;
+ mp_set_from_integer(0, z);
return;
}
@@ -1645,14 +1598,13 @@ mp_multiply_fraction(const MPNumber *x, int numerator, int denominator, MPNumber
int is, js;
if (denominator == 0) {
- mpchk();
mperr("*** ATTEMPTED DIVISION BY ZERO IN MP_MULTIPLY_FRACTION ***");
- z->sign = 0;
+ mp_set_from_integer(0, z);
return;
}
if (numerator == 0) {
- z->sign = 0;
+ mp_set_from_integer(0, z);
return;
}
@@ -1698,7 +1650,7 @@ mp_normalize(MPNumber *x, int trunc)
/* CHECK THAT SIGN = +-1 */
if (abs(x->sign) > 1) {
mperr("*** SIGN NOT 0, +1 OR -1 IN CALL TO MP_NORMALIZE, POSSIBLE OVERWRITING PROBLEM ***");
- x->sign = 0;
+ mp_set_from_integer(0, x);
return;
}
@@ -1709,7 +1661,7 @@ mp_normalize(MPNumber *x, int trunc)
/* Find how many places to shift and detect 0 fraction */
for (i = 1; i < i2 && x->fraction[i] == 0; i++);
if (i == i2) {
- x->sign = 0;
+ mp_set_from_integer(0, x);
return;
}
@@ -1777,14 +1729,10 @@ mp_normalize(MPNumber *x, int trunc)
}
/* Check for over and underflow */
- if (x->exponent > MP.m) {
+ if (x->exponent > MP.m)
mpovfl(x, "*** OVERFLOW OCCURRED IN MP_NORMALIZE ***");
- return;
- }
- if (x->exponent < -MP.m) {
+ else if (x->exponent < -MP.m)
mpunfl(x);
- return;
- }
}
@@ -1793,7 +1741,7 @@ mp_normalize(MPNumber *x, int trunc)
* (2T+6 IS ENOUGH IF N NONNEGATIVE).
*/
void
-mp_pwr_integer(const MPNumber *x, int n, MPNumber *y)
+mp_pwr_integer(const MPNumber *x, int n, MPNumber *z)
{
int n2, ns;
MPNumber t;
@@ -1801,43 +1749,41 @@ mp_pwr_integer(const MPNumber *x, int n, MPNumber *y)
n2 = n;
if (n2 < 0) {
/* N < 0 */
- mpchk();
n2 = -n2;
if (x->sign == 0) {
mperr("*** ATTEMPT TO RAISE ZERO TO NEGATIVE POWER IN CALL TO SUBROUTINE MP_PWR_INTEGER ***");
- y->sign = 0;
+ mp_set_from_integer(0, z);
return;
}
} else if (n2 == 0) {
/* N == 0, RETURN Y = 1. */
- mp_set_from_integer(1, y);
+ mp_set_from_integer(1, z);
return;
} else {
/* N > 0 */
- mpchk();
if (x->sign == 0) {
- y->sign = 0;
+ mp_set_from_integer(0, z);
return;
}
}
/* MOVE X */
- mp_set_from_mp(x, y);
+ mp_set_from_mp(x, z);
/* IF N < 0 FORM RECIPROCAL */
if (n < 0)
- mp_reciprocal(y, y);
- mp_set_from_mp(y, &t);
+ mp_reciprocal(z, z);
+ mp_set_from_mp(z, &t);
/* SET PRODUCT TERM TO ONE */
- mp_set_from_integer(1, y);
+ mp_set_from_integer(1, z);
/* MAIN LOOP, LOOK AT BITS OF N2 FROM RIGHT */
while(1) {
ns = n2;
n2 /= 2;
if (n2 << 1 != ns)
- mp_multiply(y, &t, y);
+ mp_multiply(z, &t, z);
if (n2 <= 0)
return;
@@ -1857,11 +1803,9 @@ mp_pwr(const MPNumber *x, const MPNumber *y, MPNumber *z)
{
MPNumber t;
- mpchk();
-
if (x->sign < 0) {
mperr(_("Negative X and non-integer Y not supported"));
- z->sign = 0;
+ mp_set_from_integer(0, z);
}
else if (x->sign == 0)
{
@@ -1869,7 +1813,7 @@ mp_pwr(const MPNumber *x, const MPNumber *y, MPNumber *z)
if (y->sign <= 0) {
mperr("*** X ZERO AND Y NONPOSITIVE IN CALL TO MP_PWR ***");
}
- z->sign = 0;
+ mp_set_from_integer(0, z);
}
else {
/* USUAL CASE HERE, X POSITIVE
@@ -1902,13 +1846,10 @@ mp_reciprocal(const MPNumber *x, MPNumber *z)
int ex, it0, t;
float rx;
- /* CHECK LEGALITY OF B, T, M AND MXR */
- mpchk();
-
/* MP_ADD_INTEGER REQUIRES 2T+6 WORDS. */
if (x->sign == 0) {
mperr("*** ATTEMPTED DIVISION BY ZERO IN CALL TO MP_RECIPROCAL ***");
- z->sign = 0;
+ mp_set_from_integer(0, z);
return;
}
@@ -2008,9 +1949,6 @@ mp_root(const MPNumber *x, int n, MPNumber *z)
float rx;
MPNumber t1, t2;
- /* CHECK LEGALITY OF B, T, M AND MXR */
- mpchk();
-
if (n == 1) {
mp_set_from_mp(x, z);
return;
@@ -2018,7 +1956,7 @@ mp_root(const MPNumber *x, int n, MPNumber *z)
if (n == 0) {
mperr("*** N == 0 IN CALL TO MP_ROOT ***");
- z->sign = 0;
+ mp_set_from_integer(0, z);
return;
}
@@ -2027,25 +1965,22 @@ mp_root(const MPNumber *x, int n, MPNumber *z)
/* LOSS OF ACCURACY IF NP LARGE, SO ONLY ALLOW NP <= MAX (B, 64) */
if (np > max(MP.b, 64)) {
mperr("*** ABS(N) TOO LARGE IN CALL TO MP_ROOT ***");
- z->sign = 0;
+ mp_set_from_integer(0, z);
return;
}
/* LOOK AT SIGN OF X */
if (x->sign == 0) {
/* X == 0 HERE, RETURN 0 IF N POSITIVE, ERROR IF NEGATIVE */
- z->sign = 0;
- if (n > 0)
- return;
-
- mperr("*** X == 0 AND N NEGATIVE IN CALL TO MP_ROOT ***");
- z->sign = 0;
+ mp_set_from_integer(0, z);
+ if (n <= 0)
+ mperr("*** X == 0 AND N NEGATIVE IN CALL TO MP_ROOT ***");
return;
}
if (x->sign < 0 && np % 2 == 0) {
mperr("*** X NEGATIVE AND N EVEN IN CALL TO MP_ROOT ***");
- z->sign = 0;
+ mp_set_from_integer(0, z);
return;
}
@@ -2153,14 +2088,12 @@ mp_sqrt(const MPNumber *x, MPNumber *z)
int i, i2, iy3;
MPNumber t;
- mpchk();
-
/* MP_ROOT NEEDS 4T+10 WORDS, BUT CAN OVERLAP SLIGHTLY. */
i2 = MP.t * 3 + 9;
if (x->sign < 0) {
mperr("*** X NEGATIVE IN CALL TO SUBROUTINE MP_SQRT ***");
} else if (x->sign == 0) {
- z->sign = 0;
+ mp_set_from_integer(0, z);
} else {
mp_root(x, -2, &t);
i = t.fraction[0];
[
Date Prev][
Date Next] [
Thread Prev][
Thread Next]
[
Thread Index]
[
Date Index]
[
Author Index]
