[gcalctool/gcalctool-newui2] Stop modifying MP.t
- From: Robert Ancell <rancell src gnome org>
- To: svn-commits-list gnome org
- Subject: [gcalctool/gcalctool-newui2] Stop modifying MP.t
- Date: Sun, 12 Jul 2009 04:57:26 +0000 (UTC)
commit 2bd7628b05a60fd4e8bab22810982248d8223515
Author: Robert Ancell <robert ancell gmail com>
Date: Sat Jul 11 18:41:54 2009 +0800
Stop modifying MP.t
src/mp-trigonometric.c | 29 +++--------------------
src/mp.c | 59 +++++-------------------------------------------
src/mp.h | 14 ++++-------
3 files changed, 15 insertions(+), 87 deletions(-)
---
diff --git a/src/mp-trigonometric.c b/src/mp-trigonometric.c
index 73c2ba2..a439ba5 100644
--- a/src/mp-trigonometric.c
+++ b/src/mp-trigonometric.c
@@ -90,19 +90,13 @@ mpsin1(const MPNumber *x, MPNumber *z, int do_sin)
mp_set_from_integer(0, z);
/* POWER SERIES LOOP. REDUCE T IF POSSIBLE */
- for (; ; i+= 2) {
- int t, ts;
-
- t = MP.t + t1.exponent + 2;
- if (t <= 2)
+ for (; ; i += 2) {
+ if (MP.t + t1.exponent <= 0)
break;
- t = min(t, MP.t);
/* IF I*(I+1) IS NOT REPRESENTABLE AS AN INTEGER, THE FOLLOWING
* DIVISION BY I*(I+1) HAS TO BE SPLIT UP.
*/
- ts = MP.t;
- MP.t = t;
mp_multiply(&t2, &t1, &t1);
if (i > b2) {
mp_divide_integer(&t1, -i, &t1);
@@ -110,7 +104,6 @@ mpsin1(const MPNumber *x, MPNumber *z, int do_sin)
} else {
mp_divide_integer(&t1, -i * (i + 1), &t1);
}
- MP.t = ts;
mp_add(&t1, z, z);
if (t1.sign == 0)
@@ -154,18 +147,12 @@ mp_atan1N(int n, MPNumber *z)
/* 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)
+ if (MP.t + 2 + t2.exponent - z->exponent <= 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) {
mp_multiply_fraction(&t2, -i, i + 2, &t2);
mp_divide_integer(&t2, n, &t2);
@@ -174,7 +161,6 @@ mp_atan1N(int n, MPNumber *z)
else {
mp_multiply_fraction(&t2, -i, (i + 2)*n*n, &t2);
}
- MP.t = ts;
/* ADD TO SUM */
mp_add(&t2, z, z);
@@ -545,18 +531,11 @@ mp_atan(const MPNumber *x, MPAngleUnit unit, MPNumber *z)
/* SERIES LOOP. REDUCE T IF POSSIBLE. */
for (i = 1; ; i += 2) {
- int t, ts;
-
- t = MP.t + 2 + t2.exponent;
- if (t <= 1)
+ if (MP.t + 2 + t2.exponent <= 1)
break;
- t = min(t, MP.t);
- ts = MP.t;
- MP.t = t;
mp_multiply(&t2, &t1, &t2);
mp_multiply_fraction(&t2, -i, i + 2, &t2);
- MP.t = ts;
mp_add(z, &t2, z);
if (t2.sign == 0)
diff --git a/src/mp.c b/src/mp.c
index c153a28..f2247c2 100644
--- a/src/mp.c
+++ b/src/mp.c
@@ -880,7 +880,7 @@ void
mp_epowy(const MPNumber *x, MPNumber *z)
{
float r__1;
- int i, ie, ix, ts, xs, tss;
+ int i, ie, ix, xs, tss;
float rx, rz, rlb;
MPNumber t1, t2;
@@ -951,32 +951,19 @@ mp_epowy(const MPNumber *x, MPNumber *z)
t2.sign = 0;
for (i = 2 ; ; i++) {
- int t;
-
- t = min(tss, tss + 2 + t1.exponent);
- if (t <= 2)
+ if (min(tss, tss + 2 + t1.exponent) <= 2)
break;
- ts = MP.t;
- MP.t = t;
mp_divide_integer(&t1, i * xs, &t1);
- MP.t = ts;
-
- ts = MP.t;
- MP.t = tss;
mp_add(&t2, &t1, &t2);
- MP.t = ts;
if (t1.sign == 0)
break;
}
/* RAISE E OR 1/E TO POWER IX */
- ts = MP.t;
- MP.t = tss;
if (xs > 0)
mp_add_integer(&t2, 2, &t2);
mp_pwr_integer(&t2, ix, &t2);
- MP.t = ts;
/* MULTIPLY EXPS OF INTEGER AND FRACTIONAL PARTS */
mp_multiply(z, &t2, z);
@@ -1079,18 +1066,11 @@ mpexp1(const MPNumber *x, MPNumber *z)
/* SUM SERIES, REDUCING T WHERE POSSIBLE */
for (i = 2; ; i++) {
- int t, ts;
-
- t = MP.t + 2 + t2.exponent - z->exponent;
- if (t <= 2)
+ if (MP.t + t2.exponent - z->exponent <= 0)
break;
- t = min(t, MP.t);
- ts = MP.t;
- MP.t = t;
mp_multiply(&t1, &t2, &t2);
mp_divide_integer(&t2, i, &t2);
- MP.t = ts;
mp_add(&t2, z, z);
if (t2.sign == 0)
@@ -1230,16 +1210,13 @@ mplns(const MPNumber *x, MPNumber *z)
while(1)
{
- int ts, ts2, ts3;
+ int ts2, ts3;
- ts = MP.t;
- MP.t = t;
mpexp1(z, &t3);
mp_multiply(&t2, &t3, &t1);
mp_add(&t3, &t1, &t3);
mp_add(&t2, &t3, &t3);
mp_subtract(z, &t3, z);
- MP.t = ts;
if (t >= MP.t)
break;
@@ -1890,24 +1867,12 @@ mp_reciprocal(const MPNumber *x, MPNumber *z)
/* MAIN ITERATION LOOP */
if (t <= MP.t) {
while(1) {
- int ts, ts2, ts3;
+ int ts2, ts3;
- ts = MP.t;
- MP.t = t;
mp_multiply(x, &t1, &t2);
mp_add_integer(&t2, -1, &t2);
- MP.t = ts;
-
- /* TEMPORARILY REDUCE T */
- ts = MP.t;
- MP.t = (t + it0) / 2;
mp_multiply(&t1, &t2, &t2);
- MP.t = ts;
-
- ts = MP.t;
- MP.t = t;
mp_subtract(&t1, &t2, &t1);
- MP.t = ts;
if (t >= MP.t)
break;
@@ -2029,26 +1994,14 @@ mp_root(const MPNumber *x, int n, MPNumber *z)
/* MAIN ITERATION LOOP */
while(1) {
- int ts, ts2, ts3;
+ int ts2, ts3;
- ts = MP.t;
- MP.t = t;
mp_pwr_integer(&t1, np, &t2);
mp_multiply(x, &t2, &t2);
mp_add_integer(&t2, -1, &t2);
- MP.t = ts;
-
- /* TEMPORARILY REDUCE T */
- ts = MP.t;
- MP.t = (t + it0) / 2;
mp_multiply(&t1, &t2, &t2);
mp_divide_integer(&t2, np, &t2);
- MP.t = ts;
-
- ts = MP.t;
- MP.t = t;
mp_subtract(&t1, &t2, &t1);
- MP.t = ts;
/* FOLLOWING LOOP ALMOST DOUBLES T (POSSIBLE BECAUSE
* NEWTONS METHOD HAS 2ND ORDER CONVERGENCE).
diff --git a/src/mp.h b/src/mp.h
index 2e8433f..a58a432 100644
--- a/src/mp.h
+++ b/src/mp.h
@@ -43,7 +43,10 @@
/* Size of the multiple precision values */
#define MP_SIZE 1000
-/* Object for a high precision floating point number representation */
+/* Object for a high precision floating point number representation
+ *
+ * x = sign * (MP.b^(exponent-1) + MP.b^(exponent-2) + ...)
+ */
typedef struct
{
/* Sign (+1, -1) or 0 for the value zero */
@@ -52,14 +55,7 @@ typedef struct
/* Exponent (to base MP.b) */
int exponent;
- /* Normalized fraction
- *
- * If exponent > 0, contains 'exponent' elements
- * x = sign * (fraction[0]*MP.b^(exponent-1) + fraction[1]*MP.b^(exponent-2) + ...)
- *
- * If exponent <= 0, contains 1-'exponent' elements
- * x = sign * (... + fraction[0]*MP.b^(exponent-2) + fraction[1]*MP.b(^exponent-1))
- */
+ /* Normalized fraction */
int fraction[MP_SIZE];
} MPNumber;
[
Date Prev][
Date Next] [
Thread Prev][
Thread Next]
[
Thread Index]
[
Date Index]
[
Author Index]