[gcalctool/gcalctool-newui2] ...
- From: Robert Ancell <rancell src gnome org>
- To: svn-commits-list gnome org
- Subject: [gcalctool/gcalctool-newui2] ...
- Date: Mon, 13 Jul 2009 02:38:53 +0000 (UTC)
commit 874264f5424fa2d2834e5de3aafd6aa05ceea4e6
Author: Robert Ancell <robert ancell gmail com>
Date: Sun Jul 12 13:30:28 2009 +0800
...
src/calctool.c | 1 -
src/mp-convert.c | 1 +
src/mp-equation.c | 1 +
src/mp-internal.h | 12 +--
src/mp-trigonometric.c | 17 +---
src/mp.c | 245 ++++--------------------------------------------
6 files changed, 25 insertions(+), 252 deletions(-)
---
diff --git a/src/calctool.c b/src/calctool.c
index 27a47fa..daffe2f 100644
--- a/src/calctool.c
+++ b/src/calctool.c
@@ -206,7 +206,6 @@ init_state(void)
gboolean use_default_digits = FALSE;
acc = MAX_DIGITS + 12; /* MP internal accuracy. */
- mp_init(acc);
v->error = FALSE; /* No calculator error initially. */
diff --git a/src/mp-convert.c b/src/mp-convert.c
index eccaf4e..875400a 100644
--- a/src/mp-convert.c
+++ b/src/mp-convert.c
@@ -2,6 +2,7 @@
/* $Header$
*
* Copyright (c) 1987-2008 Sun Microsystems, Inc. All Rights Reserved.
+ * Copyright (c) 2008-2009 Robert Ancell
*
* This program is free software; you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
diff --git a/src/mp-equation.c b/src/mp-equation.c
index fada9e0..fa52d22 100644
--- a/src/mp-equation.c
+++ b/src/mp-equation.c
@@ -2,6 +2,7 @@
/* $Header$
*
* Copyright (C) 2004-2008 Sami Pietila
+ * Copyright (c) 2008-2009 Robert Ancell
*
* This program is free software; you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
diff --git a/src/mp-internal.h b/src/mp-internal.h
index f6ba4ea..1de26c6 100644
--- a/src/mp-internal.h
+++ b/src/mp-internal.h
@@ -2,6 +2,7 @@
/* $Header$
*
* Copyright (c) 1987-2008 Sun Microsystems, Inc. All Rights Reserved.
+ * Copyright (c) 2008-2009 Robert Ancell
*
* This program is free software; you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
@@ -32,17 +33,8 @@
#define min(a, b) ((a) <= (b) ? (a) : (b))
#define max(a, b) ((a) >= (b) ? (a) : (b))
-/* Evil global variables that must be removed */
-struct {
- /* 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);
+//MP.t = (int) ((float) (accuracy) * log((float)10.) / log((float) MP_BASE) + (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;
diff --git a/src/mp-trigonometric.c b/src/mp-trigonometric.c
index e809c82..c46f2c4 100644
--- a/src/mp-trigonometric.c
+++ b/src/mp-trigonometric.c
@@ -2,6 +2,7 @@
/* $Header$
*
* Copyright (c) 1987-2008 Sun Microsystems, Inc. All Rights Reserved.
+ * Copyright (c) 2008-2009 Robert Ancell
*
* This program is free software; you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
@@ -645,19 +646,9 @@ mp_cosh(const MPNumber *x, MPNumber *z)
}
mp_abs(x, &t);
-
- /* IF ABS(X) TOO LARGE MP_EPOWY WILL PRINT ERROR MESSAGE
- * INCREASE M TO AVOID OVERFLOW WHEN COSH(X) REPRESENTABLE
- */
- MP.m += 2;
mp_epowy(&t, &t);
mp_reciprocal(&t, z);
mp_add(&t, z, z);
-
- /* RESTORE M. IF RESULT OVERFLOWS OR UNDERFLOWS, mp_divide_integer WILL
- * ACT ACCORDINGLY.
- */
- MP.m -= 2;
mp_divide_integer(z, 2, z);
}
@@ -693,15 +684,9 @@ mp_sinh(const MPNumber *x, MPNumber *z)
* INCREASE M TO AVOID OVERFLOW IF SINH(X) REPRESENTABLE
*/
else {
- MP.m += 2;
mp_epowy(&t2, &t2);
mp_reciprocal(&t2, z);
mp_subtract(&t2, z, z);
-
- /* RESTORE M. IF RESULT OVERFLOWS OR UNDERFLOWS, mp_divide_integer AT
- * STATEMENT 30 WILL ACT ACCORDINGLY.
- */
- MP.m -= 2;
}
/* DIVIDE BY TWO AND RESTORE SIGN */
diff --git a/src/mp.c b/src/mp.c
index 45e9a53..94e8e74 100644
--- a/src/mp.c
+++ b/src/mp.c
@@ -2,6 +2,7 @@
/* $Header$
*
* Copyright (c) 1987-2008 Sun Microsystems, Inc. All Rights Reserved.
+ * Copyright (c) 2008-2009 Robert Ancell
*
* This program is free software; you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
@@ -28,63 +29,8 @@
#include "mp-internal.h"
#include "calctool.h" // FIXME: Required for doerr() and MAXLINE
-/* 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;
-
- it = MP_BASE - 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);
-
- /* 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 *z)
-{
- /* 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.
- */
- mp_set_from_integer(0, z);
-}
+// FIXME: Re-add overflow and underflow detection
/* ROUTINE CALLED BY MP_DIVIDE AND MP_SQRT TO ENSURE THAT
@@ -122,84 +68,6 @@ mpext(int i, int j, MPNumber *x)
}
-/* SETS BASE (B) AND NUMBER OF DIGITS (T) TO GIVE THE
- * EQUIVALENT OF AT LEAST IDECPL DECIMAL DIGITS.
- * IDECPL SHOULD BE POSITIVE.
- * ITMAX2 IS THE DIMENSION OF ARRAYS USED FOR MP NUMBERS,
- * SO AN ERROR OCCURS IF THE COMPUTED T EXCEEDS ITMAX2 - 2.
- * MPSET ALSO SETS
- * MXR = MAXDR (DIMENSION OF R IN COMMON, >= T+4), AND
- * M = (W-1)/4 (MAXIMUM ALLOWABLE EXPONENT),
- * WHERE W IS THE LARGEST INTEGER OF THE FORM 2**K-1 WHICH IS
- * REPRESENTABLE IN THE MACHINE, K <= 47
- * THE COMPUTED B AND T SATISFY THE CONDITIONS
- * (T-1)*LN(B)/LN(10) >= IDECPL AND 8*B*B-1 <= W .
- * APPROXIMATELY MINIMAL T AND MAXIMAL B SATISFYING
- * THESE CONDITIONS ARE CHOSEN.
- * PARAMETERS IDECPL, ITMAX2 AND MAXDR ARE INTEGERS.
- * BEWARE - MPSET WILL CAUSE AN INTEGER OVERFLOW TO OCCUR
- * ****** IF WORDLENGTH IS LESS THAN 48 BITS.
- * IF THIS IS NOT ALLOWABLE, CHANGE THE DETERMINATION
- * OF W (DO 30 ... TO 30 W = WN) OR SET B, T, M,
- * AND MXR WITHOUT CALLING MPSET.
- * FIRST SET MXR
- */
-void
-mp_init(int accuracy)
-{
- 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
- * 47 TIMES AT MOST. IF GENUINE INTEGER WORDLENGTH
- * EXCEEDS 47 BITS THIS RESTRICTION CAN BE RELAXED.
- */
- w = 0;
- k = 0;
- for (i = 1; i <= 47; ++i) {
- int w2, wn;
-
- /* INTEGER OVERFLOW WILL EVENTUALLY OCCUR HERE
- * IF WORDLENGTH < 48 BITS
- */
- w2 = w + w;
- wn = w2 + 1;
-
- /* APPARENTLY REDUNDANT TESTS MAY BE NECESSARY ON SOME
- * MACHINES, DEPENDING ON HOW INTEGER OVERFLOW IS HANDLED
- */
- if (wn <= w || wn <= w2 || wn <= 0)
- break;
- k = i;
- w = wn;
- }
-
- /* SET MAXIMUM EXPONENT TO (W-1)/4 */
- MP.m = w / 4;
-
- /* B IS THE LARGEST POWER OF 2 SUCH THAT (8*B*B-1) <= W */
- /*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;
- while (MP.b % (10 * b) != MP.b)
- b *= 10;
- MP.b = b;*/
-}
-
-
void
mp_get_eulers(MPNumber *z)
{
@@ -519,22 +387,6 @@ mp_integer_component(const MPNumber *x, MPNumber *z)
}
}
-/* COMPARES MP NUMBER X WITH REAL NUMBER RI, RETURNING
- * +1 IF X > RI,
- * 0 IF X == RI,
- * -1 IF X < RI
- * DIMENSION OF R IN COMMON AT LEAST 2T+6
- * CHECK LEGALITY OF B, T, M AND MXR
- */
-static int
-mp_compare_mp_to_float(const MPNumber *x, float ri)
-{
- MPNumber t;
-
- /* CONVERT RI TO MULTIPLE-PRECISION AND COMPARE */
- mp_set_from_float(ri, &t);
- return mp_compare_mp_to_mp(x, &t);
-}
/* COMPARES THE MULTIPLE-PRECISION NUMBERS X AND Y,
* RETURNING +1 IF X > Y,
@@ -605,9 +457,6 @@ mp_divide(const MPNumber *x, const MPNumber *y, MPNumber *z)
return;
}
- /* INCREASE M TO AVOID OVERFLOW IN MP_RECIPROCAL */
- MP.m += 2;
-
/* FORM RECIPROCAL OF Y */
mp_reciprocal(y, &t);
@@ -621,15 +470,7 @@ mp_divide(const MPNumber *x, const MPNumber *y, MPNumber *z)
iz3 = z->fraction[0];
mpext(i, iz3, z);
- /* RESTORE M, CORRECT EXPONENT AND RETURN */
- MP.m -= 2;
z->exponent += ie;
-
- /* Check for overflow or underflow */
- if (z->exponent < -MP.m)
- mpunfl(z);
- else if (z->exponent > MP.m)
- mpovfl(z, "*** OVERFLOW OCCURRED IN MP_DIVIDE ***");
}
@@ -660,11 +501,6 @@ mp_divide_integer(const MPNumber *x, int iy, MPNumber *z)
/* CHECK FOR DIVISION BY B */
if (iy == MP_BASE) {
- if (x->exponent <= -MP.m)
- {
- mpunfl(z);
- return;
- }
mp_set_from_mp(x, z);
z->exponent--;
return;
@@ -855,7 +691,7 @@ void
mp_epowy(const MPNumber *x, MPNumber *z)
{
float r__1;
- int i, ie, ix, xs, tss;
+ int i, ix, xs, tss;
float rx, rz, rlb;
MPNumber t1, t2;
@@ -876,18 +712,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_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);
- return;
- }
-
- if (mp_compare_mp_to_float(x, (float) MP.m * rlb) > 0) {
- /* OVERFLOW HERE */
- mpovfl(z, "*** OVERFLOW IN SUBROUTINE MP_EPOWY ***");
- return;
- }
+ rlb = log((float)MP_BASE) * 1.01f;
/* NOW SAFE TO CONVERT X TO REAL */
rx = mp_cast_to_float(x);
@@ -899,9 +724,9 @@ mp_epowy(const MPNumber *x, MPNumber *z)
/* IF ABS(X) > M POSSIBLE THAT INT(X) OVERFLOWS,
* SO DIVIDE BY 32.
*/
- if (fabs(rx) > (float) MP.m) {
+ /*if (fabs(rx) > (float)MP.m) {
mp_divide_integer(&t2, 32, &t2);
- }
+ }*/
/* GET FRACTIONAL AND INTEGER PARTS OF ABS(X) */
ix = mp_cast_to_int(&t2);
@@ -944,34 +769,24 @@ mp_epowy(const MPNumber *x, MPNumber *z)
mp_multiply(z, &t2, z);
/* MUST CORRECT RESULT IF DIVIDED BY 32 ABOVE. */
- if (fabs(rx) > (float) MP.m && z->sign != 0) {
+ /*if (fabs(rx) > (float)MP.m && z->sign != 0) {
for (i = 1; i <= 5; ++i) {
- /* SAVE EXPONENT TO AVOID OVERFLOW IN MP_MULTIPLY */
+ // SAVE EXPONENT TO AVOID OVERFLOW IN MP_MULTIPLY
ie = z->exponent;
z->exponent = 0;
mp_multiply(z, z, z);
z->exponent += ie << 1;
-
- /* CHECK FOR UNDERFLOW AND OVERFLOW */
- if (z->exponent < -MP.m) {
- mpunfl(z);
- return;
- }
- if (z->exponent > MP.m) {
- mpovfl(z, "*** OVERFLOW IN SUBROUTINE MP_EPOWY ***");
- return;
- }
}
- }
+ }*/
/* CHECK THAT RELATIVE ERROR LESS THAN 0.01 UNLESS ABS(X) LARGE
* (WHEN EXP MIGHT OVERFLOW OR UNDERFLOW)
*/
- if (fabs(rx) > (float)10.0)
+ if (fabs(rx) > 10.0f)
return;
rz = mp_cast_to_float(z);
- if ((r__1 = rz - exp(rx), fabs(r__1)) < rz * (float)0.01)
+ if ((r__1 = rz - exp(rx), fabs(r__1)) < rz * 0.01f)
return;
/* THE FOLLOWING MESSAGE MAY INDICATE THAT
@@ -1013,11 +828,10 @@ mpexp1(const MPNumber *x, MPNumber *z)
}
mp_set_from_mp(x, &t1);
- rlb = log((float) MP_BASE);
+ 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 *
- (float)1.44 * rlb);
+ q = (int)(sqrt((float)MP_T * 0.48f * rlb) + (float) x->exponent * 1.44f * rlb);
/* HALVE Q TIMES */
if (q > 0) {
@@ -1268,7 +1082,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_BASE);
+ 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 */
@@ -1444,14 +1258,9 @@ mpmul2(const MPNumber *x, int iy, MPNumber *z, int trunc)
/* CHECK FOR MULTIPLICATION BY B */
if (iy == MP_BASE) {
- if (x->exponent < MP.m) {
- mp_set_from_mp(x, z);
- z->sign = z->sign;
- z->exponent = x->exponent + 1;
- }
- else {
- mpovfl(z, "*** OVERFLOW OCCURRED IN MPMUL2 ***");
- }
+ mp_set_from_mp(x, z);
+ z->sign = z->sign;
+ z->exponent = x->exponent + 1;
return;
}
}
@@ -1689,12 +1498,6 @@ mp_normalize(MPNumber *x, int trunc)
}
}
}
-
- /* Check for over and underflow */
- if (x->exponent > MP.m)
- mpovfl(x, "*** OVERFLOW OCCURRED IN MP_NORMALIZE ***");
- else if (x->exponent < -MP.m)
- mpunfl(x);
}
@@ -1817,9 +1620,6 @@ mp_reciprocal(const MPNumber *x, MPNumber *z)
ex = x->exponent;
- /* TEMPORARILY INCREASE M TO AVOID OVERFLOW */
- MP.m += 2;
-
/* SET EXPONENT TO ZERO SO RX NOT TOO LARGE OR SMALL. */
/* work-around to avoid touching X */
mp_set_from_mp(x, &tmp_x);
@@ -1827,7 +1627,7 @@ mp_reciprocal(const MPNumber *x, MPNumber *z)
rx = mp_cast_to_float(&tmp_x);
/* USE SINGLE-PRECISION RECIPROCAL AS FIRST APPROXIMATION */
- mp_set_from_float((float)1. / rx, &t1);
+ mp_set_from_float(1.0f / rx, &t1);
/* CORRECT EXPONENT OF FIRST APPROXIMATION */
t1.exponent -= ex;
@@ -1874,11 +1674,6 @@ mp_reciprocal(const MPNumber *x, MPNumber *z)
/* MOVE RESULT TO Y AND RETURN AFTER RESTORING T */
mp_set_from_mp(&t1, z);
-
- /* RESTORE M AND CHECK FOR OVERFLOW (UNDERFLOW IMPOSSIBLE) */
- MP.m -= 2;
- if (z->exponent > MP.m)
- mpovfl(z, "*** OVERFLOW OCCURRED IN MP_RECIPROCAL ***");
}
@@ -1946,8 +1741,8 @@ 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_BASE) -
- log((fabs(rx)))) / (float) np);
+ r__1 = exp(((float)(np * ex - x->exponent) * log((float)MP_BASE) -
+ log((fabs(rx)))) / (float)np);
mp_set_from_float(r__1, &t1);
/* SIGN OF APPROXIMATION SAME AS THAT OF X */
[
Date Prev][
Date Next] [
Thread Prev][
Thread Next]
[
Thread Index]
[
Date Index]
[
Author Index]
