[gcalctool] Tidy up mp.c

[Robert Ancell <rancell src gnome org>]

commit 70c27cb140830daf3d65c3a2b2656dfaf902c4c2
Author: Robert Ancell <robert ancell gmail com>
Date:   Mon Nov 16 16:40:36 2009 -0600

    Tidy up mp.c

 src/mp-trigonometric.c |   39 +++++++++++++++++++++------------------
 src/mp.c               |   18 ++++++++++--------
 src/mp.h               |    6 +++---
 3 files changed, 34 insertions(+), 29 deletions(-)
---
diff --git a/src/mp-trigonometric.c b/src/mp-trigonometric.c
index 5e278d9..b1e6f23 100644
--- a/src/mp-trigonometric.c
+++ b/src/mp-trigonometric.c
@@ -129,7 +129,7 @@ mpsin1(const MPNumber *x, MPNumber *z, int do_sin)
 
     /* Taylor series */
     /* POWER SERIES LOOP.  REDUCE T IF POSSIBLE */
-    b2 = max(MP_BASE, 64) << 1;
+    b2 = 2 * max(MP_BASE, 64);
     do {
         if (MP_T + t1.exponent <= 0)
             break;
@@ -250,7 +250,7 @@ mp_cos(const MPNumber *x, MPAngleUnit unit, MPNumber *z)
 
     convert_to_radians(x, unit, z);
 
-    /* Use power series if -1 >= x >= 1 */
+    /* Use power series if |x| <= 1 */
     mp_abs(z, z);
     if (mp_compare_mp_to_int(z, 1) <= 0) {
         mpsin1(z, z, 0);
@@ -291,14 +291,14 @@ mp_asin(const MPNumber *x, MPAngleUnit unit, MPNumber *z)
 {
     MPNumber t1, t2;
 
-    /* asin(0) = 0 */
+    /* asin�¹(0) = 0 */
     if (mp_is_zero(x)) {
         mp_set_from_integer(0, z);
         return;
     }
 
+    /* sinâ?»Â¹(x) = tanâ?»Â¹(x / â??(1 - x²)), |x| < 1 */
     if (x->exponent <= 0) {
-        /* HERE ABS(X) < 1,  SO USE ARCTAN(X/SQRT(1 - X^2)) */
         mp_set_from_integer(1, &t1);
         mp_set_from_mp(&t1, &t2);
         mp_subtract(&t1, x, &t1);
@@ -310,18 +310,18 @@ mp_asin(const MPNumber *x, MPAngleUnit unit, MPNumber *z)
         return;
     }
 
-    /* HERE ABS(X) >= 1.  SEE IF X == +-1 */
+    /* sinâ?»Â¹(1) = Ï?/2, sinâ?»Â¹(-1) = -Ï?/2 */
     mp_set_from_integer(x->sign, &t2);
-    if (!mp_is_equal(x, &t2)) {
-        /* Translators: Error displayed when inverse sine value is undefined */
-        mperr(_("Inverse sine not defined for values outside [-1, 1]"));
+    if (mp_is_equal(x, &t2)) {
+        mp_get_pi(z);
+        mp_divide_integer(z, 2 * t2.sign, z);
+        convert_from_radians(z, unit, z);
+        return;
     }
 
-    /* X == +-1 SO RETURN +-PI/2 */
-    mp_get_pi(z);
-    mp_divide_integer(z, t2.sign << 1, z);
-
-    convert_from_radians(z, unit, z);
+    /* Translators: Error displayed when inverse sine value is undefined */
+    mperr(_("Inverse sine not defined for values outside [-1, 1]"));
+    mp_set_from_integer(0, z);
 }
 
 
@@ -346,6 +346,7 @@ mp_acos(const MPNumber *x, MPAngleUnit unit, MPNumber *z)
     } else if (mp_is_equal(x, &MPn1)) {
         mp_set_from_mp(&pi, z);
     } else {
+        /* cosâ?»Â¹(x) = tanâ?»Â¹(â??(1 - x²) / x) */
         mp_multiply(x, x, &t2);
         mp_subtract(&t1, &t2, &t2);
         mp_sqrt(&t2, &t2);
@@ -382,10 +383,12 @@ mp_atan(const MPNumber *x, MPAngleUnit unit, MPNumber *z)
     q = 1;
     while (t2.exponent >= 0)
     {
-        if (t2.exponent == 0 && (t2.fraction[0] + 1) << 1 <= MP_BASE)
+        if (t2.exponent == 0 && 2 * (t2.fraction[0] + 1) <= MP_BASE)
             break;
 
-        q <<= 1;
+        q *= 2;
+        
+        /* t = t / (â??(t² + 1) + 1) */
         mp_multiply(&t2, &t2, z);
         mp_add_integer(z, 1, z);
         mp_sqrt(z, z);
@@ -410,7 +413,7 @@ mp_atan(const MPNumber *x, MPAngleUnit unit, MPNumber *z)
             break;
     }
 
-    /* CORRECT FOR ARGUMENT REDUCTION, AND EXIT */
+    /* CORRECT FOR ARGUMENT REDUCTION */
     mp_multiply_integer(z, q, z);
 
     /*  CHECK THAT RELATIVE ERROR LESS THAN 0.01 UNLESS EXPONENT
@@ -534,7 +537,7 @@ mp_asinh(const MPNumber *x, MPNumber *z)
 {
     MPNumber t;
 
-    /* asinh(x) = ln(x + sqrt(x^2 + 1)) */
+    /* sinhâ?»Â¹(x) = ln(x + â??(x² + 1)) */
     mp_multiply(x, x, &t);
     mp_add_integer(&t, 1, &t);
     mp_sqrt(&t, &t);
@@ -557,7 +560,7 @@ mp_acosh(const MPNumber *x, MPNumber *z)
         return;
     }
 
-    /* acosh(x) = ln(x + sqrt(x^2 - 1)) */
+    /* coshâ?»Â¹(x) = ln(x + â??(x² - 1)) */
     mp_multiply(x, x, &t);
     mp_add_integer(&t, -1, &t);
     mp_sqrt(&t, &t);
diff --git a/src/mp.c b/src/mp.c
index 92e8f40..52ad07c 100644
--- a/src/mp.c
+++ b/src/mp.c
@@ -167,7 +167,7 @@ mp_arg(const MPNumber *x, MPAngleUnit unit, MPNumber *z)
 void
 mp_conjugate(const MPNumber *x, MPNumber *z)
 {
-    //mp_set_from_mp(x, z);
+    mp_set_from_mp(x, z);
     //z->im_sign = -z->im_sign;
 }
 
@@ -537,18 +537,14 @@ mp_divide(const MPNumber *x, const MPNumber *y, MPNumber *z)
         return;
     }
 
-    /* FORM RECIPROCAL OF Y */
+    /* z = x � y�¹ */
+    /* FIXME: Set exponent to zero to avoid overflow in mp_multiply??? */
     mp_reciprocal(y, &t);
-
-    /* SET EXPONENT OF R(I2) TO ZERO TO AVOID OVERFLOW IN MP_MULTIPLY */
     ie = t.exponent;
     t.exponent = 0;
     i = t.fraction[0];
-
-    /* MULTIPLY BY X */
     mp_multiply(x, &t, z);
     mp_ext(i, z->fraction[0], z);
-
     z->exponent += ie;
 }
 
