Port to GNU MPFR/MPC Library

For further information please visit:
https://www.mpfr.org/
http://www.multiprecision.org/mpc

1 files changed, 28 insertions, 408 deletions

diff --git a/src/mp-convert.c b/src/mp-convert.c
index f41c874..1770853 100644
--- a/src/mp-convert.c
+++ b/src/mp-convert.c
@@ -17,226 +17,38 @@
 #include <langinfo.h>
 
 #include "mp.h"
-#include "mp-private.h"
 
 void
 mp_set_from_mp(const MPNumber *x, MPNumber *z)
 {
-    if (z != x)
-        memcpy(z, x, sizeof(MPNumber));
-}
-
-
-void
-mp_set_from_float(float rx, MPNumber *z)
-{
-    int i, k, ib, ie, tp;
-    float rj;
-
-    mp_set_from_integer(0, z);
-
-    /* CHECK SIGN */
-    if (rx < 0.0f) {
-        z->sign = -1;
-        rj = -(double)(rx);
-    } else if (rx > 0.0f) {
-        z->sign = 1;
-        rj = rx;
-    } else {
-        /* IF RX = 0E0 RETURN 0 */
-        mp_set_from_integer(0, z);
-        return;
-    }
-
-    /* INCREASE IE AND DIVIDE RJ BY 16. */
-    ie = 0;
-    while (rj >= 1.0f) {
-        ++ie;
-        rj *= 0.0625f;
-    }
-    while (rj < 0.0625f) {
-        --ie;
-        rj *= 16.0f;
-    }
-
-    /*  NOW RJ IS DY DIVIDED BY SUITABLE POWER OF 16.
-     *  SET EXPONENT TO 0
-     */
-    z->exponent = 0;
-
-    /* CONVERSION LOOP (ASSUME SINGLE-PRECISION OPS. EXACT) */
-    for (i = 0; i < MP_T + 4; i++) {
-        rj *= (float) MP_BASE;
-        z->fraction[i] = (int) rj;
-        rj -= (float) z->fraction[i];
-    }
-
-    /* NORMALIZE RESULT */
-    mp_normalize(z);
-
-    /* Computing MAX */
-    ib = max(MP_BASE * 7 * MP_BASE, 32767) / 16;
-    tp = 1;
-
-    /* NOW MULTIPLY BY 16**IE */
-    if (ie < 0)  {
-        k = -ie;
-        for (i = 1; i <= k; ++i) {
-            tp <<= 4;
-            if (tp <= ib && tp != MP_BASE && i < k)
-                continue;
-            mp_divide_integer(z, tp, z);
-            tp = 1;
-        }
-    } else if (ie > 0)  {
-        for (i = 1; i <= ie; ++i) {
-            tp <<= 4;
-            if (tp <= ib && tp != MP_BASE && i < ie)
-                continue;
-            mp_multiply_integer(z, tp, z);
-            tp = 1;
-        }
-    }
+    mpc_set(z->num, x->num, MPC_RNDNN);
 }
 
 
 void
 mp_set_from_double(double dx, MPNumber *z)
 {
-    int i, k, ib, ie, tp;
-    double dj;
-
-    mp_set_from_integer(0, z);
-
-    /* CHECK SIGN */
-    if (dx < 0.0)  {
-        z->sign = -1;
-        dj = -dx;
-    } else if (dx > 0.0)  {
-        z->sign = 1;
-        dj = dx;
-    } else {
-        mp_set_from_integer(0, z);
-        return;
-    }
-
-    /* INCREASE IE AND DIVIDE DJ BY 16. */
-    for (ie = 0; dj >= 1.0; ie++)
-      dj *= 1.0/16.0;
-
-    for ( ; dj < 1.0/16.0; ie--)
-      dj *= 16.0;
-
-    /*  NOW DJ IS DY DIVIDED BY SUITABLE POWER OF 16
-     *  SET EXPONENT TO 0
-     */
-    z->exponent = 0;
-
-    /* CONVERSION LOOP (ASSUME DOUBLE-PRECISION OPS. EXACT) */
-    for (i = 0; i < MP_T + 4; i++) {
-        dj *= (double) MP_BASE;
-        z->fraction[i] = (int) dj;
-        dj -= (double) z->fraction[i];
-    }
-
-    /* NORMALIZE RESULT */
-    mp_normalize(z);
-
-    /* Computing MAX */
-    ib = max(MP_BASE * 7 * MP_BASE, 32767) / 16;
-    tp = 1;
-
-    /* NOW MULTIPLY BY 16**IE */
-    if (ie < 0) {
-        k = -ie;
-        for (i = 1; i <= k; ++i) {
-            tp <<= 4;
-            if (tp <= ib && tp != MP_BASE && i < k)
-                continue;
-            mp_divide_integer(z, tp, z);
-            tp = 1;
-        }
-    } else if (ie > 0) {
-        for (i = 1; i <= ie; ++i) {
-            tp <<= 4;
-            if (tp <= ib && tp != MP_BASE && i < ie)
-                continue;
-            mp_multiply_integer(z, tp, z);
-            tp = 1;
-        }
-    }
+    mpc_set_d(z->num, dx, MPC_RNDNN);
 }
 
 
 void
 mp_set_from_integer(int64_t x, MPNumber *z)
 {
-    int i;
-
-    memset(z, 0, sizeof(MPNumber));
-
-    if (x == 0) {
-        z->sign = 0;
-        return;
-    }
-
-    if (x < 0) {
-        x = -x;
-        z->sign = -1;
-    }
-    else if (x > 0)
-        z->sign = 1;
-
-    while (x != 0) {
-        z->fraction[z->exponent] = x % MP_BASE;
-        z->exponent++;
-        x /= MP_BASE;
-    }
-    for (i = 0; i < z->exponent / 2; i++) {
-        int t = z->fraction[i];
-        z->fraction[i] = z->fraction[z->exponent - i - 1];
-        z->fraction[z->exponent - i - 1] = t;
-    }
+    mpc_set_si(z->num, x, MPC_RNDNN);
 }
 
 
 void
 mp_set_from_unsigned_integer(uint64_t x, MPNumber *z)
 {
-    int i;
-
-    mp_set_from_integer(0, z);
-
-    if (x == 0) {
-        z->sign = 0;
-        return;
-    }
-    z->sign = 1;
-
-    while (x != 0) {
-        z->fraction[z->exponent] = x % MP_BASE;
-        x = x / MP_BASE;
-        z->exponent++;
-    }
-    for (i = 0; i < z->exponent / 2; i++) {
-        int t = z->fraction[i];
-        z->fraction[i] = z->fraction[z->exponent - i - 1];
-        z->fraction[z->exponent - i - 1] = t;
-    }
+    mpc_set_ui(z->num, x, MPC_RNDNN);
 }
 
 
 void
 mp_set_from_fraction(int64_t numerator, int64_t denominator, MPNumber *z)
 {
-    mp_gcd(&numerator, &denominator);
-
-    if (denominator == 0) {
-        mperr("*** J == 0 IN CALL TO MP_SET_FROM_FRACTION ***\n");
-        mp_set_from_integer(0, z);
-        return;
-    }
-
     if (denominator < 0) {
         numerator = -numerator;
         denominator = -denominator;
@@ -251,252 +63,54 @@ mp_set_from_fraction(int64_t numerator, int64_t denominator, MPNumber *z)
 void
 mp_set_from_polar(const MPNumber *r, MPAngleUnit unit, const MPNumber *theta, MPNumber *z)
 {
-    MPNumber x, y;
+    MPNumber x = mp_new();
+    MPNumber y = mp_new();
 
     mp_cos(theta, unit, &x);
     mp_multiply(&x, r, &x);
     mp_sin(theta, unit, &y);
     mp_multiply(&y, r, &y);
     mp_set_from_complex(&x, &y, z);
-}
 
+    mp_clear(&x);
+    mp_clear(&y);
+}
 
 void
 mp_set_from_complex(const MPNumber *x, const MPNumber *y, MPNumber *z)
 {
-    /* NOTE: Do imaginary component first as z may be x or y */
-    z->im_sign = y->sign;
-    z->im_exponent = y->exponent;
-    memcpy(z->im_fraction, y->fraction, sizeof(int) * MP_SIZE);
-
-    z->sign = x->sign;
-    z->exponent = x->exponent;
-    if (z != x)
-        memcpy(z->fraction, x->fraction, sizeof(int) * MP_SIZE);
+    mpc_set_fr_fr(z->num, mpc_realref(x->num), mpc_realref(y->num), MPC_RNDNN);
 }
 
-
 void
 mp_set_from_random(MPNumber *z)
 {
     mp_set_from_double(drand48(), z);
 }
 
-
 int64_t
-mp_cast_to_int(const MPNumber *x)
+mp_to_integer(const MPNumber *x)
 {
-    int i;
-    int64_t z = 0, v;
-
-    /* |x| <= 1 */
-    if (x->sign == 0 || x->exponent <= 0)
-        return 0;
-
-    /* Multiply digits together */
-    for (i = 0; i < x->exponent; i++) {
-        int64_t t;
-
-        t = z;
-        z = z * MP_BASE + x->fraction[i];
-
-        /* Check for overflow */
-        if (z <= t)
-            return 0;
-    }
-
-    /* Validate result */
-    v = z;
-    for (i = x->exponent - 1; i >= 0; i--) {
-        int64_t digit;
-
-        /* Get last digit */
-        digit = v - (v / MP_BASE) * MP_BASE;
-        if (x->fraction[i] != digit)
-            return 0;
-
-        v /= MP_BASE;
-    }
-    if (v != 0)
-        return 0;
-
-    return x->sign * z;
+    return mpfr_get_si(mpc_realref(x->num), MPFR_RNDN);
 }
 
 
 uint64_t
-mp_cast_to_unsigned_int(const MPNumber *x)
-{
-    int i;
-    uint64_t z = 0, v;
-
-    /* x <= 1 */
-    if (x->sign <= 0 || x->exponent <= 0)
-        return 0;
-
-    /* Multiply digits together */
-    for (i = 0; i < x->exponent; i++) {
-        uint64_t t;
-
-        t = z;
-        z = z * MP_BASE + x->fraction[i];
-
-        /* Check for overflow */
-        if (z <= t)
-            return 0;
-    }
-
-    /* Validate result */
-    v = z;
-    for (i = x->exponent - 1; i >= 0; i--) {
-        uint64_t digit;
-
-        /* Get last digit */
-        digit = v - (v / MP_BASE) * MP_BASE;
-        if (x->fraction[i] != digit)
-            return 0;
-
-        v /= MP_BASE;
-    }
-    if (v != 0)
-        return 0;
-
-    return z;
-}
-
-
-static double
-mppow_ri(float ap, int bp)
+mp_to_unsigned_integer(const MPNumber *x)
 {
-    double pow;
-
-    if (bp == 0)
-        return 1.0;
-
-    if (bp < 0) {
-        if (ap == 0)
-            return 1.0;
-        bp = -bp;
-        ap = 1 / ap;
-    }
-
-    pow = 1.0;
-    for (;;) {
-        if (bp & 01)
-            pow *= ap;
-        if (bp >>= 1)
-            ap *= ap;
-        else
-            break;
-    }
-
-    return pow;
+    return mpfr_get_ui(mpc_realref(x->num), MPFR_RNDN);
 }
 
-
 float
-mp_cast_to_float(const MPNumber *x)
-{
-    int i;
-    float rz = 0.0;
-
-    if (mp_is_zero(x))
-        return 0.0;
-
-    for (i = 0; i < MP_T; i++) {
-        rz = (float) MP_BASE * rz + (float)x->fraction[i];
-
-        /* CHECK IF FULL SINGLE-PRECISION ACCURACY ATTAINED */
-        if (rz + 1.0f <= rz)
-            break;
-    }
-
-    /* NOW ALLOW FOR EXPONENT */
-    rz *= mppow_ri((float) MP_BASE, x->exponent - i - 1);
-
-    /* CHECK REASONABLENESS OF RESULT */
-    /* LHS SHOULD BE <= 0.5, BUT ALLOW FOR SOME ERROR IN ALOG */
-    if (rz <= (float)0. ||
-        fabs((float) x->exponent - (log(rz) / log((float) MP_BASE) + (float).5)) > (float).6) {
-        /*  FOLLOWING MESSAGE INDICATES THAT X IS TOO LARGE OR SMALL -
-         *  TRY USING MPCMRE INSTEAD.
-         */
-        mperr("*** FLOATING-POINT OVER/UNDER-FLOW IN MP_CAST_TO_FLOAT ***\n");
-        return 0.0;
-    }
-
-    if (x->sign < 0)
-        rz = -(double)(rz);
-
-    return rz;
-}
-
-
-static double
-mppow_di(double ap, int bp)
+mp_to_float(const MPNumber *x)
 {
-    double pow = 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;
-        }
-    }
-
-    return(pow);
+    return mpfr_get_flt(mpc_realref(x->num), MPFR_RNDN);  
 }
 
