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)