@@ -730,6 +726,9 @@ int
 mp_is_integer(const MPNumber *x)
 {
     MPNumber t1, t2, t3;
+    
+    if (mp_is_complex(x))
+        return 0;
 
     /* This fix is required for 1/3 repiprocal not being detected as an integer */
     /* Multiplication and division by 10000 is used to get around a
@@ -766,7 +765,10 @@ mp_is_integer(const MPNumber *x)
 int
 mp_is_natural(const MPNumber *x)
 {
-    return x->sign > 0 && mp_is_integer(x);
+    if (mp_is_complex(x))
+        return 0;
+    else
+        return x->sign > 0 && mp_is_integer(x);
 }
 
 
diff --git a/src/mp.h b/src/mp.h
index 1467b39..d771091 100644
--- a/src/mp.h
+++ b/src/mp.h
@@ -54,13 +54,13 @@
 typedef struct
 {
    /* Sign (+1, -1) or 0 for the value zero */
-   int sign;
+   int sign; //, im_sign;
 
    /* Exponent (to base MP_BASE) */
-   int exponent;
+   int exponent; //, im_exponent;
 
    /* Normalized fraction */
-   int fraction[MP_SIZE];
+   int fraction[MP_SIZE]; //, im_fraction[MP_SIZE];
 } MPNumber;
 
 typedef enum