-
 double
-mp_cast_to_double(const MPNumber *x)
+mp_to_double(const MPNumber *x)
 {
-    int i, tm = 0;
-    double d__1, dz2, ret_val = 0.0;
-
-    if (mp_is_zero(x))
-        return 0.0;
-
-    for (i = 0; i < MP_T; i++) {
-        ret_val = (double) MP_BASE * ret_val + (double) x->fraction[i];
-        tm = i;
-
-        /* CHECK IF FULL DOUBLE-PRECISION ACCURACY ATTAINED */
-        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.0)
-            break;
-    }
-
-    /* NOW ALLOW FOR EXPONENT */
-    ret_val *= mppow_di((double) MP_BASE, x->exponent - tm - 1);
-
-    /* CHECK REASONABLENESS OF RESULT. */
-    /* LHS SHOULD BE .LE. 0.5 BUT ALLOW FOR SOME ERROR IN DLOG */
-    if (ret_val <= 0. ||
-        ((d__1 = (double) ((float) x->exponent) - (log(ret_val) / log((double)
-                ((float) MP_BASE)) + .5), abs(d__1)) > .6)) {
-        /*  FOLLOWING MESSAGE INDICATES THAT X IS TOO LARGE OR SMALL -
-         *  TRY USING MPCMDE INSTEAD.
-         */
-        mperr("*** FLOATING-POINT OVER/UNDER-FLOW IN MP_CAST_TO_DOUBLE ***\n");
-        return 0.0;
-    }
-    else
-    {
-        if (x->sign < 0)
-            ret_val = -ret_val;
-        return ret_val;
-    }
+    return mpfr_get_d(mpc_realref(x->num), MPFR_RNDN);  
 }
 
 static int
