[gcalctool] Tidy up MP number indexing

[Robert Ancell <rancell src gnome org>]

commit d3066c14da8f91aca87b74539898a3a518a32f63
Author: Robert Ancell <robert ancell gmail com>
Date:   Wed May 6 11:40:12 2009 +1000

    Tidy up MP number indexing
---
 gcalctool/mp-convert.c |   87 +++++++++++++++++++++++++----------------------
 gcalctool/mp.c         |   57 +++++++++++++------------------
 2 files changed, 70 insertions(+), 74 deletions(-)

diff --git a/gcalctool/mp-convert.c b/gcalctool/mp-convert.c
index dcf1c18..357c58f 100644
--- a/gcalctool/mp-convert.c
+++ b/gcalctool/mp-convert.c
@@ -286,19 +286,20 @@ mp_set_from_fraction(int i, int j, MPNumber *q)
 int
 mp_cast_to_int(const MPNumber *x)
 {
-    int i, j, k, j1, x2, kx, xs, izs, ret_val = 0;
+    int i, j, x2, xs, ret_val = 0;
 
+    /* RETURN 0 IF X = 0 OR IF NUMBER FRACTION */
     xs = x->sign;
-    /* RETURN 0 IF X = 0 OR IF NUMBER FRACTION */    
-    if (xs == 0  ||  x->exponent <= 0)
+    if (xs == 0 || x->exponent <= 0)
         return 0;
 
     x2 = x->exponent;
-    for (i = 1; i <= x2; ++i) {
+    for (i = 0; i < x2; i++) {
+        int izs;
         izs = ret_val;
         ret_val = MP.b * ret_val;
-        if (i <= MP.t)
-            ret_val += x->fraction[i - 1];
+        if (i < MP.t)
+            ret_val += x->fraction[i];
 
         /* CHECK FOR SIGNS OF INTEGER OVERFLOW */
         if (ret_val <= 0 || ret_val <= izs)
@@ -309,12 +310,13 @@ mp_cast_to_int(const MPNumber *x)
      *  HAVE OCCURRED).
      */
     j = ret_val;
-    for (i = 1; i <= x2; ++i) {
+    for (i = x2 - 1; i >= 0; i--) {
+        int j1, kx;
+        
         j1 = j / MP.b;
-        k = x2 + 1 - i;
         kx = 0;
-        if (k <= MP.t)
-            kx = x->fraction[k - 1];
+        if (i < MP.t)
+            kx = x->fraction[i];
         if (kx != j - MP.b * j1)
             return 0;
         j = j1;
@@ -323,8 +325,7 @@ mp_cast_to_int(const MPNumber *x)
         return 0;
 
     /* RESULT CORRECT SO RESTORE SIGN AND RETURN */
-    ret_val = xs * ret_val;
-    return ret_val;
+    return xs * ret_val;
 
     /* Old comment about returning zero: */
     /*  HERE OVERFLOW OCCURRED (OR X WAS UNNORMALIZED), SO
@@ -335,22 +336,29 @@ mp_cast_to_int(const MPNumber *x)
 static double
 mppow_ri(float ap, int bp)
 {
-    double pow = 1.0;
+    double pow;
+    
+    if (bp == 0)
+        return 1.0;
 
-    if (bp != 0) { 
-        if (bp < 0) {
-            if (ap == 0) return(pow);
-            bp = -bp;
-            ap = 1/ap;
-        }
-        for (;;) { 
-            if (bp & 01)  pow *= ap;
-            if (bp >>= 1) ap *= ap;
-            else break;
-        }
+    if (bp < 0) {
+        if (ap == 0)
+            return 1.0;
+        bp = -bp;
+        ap = 1 / ap;
     }
-
-    return(pow);
+    
+    pow = 1.0;
+    for (;;) { 
+        if (bp & 01)
+            pow *= ap;
+        if (bp >>= 1)
+            ap *= ap;
+        else
+            break;
+    }
+    
+    return pow;
 }
 
 /*  CONVERTS MULTIPLE-PRECISION X TO SINGLE-PRECISION.
@@ -361,28 +369,24 @@ mppow_ri(float ap, int bp)
 float
 mp_cast_to_float(const MPNumber *x)
 {
-    float rz = 0.0;
-
-    int i, tm = 0;
-    float rb, rz2;
+    int i;
+    float rb, rz = 0.0;
     
     mpchk(1, 4);
     if (x->sign == 0)
         return 0.0;
 
     rb = (float) MP.b;
-    for (i = 1; i <= MP.t; ++i) {
-        rz = rb * rz + (float) x->fraction[i - 1];
-        tm = i;
+    for (i = 0; i < MP.t; i++) {
+        rz = rb * rz + (float)x->fraction[i];
 
         /* CHECK IF FULL SINGLE-PRECISION ACCURACY ATTAINED */
-        rz2 = rz + (float) 1.0;
-        if (rz2 <= rz)
+        if (rz + 1.0f <= rz)
             break;
     }
 
     /* NOW ALLOW FOR EXPONENT */
-    rz *= mppow_ri(rb, x->exponent - tm);
+    rz *= mppow_ri(rb, x->exponent - i - 1);
 
     /* CHECK REASONABLENESS OF RESULT */
     /* LHS SHOULD BE <= 0.5, BUT ALLOW FOR SOME ERROR IN ALOG */
@@ -397,6 +401,7 @@ mp_cast_to_float(const MPNumber *x)
 
     if (x->sign < 0)
         rz = -(double)(rz);
+
     return rz;
 }
 
@@ -439,17 +444,17 @@ mp_cast_to_double(const MPNumber *x)
         return 0.0;
 
     db = (double) MP.b;
-    for (i = 1; i <= MP.t; ++i) {
-        ret_val = db * ret_val + (double) x->fraction[i - 1];
+    for (i = 0; i < MP.t; i++) {
+        ret_val = db * ret_val + (double) x->fraction[i];
         tm = i;
 
         /* CHECK IF FULL DOUBLE-PRECISION ACCURACY ATTAINED */
-        dz2 = ret_val + 1.;
+        dz2 = ret_val + 1.0;
 
         /*  TEST BELOW NOT ALWAYS EQUIVALENT TO - IF (DZ2.LE.DZ) GO TO 20,
          *  FOR EXAMPLE ON CYBER 76.
          */
-        if (dz2 - ret_val <= 0.)
+        if (dz2 - ret_val <= 0.0)
             break;
     }
 
diff --git a/gcalctool/mp.c b/gcalctool/mp.c
index cf8d6d7..f887457 100644
--- a/gcalctool/mp.c
+++ b/gcalctool/mp.c
@@ -649,10 +649,10 @@ mpdivi(const MPNumber *x, int iy, MPNumber *z)
     if (iy < b2) {
         /* LOOK FOR FIRST NONZERO DIGIT IN QUOTIENT */
         do {
-            ++i;
             c = MP.b * c;
-            if (i <= MP.t)
-                c += x->fraction[i - 1];
+            if (i < MP.t)
+                c += x->fraction[i];
+            i++;
             r1 = c / iy;
             if (r1 < 0)
                 goto L210;
@@ -665,11 +665,11 @@ mpdivi(const MPNumber *x, int iy, MPNumber *z)
         kh = 2;
         if (i < MP.t) {
             kh = MP.t + 1 - i;
-            for (k = 2; k <= kh; ++k) {
-                ++i;
-                c += x->fraction[i - 1];
-                MP.r[k - 1] = c / iy;
-                c = MP.b * (c - iy * MP.r[k - 1]);
+            for (k = 1; k < kh; k++) {
+                c += x->fraction[i];
+                MP.r[k] = c / iy;
+                c = MP.b * (c - iy * MP.r[k]);
+                i++;
             }
             if (c < 0)
                 goto L210;
@@ -690,33 +690,25 @@ mpdivi(const MPNumber *x, int iy, MPNumber *z)
     }
     
     /* HERE NEED SIMULATED DOUBLE-PRECISION DIVISION */
-    c2 = 0;
     j1 = iy / MP.b;
-    j2 = iy - j1 * MP.b;
-    j11 = j1 + 1;
 
     /* LOOK FOR FIRST NONZERO DIGIT */
-    while(1) {
-        i++;
+    c2 = 0;
+    j2 = iy - j1 * MP.b;
+    do {
         c = MP.b * c + c2;
-        c2 = 0;
-        if (i <= MP.t)
-            c2 = x->fraction[i - 1];
         i__1 = c - j1;
-        if (i__1 < 0)
-            continue;
-        else if (i__1 == 0) {
-            if (c2 < j2)
-                continue;
-        }
-        break;
-    }
+        c2 = i < MP.t ? x->fraction[i] : 0;
+        i++;
+    } while (i__1 < 0 || (i__1 == 0 && c2 < j2));
 
     /* COMPUTE T+4 QUOTIENT DIGITS */
     re = re + 1 - i;
+    i--;
     k = 1;
 
     /* MAIN LOOP FOR LARGE ABS(IY) CASE */
+    j11 = j1 + 1;
     while(1) {
         /* GET APPROXIMATE QUOTIENT FIRST */
         ir = c / j11;
@@ -736,8 +728,9 @@ mpdivi(const MPNumber *x, int iy, MPNumber *z)
             iq += iy;
         }
 
-        if (i <= MP.t)
-            iq += x->fraction[i - 1];
+        if (i < MP.t)
+            iq += x->fraction[i];
+        i++;
         iqj = iq / iy;
 
         /* R(K) = QUOTIENT, C = REMAINDER */
@@ -752,7 +745,6 @@ mpdivi(const MPNumber *x, int iy, MPNumber *z)
             mp_get_normalized_register(rs, &re, z, 0);
             return;
         }
-        ++i;
     }
 
 L210:
@@ -1356,8 +1348,8 @@ mpmaxr(MPNumber *x)
     it = MP.b - 1;
 
     /* SET FRACTION DIGITS TO B-1 */
-    for (i = 1; i <= MP.t; i++)
-        x->fraction[i - 1] = it;
+    for (i = 0; i < MP.t; i++)
+        x->fraction[i] = it;
 
     /* SET SIGN AND EXPONENT */
     x->sign = 1;
@@ -1395,7 +1387,6 @@ void
 mpmul(const MPNumber *x, const MPNumber *y, MPNumber *z)
 {
     int i__1;
-    
     int c, i, j, i2, j1, re, ri, xi, rs, i2p;
 
     mpchk(1, 4);
@@ -1420,15 +1411,15 @@ mpmul(const MPNumber *x, const MPNumber *y, MPNumber *z)
     /* PERFORM MULTIPLICATION */
     c = 8;
     i__1 = MP.t;
-    for (i = 1; i <= i__1; ++i) {
-        xi = x->fraction[i - 1];
+    for (i = 0; i < i__1; i++) {
+        xi = x->fraction[i];
 
         /* FOR SPEED, PUT THE NUMBER WITH MANY ZEROS FIRST */
         if (xi == 0)
             continue;
 
         /* Computing MIN */
-        mpmlp(&MP.r[i], y->fraction, xi, min(MP.t, i2 - i));
+        mpmlp(&MP.r[i+1], y->fraction, xi, min(MP.t, i2 - i - 1));
         --c;
         if (c > 0)
             continue;