[gcalctool] Remove shared mpmul2 function, rewrite mp_xpowy_integer



commit 817af0f189a0668710f49e797f35d47ef67c497e
Author: Robert Ancell <robert ancell gmail com>
Date:   Mon Jul 13 16:12:23 2009 +1000

    Remove shared mpmul2 function, rewrite mp_xpowy_integer

 src/mp-convert.c  |   35 +++++++++++++--------------
 src/mp-internal.h |    1 -
 src/mp.c          |   68 +++++++++++++++++++---------------------------------
 3 files changed, 42 insertions(+), 62 deletions(-)
---
diff --git a/src/mp-convert.c b/src/mp-convert.c
index d09473c..3df34ce 100644
--- a/src/mp-convert.c
+++ b/src/mp-convert.c
@@ -213,30 +213,29 @@ mp_set_from_double(double dx, MPNumber *z)
  *  CHECK LEGALITY OF B, T, M AND MXR
  */
 void
-mp_set_from_integer(int ix, MPNumber *z)
+mp_set_from_integer(int x, MPNumber *z)
 {
     memset(z, 0, sizeof(MPNumber));
 
-    if (ix == 0) {
-        z->exponent = 1;
-        return;
-    }
-
-    if (ix < 0) {
-        ix = -ix;
+    if (x < 0) {
+        x = -x;
         z->sign = -1;
     }
-    else
+    else if (x > 0)
         z->sign = 1;
-
-    /* SET EXPONENT TO T */
-    z->exponent = MP_T;
-
-    /* INSERT IX */
-    z->fraction[MP_T - 1] = ix;
-
-    /* NORMALIZE BY CALLING MPMUL2 */
-    mpmul2(z, 1, z, 1);
+    else
+        z->sign = 0; /* Optimisation for indicating zero */
+
+    z->exponent = 1;
+    z->fraction[0] = x;
+    while (z->fraction[0] >= MP_BASE) {
+        int i;
+        for (i = z->exponent; i >= 0; i--)
+            z->fraction[i] = z->fraction[i-1];
+        z->fraction[0] = z->fraction[1] / MP_BASE;
+        z->fraction[1] = z->fraction[1] % MP_BASE;
+        z->exponent++;
+    }
 }
 
 /* CONVERTS THE RATIONAL NUMBER I/J TO MULTIPLE PRECISION Q. */
diff --git a/src/mp-internal.h b/src/mp-internal.h
index d47695b..d289c1a 100644
--- a/src/mp-internal.h
+++ b/src/mp-internal.h
@@ -43,7 +43,6 @@
 
 void mperr(const char *format, ...) __attribute__((format(printf, 1, 2)));
 void mpgcd(int *, int *);
-void mpmul2(const MPNumber *, int, MPNumber *, int);
 void mp_normalize(MPNumber *, int trunc);
 
 #endif /* MP_INTERNAL_H */
diff --git a/src/mp.c b/src/mp.c
index 9400f3e..69cf983 100644
--- a/src/mp.c
+++ b/src/mp.c
@@ -1232,7 +1232,7 @@ mp_multiply(const MPNumber *x, const MPNumber *y, MPNumber *z)
  *  EVEN IF SOME DIGITS ARE GREATER THAN B-1.
  *  RESULT IS ROUNDED IF TRUNC == 0, OTHERWISE TRUNCATED.
  */
-void
+static void
 mpmul2(const MPNumber *x, int iy, MPNumber *z, int trunc)
 {
     int c, i, c1, c2, j1, j2;
@@ -1380,7 +1380,7 @@ mp_multiply_fraction(const MPNumber *x, int numerator, int denominator, MPNumber
         mp_divide_integer(x, is * js, z);
     } else {
         mp_divide_integer(x, js, z);
-        mpmul2(z, is, z, 0);
+        mp_multiply_integer(z, is, z);
     }
 }
 
@@ -1512,8 +1512,8 @@ mp_pwr(const MPNumber *x, const MPNumber *y, MPNumber *z)
         return;
     }
     
-    /* 0^0 = 1 */
-    if (x->sign == 0 && y->sign == 0) {
+    /* x^0 = 1 */
+    if (y->sign == 0) {
         mp_set_from_integer(1, z);
         return;
     }
@@ -1789,6 +1789,7 @@ mp_factorial(const MPNumber *x, MPNumber *z)
     }
     if (!mp_is_natural(x)) {
         mperr("Cannot take factorial of non-natural number");
+        return;
     }
 
     /* Convert to integer - if couldn't be converted then the factorial would be too big anyway */
@@ -1840,50 +1841,31 @@ mp_xpowy(const MPNumber *x, const MPNumber *y, MPNumber *z)
 void
 mp_xpowy_integer(const MPNumber *x, int n, MPNumber *z)
 {
-    int n2, ns;
+    int i;
     MPNumber t;
-   
-    n2 = n;
-    if (n2 < 0) {
-        /* N < 0 */
-        n2 = -n2;
-        if (x->sign == 0) {
-            mperr("*** ATTEMPT TO RAISE ZERO TO NEGATIVE POWER IN CALL TO SUBROUTINE mp_xpowy_integer ***");
-            mp_set_from_integer(0, z);
-            return;
-        }
-    } else if (n2 == 0) {
-        /* N == 0, RETURN Y = 1. */
-        mp_set_from_integer(1, z);
+    
+    /* x^-n invalid */
+    if (x->sign == 0 && n < 0) {
+        mperr("*** ATTEMPT TO RAISE ZERO TO NEGATIVE POWER IN CALL TO SUBROUTINE mp_xpowy_integer ***");
+        return;
+    }
+    
+    /* 0^n = 0 */
+    if (x->sign == 0) {
+        mp_set_from_integer(0, z);
         return;
-    } else {
-        /* N > 0 */
-        if (x->sign == 0) {
-            mp_set_from_integer(0, z);
-            return;
-        }
     }
 
-    /* MOVE X */
-    mp_set_from_mp(x, z);
-
-    /* IF N < 0 FORM RECIPROCAL */
-    if (n < 0)
-        mp_reciprocal(z, z);
-    mp_set_from_mp(z, &t);
+    if (n < 0) {
+        mp_reciprocal(x, &t);
+        n = -n;
+    }
+    else
+        mp_set_from_mp(x, &t);
 
-    /* SET PRODUCT TERM TO ONE */
+    /* Multply x n times */
     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(z, &t, z);
-        if (n2 <= 0)
-            return;
-        
-        mp_multiply(&t, &t, &t);
+    for (i = 0; i < n; i++) {
+        mp_multiply(z, &t, z);
     }
 }



[Date Prev][Date Next]   [Thread Prev][Thread Next]   [Thread Index] [Date Index] [Author Index]