@@ -586,7 +200,6 @@ set_from_sexagesimal(const char *str, int length, MPNumber *z)
 {
     int degrees = 0, minutes = 0;
     char seconds[length+1];
-    MPNumber t;
     int n_matched;
 
     seconds[0] = '\0';
@@ -594,6 +207,7 @@ set_from_sexagesimal(const char *str, int length, MPNumber *z)
 
     if (n_matched < 1)
         return true;
+    MPNumber t = mp_new();
     mp_set_from_integer(degrees, z);
     if (n_matched > 1) {
         mp_set_from_integer(minutes, &t);
@@ -605,6 +219,7 @@ set_from_sexagesimal(const char *str, int length, MPNumber *z)
         mp_divide_integer(&t, 3600, &t);
         mp_add(z, &t, z);
     }
+    mp_clear(&t);
 
     return false;
 }
@@ -674,9 +289,10 @@ mp_set_from_string(const char *str, int default_base, MPNumber *z)
         }
     }
     if (fractions[i] != NULL) {
-        MPNumber fraction;
+        MPNumber fraction = mp_new();
         mp_set_from_fraction(numerators[i], denominators[i], &fraction);
         mp_add(z, &fraction, z);
+        mp_clear(&fraction);
     }
 
     if (*c == '.') {
@@ -686,7 +302,8 @@ mp_set_from_string(const char *str, int default_base, MPNumber *z)
 
     /* Convert fractional part */
     if (has_fraction) {
-        MPNumber numerator, denominator;
+        MPNumber numerator = mp_new();
+        MPNumber denominator = mp_new();
 
         mp_set_from_integer(0, &numerator);
         mp_set_from_integer(1, &denominator);
@@ -697,6 +314,8 @@ mp_set_from_string(const char *str, int default_base, MPNumber *z)
         }
         mp_divide(&numerator, &denominator, &numerator);
         mp_add(z, &numerator, z);
+        mp_clear(&numerator);
+        mp_clear(&denominator);
     }
 
     if (c != end) {
@@ -704,10 +323,11 @@ mp_set_from_string(const char *str, int default_base, MPNumber *z)
     }
 
     if (multiplier != 0) {
-        MPNumber t;
+        MPNumber t = mp_new();
         mp_set_from_integer(10, &t);
         mp_xpowy_integer(&t, multiplier, &t);
         mp_multiply(z, &t, z);
+        mp_clear(&t);
     }
 
     if (negate == 1)