[gcalctool/gcalctool-newui2] Tidy up
- From: Robert Ancell <rancell src gnome org>
- To: svn-commits-list gnome org
- Subject: [gcalctool/gcalctool-newui2] Tidy up
- Date: Mon, 13 Jul 2009 02:38:58 +0000 (UTC)
commit b19aa3358a94f70e34712cbda508f67f38581026
Author: Robert Ancell <robert ancell gmail com>
Date: Mon Jul 13 11:41:02 2009 +1000
Tidy up
src/mp-binary.c | 2 +-
src/mp-convert.c | 11 +-
src/mp-equation.h | 1 +
src/mp-internal.h | 2 -
src/mp-trigonometric.c | 552 +++++++++++++++++++-----------------------------
src/mp.c | 458 +++++++++++++++++++---------------------
src/mp.h | 59 +++---
7 files changed, 479 insertions(+), 606 deletions(-)
---
diff --git a/src/mp-binary.c b/src/mp-binary.c
index 740c86b..e3499c0 100644
--- a/src/mp-binary.c
+++ b/src/mp-binary.c
@@ -68,7 +68,7 @@ mp_is_overflow (const MPNumber *x, int wordlen)
{
MPNumber tmp1, tmp2;
mp_set_from_integer(2, &tmp1);
- mp_pwr_integer(&tmp1, wordlen, &tmp2);
+ mp_xpowy_integer(&tmp1, wordlen, &tmp2);
return mp_is_greater_than (&tmp2, x);
}
diff --git a/src/mp-convert.c b/src/mp-convert.c
index 875400a..db51dc0 100644
--- a/src/mp-convert.c
+++ b/src/mp-convert.c
@@ -217,8 +217,10 @@ mp_set_from_integer(int ix, MPNumber *z)
{
memset(z, 0, sizeof(MPNumber));
- if (ix == 0)
+ if (ix == 0) {
+ z->exponent = 1;
return;
+ }
if (ix < 0) {
ix = -ix;
@@ -230,9 +232,6 @@ mp_set_from_integer(int ix, MPNumber *z)
/* SET EXPONENT TO T */
z->exponent = MP_T;
- /* CLEAR FRACTION */
- memset(z->fraction, 0, (MP_T-1)*sizeof(int));
-
/* INSERT IX */
z->fraction[MP_T - 1] = ix;
@@ -488,7 +487,7 @@ mp_cast_to_string(const MPNumber *x, int base, int accuracy, int trim_zeroes, ch
/* Add rounding factor */
mp_set_from_integer(base, &MPbase);
- mp_pwr_integer(&MPbase, -(accuracy+1), &temp);
+ mp_xpowy_integer(&MPbase, -(accuracy+1), &temp);
mp_multiply_integer(&temp, base, &temp);
mp_divide_integer(&temp, 2, &temp);
mp_add(&number, &temp, &number);
@@ -729,7 +728,7 @@ mp_set_from_string(const char *str, MPNumber *z)
if (multiplier != 0) {
MPNumber t;
mp_set_from_integer(10, &t);
- mp_pwr_integer(&t, multiplier, &t);
+ mp_xpowy_integer(&t, multiplier, &t);
mp_multiply(z, &t, z);
}
diff --git a/src/mp-equation.h b/src/mp-equation.h
index a5203ce..41f9d9c 100644
--- a/src/mp-equation.h
+++ b/src/mp-equation.h
@@ -2,6 +2,7 @@
/* $Header$
*
* Copyright (C) 2004-2008 Sami Pietila
+ * Copyright (c) 2008-2009 Robert Ancell
*
* This program is free software; you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
diff --git a/src/mp-internal.h b/src/mp-internal.h
index 1de26c6..d47695b 100644
--- a/src/mp-internal.h
+++ b/src/mp-internal.h
@@ -45,7 +45,5 @@ 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);
-void mpexp1(const MPNumber *, MPNumber *);
-void mpmulq(const MPNumber *, int, int, MPNumber *);
#endif /* MP_INTERNAL_H */
diff --git a/src/mp-trigonometric.c b/src/mp-trigonometric.c
index c46f2c4..29e4ed2 100644
--- a/src/mp-trigonometric.c
+++ b/src/mp-trigonometric.c
@@ -32,144 +32,15 @@
* +1 IF X > I,
* 0 IF X == I,
* -1 IF X < I
- * DIMENSION OF R IN COMMON AT LEAST 2T+6
- * CHECK LEGALITY OF B, T, M AND MXR
*/
static int
mp_compare_mp_to_int(const MPNumber *x, int i)
{
- MPNumber t;
-
- /* CONVERT I TO MULTIPLE-PRECISION AND COMPARE */
+ MPNumber t;
mp_set_from_integer(i, &t);
return mp_compare_mp_to_mp(x, &t);
}
-/* COMPUTES Z = SIN(X) IF DO_SIN != 0, Z = COS(X) IF DO_SIN == 0,
- * USING TAYLOR SERIES. ASSUMES ABS(X) <= 1.
- * X AND Y ARE MP NUMBERS, IS AN INTEGER.
- * TIME IS O(M(T)T/LOG(T)). THIS COULD BE REDUCED TO
- * O(SQRT(T)M(T)) AS IN MPEXP1, BUT NOT WORTHWHILE UNLESS
- * T IS VERY LARGE. ASYMPTOTICALLY FASTER METHODS ARE
- * DESCRIBED IN THE REFERENCES GIVEN IN COMMENTS
- * TO MP_ATAN AND MPPIGL.
- * DIMENSION OF R IN CALLING PROGRAM MUST BE AT LEAST 3T+8
- * CHECK LEGALITY OF B, T, M AND MXR
- */
-static void
-mpsin1(const MPNumber *x, MPNumber *z, int do_sin)
-{
- int i, b2;
- MPNumber t1, t2;
-
- /* SIN(0) = 0, COS(0) = 1 */
- if (x->sign == 0) {
- if (do_sin == 0)
- mp_set_from_integer(1, z);
- else
- mp_set_from_integer(0, z);
- return;
- }
-
- b2 = max(MP_BASE, 64) << 1;
- mp_multiply(x, x, &t2);
- if (mp_compare_mp_to_int(&t2, 1) > 0) {
- mperr("*** ABS(X) > 1 IN CALL TO MPSIN1 ***");
- }
-
- if (do_sin == 0)
- mp_set_from_integer(1, &t1);
- else
- mp_set_from_mp(x, &t1);
-
- i = 1;
- if (do_sin != 0) {
- mp_set_from_mp(&t1, z);
- i = 2;
- }
- else
- mp_set_from_integer(0, z);
-
- /* POWER SERIES LOOP. REDUCE T IF POSSIBLE */
- for (; ; i += 2) {
- if (MP_T + t1.exponent <= 0)
- break;
-
- /* IF I*(I+1) IS NOT REPRESENTABLE AS AN INTEGER, THE FOLLOWING
- * DIVISION BY I*(I+1) HAS TO BE SPLIT UP.
- */
- mp_multiply(&t2, &t1, &t1);
- if (i > b2) {
- mp_divide_integer(&t1, -i, &t1);
- mp_divide_integer(&t1, i + 1, &t1);
- } else {
- mp_divide_integer(&t1, -i * (i + 1), &t1);
- }
- mp_add(&t1, z, z);
-
- if (t1.sign == 0)
- break;
- }
-
- if (do_sin == 0)
- mp_add_integer(z, 1, z);
-}
-
-/* COMPUTES MP Z = ARCTAN(1/N), ASSUMING INTEGER N > 1.
- * USES SERIES ARCTAN(X) = X - X**3/3 + X**5/5 - ...
- * DIMENSION OF R IN CALLING PROGRAM MUST BE
- * AT LEAST 2T+6
- * CHECK LEGALITY OF B, T, M AND MXR
- */
-static void
-mp_atan1N(int n, MPNumber *z)
-{
- int i, b2, id;
- MPNumber t2;
-
- if (n <= 1) {
- mperr("*** N <= 1 IN CALL TO MP_ATAN1N ***");
- mp_set_from_integer(0, z);
- return;
- }
-
- /* SET SUM TO X = 1/N */
- mp_set_from_fraction(1, n, z);
-
- /* SET ADDITIVE TERM TO X */
- mp_set_from_mp(z, &t2);
-
- /* ASSUME AT LEAST 16-BIT WORD. */
- b2 = max(MP_BASE, 64);
- if (n < b2)
- id = b2 * 7 * b2 / (n * n);
- else
- id = 0;
-
- /* MAIN LOOP. FIRST REDUCE T IF POSSIBLE */
- for (i = 1; ; i += 2) {
- if (MP_T + 2 + t2.exponent - z->exponent <= 1)
- break;
-
- /* IF (I+2)*N**2 IS NOT REPRESENTABLE AS AN INTEGER THE DIVISION
- * FOLLOWING HAS TO BE PERFORMED IN SEVERAL STEPS.
- */
- if (i >= id) {
- mp_multiply_fraction(&t2, -i, i + 2, &t2);
- mp_divide_integer(&t2, n, &t2);
- mp_divide_integer(&t2, n, &t2);
- }
- else {
- mp_multiply_fraction(&t2, -i, (i + 2)*n*n, &t2);
- }
-
- /* ADD TO SUM */
- mp_add(&t2, z, z);
- if (t2.sign == 0)
- break;
- }
-}
-
/* Convert x to radians */
static void
@@ -197,6 +68,14 @@ convert_to_radians(const MPNumber *x, MPAngleUnit unit, MPNumber *z)
}
}
+
+void
+mp_get_pi(MPNumber *z)
+{
+ mp_set_from_string("3.1415926535897932384626433832795028841971693993751058209749445923078164062862089986280348253421170679", z);
+}
+
+
static void
convert_from_radians(const MPNumber *x, MPAngleUnit unit, MPNumber *z)
{
@@ -222,28 +101,63 @@ convert_from_radians(const MPNumber *x, MPAngleUnit unit, MPNumber *z)
}
}
-/* SETS MP Z = PI TO THE AVAILABLE PRECISION.
- * USES PI/4 = 4.ARCTAN(1/5) - ARCTAN(1/239).
- * TIME IS O(T**2).
- * DIMENSION OF R MUST BE AT LEAST 3T+8
- * CHECK LEGALITY OF B, T, M AND MXR
+
+/* COMPUTES Z = SIN(X) IF DO_SIN != 0, Z = COS(X) IF DO_SIN == 0,
+ * USING TAYLOR SERIES. ASSUMES ABS(X) <= 1.
*/
-void
-mp_get_pi(MPNumber *z)
+static void
+mpsin1(const MPNumber *x, MPNumber *z, int do_sin)
{
- float prec;
- MPNumber t;
+ int i, b2;
+ MPNumber t1, t2;
+
+ /* sin(0) = 0, cos(0) = 1 */
+ if (x->sign == 0) {
+ if (do_sin == 0)
+ mp_set_from_integer(1, z);
+ else
+ mp_set_from_integer(0, z);
+ return;
+ }
- mp_atan1N(5, &t);
- mp_multiply_integer(&t, 4, &t);
- mp_atan1N(239, z);
- mp_subtract(&t, z, z);
- mp_multiply_integer(z, 4, z);
+ mp_multiply(x, x, &t2);
+ if (mp_compare_mp_to_int(&t2, 1) > 0) {
+ mperr("*** ABS(X) > 1 IN CALL TO MPSIN1 ***");
+ }
+
+ if (do_sin == 0) {
+ mp_set_from_integer(1, &t1);
+ mp_set_from_integer(0, z);
+ i = 1;
+ } else {
+ mp_set_from_mp(x, &t1);
+ mp_set_from_mp(&t1, z);
+ i = 2;
+ }
+
+ /* POWER SERIES LOOP. REDUCE T IF POSSIBLE */
+ b2 = max(MP_BASE, 64) << 1;
+ do {
+ if (MP_T + t1.exponent <= 0)
+ break;
- /* FOLLOWING MESSAGE MAY INDICATE THAT B**(T-1) IS TOO SMALL */
- prec = fabs(mp_cast_to_float(z) - 3.1416);
- if (prec >= 0.01)
- mperr("*** ERROR OCCURRED IN MP_GET_PI, RESULT INCORRECT ***");
+ /* IF I*(I+1) IS NOT REPRESENTABLE AS AN INTEGER, THE FOLLOWING
+ * DIVISION BY I*(I+1) HAS TO BE SPLIT UP.
+ */
+ mp_multiply(&t2, &t1, &t1);
+ if (i > b2) {
+ mp_divide_integer(&t1, -i, &t1);
+ mp_divide_integer(&t1, i + 1, &t1);
+ } else {
+ mp_divide_integer(&t1, -i * (i + 1), &t1);
+ }
+ mp_add(&t1, z, z);
+
+ i += 2;
+ } while (t1.sign != 0);
+
+ if (do_sin == 0)
+ mp_add_integer(z, 1, z);
}
@@ -258,73 +172,69 @@ mp_sin(const MPNumber *x, MPAngleUnit unit, MPNumber *z)
{
int ie, xs;
float rx = 0.0;
- MPNumber t1, t2;
+ MPNumber x_radians;
- convert_to_radians(x, unit, &t1);
-
- if (t1.sign == 0) {
+ /* sin(0) = 0 */
+ if (x->sign == 0) {
mp_set_from_integer(0, z);
return;
}
- xs = t1.sign;
- ie = abs(t1.exponent);
+ convert_to_radians(x, unit, &x_radians);
+
+ xs = x_radians.sign;
+ ie = abs(x_radians.exponent);
if (ie <= 2)
- rx = mp_cast_to_float(&t1);
+ rx = mp_cast_to_float(&x_radians);
- mp_abs(&t1, &t1);
+ mp_abs(&x_radians, &x_radians);
/* USE MPSIN1 IF ABS(X) <= 1 */
- if (mp_compare_mp_to_int(&t1, 1) <= 0)
+ if (mp_compare_mp_to_int(&x_radians, 1) <= 0)
{
- mpsin1(&t1, z, 1);
+ mpsin1(&x_radians, z, 1);
}
- /* FIND ABS(X) MODULO 2PI (IT WOULD SAVE TIME IF PI WERE
- * PRECOMPUTED AND SAVED IN COMMON).
- * FOR INCREASED ACCURACY COMPUTE PI/4 USING MP_ATAN1N
- */
+ /* FIND ABS(X) MODULO 2PI */
else {
- mp_atan1N(5, &t2);
- mp_multiply_integer(&t2, 4, &t2);
- mp_atan1N(239, z);
- mp_subtract(&t2, z, z);
- mp_divide(&t1, z, &t1);
- mp_divide_integer(&t1, 8, &t1);
- mp_fractional_component(&t1, &t1);
+ mp_get_pi(z);
+ mp_divide_integer(z, 4, z);
+ mp_divide(&x_radians, z, &x_radians);
+ mp_divide_integer(&x_radians, 8, &x_radians);
+ mp_fractional_component(&x_radians, &x_radians);
/* SUBTRACT 1/2, SAVE SIGN AND TAKE ABS */
- mp_add_fraction(&t1, -1, 2, &t1);
- xs = -xs * t1.sign;
+ mp_add_fraction(&x_radians, -1, 2, &x_radians);
+ xs = -xs * x_radians.sign;
if (xs == 0) {
mp_set_from_integer(0, z);
return;
}
- t1.sign = 1;
- mp_multiply_integer(&t1, 4, &t1);
+ x_radians.sign = 1;
+ mp_multiply_integer(&x_radians, 4, &x_radians);
/* IF NOT LESS THAN 1, SUBTRACT FROM 2 */
- if (t1.exponent > 0)
- mp_add_integer(&t1, -2, &t1);
+ if (x_radians.exponent > 0)
+ mp_add_integer(&x_radians, -2, &x_radians);
- if (t1.sign == 0) {
+ if (x_radians.sign == 0) {
mp_set_from_integer(0, z);
return;
}
- t1.sign = 1;
- mp_multiply_integer(&t1, 2, &t1);
+ x_radians.sign = 1;
+ mp_multiply_integer(&x_radians, 2, &x_radians);
/* NOW REDUCED TO FIRST QUADRANT, IF LESS THAN PI/4 USE
* POWER SERIES, ELSE COMPUTE COS OF COMPLEMENT
*/
- if (t1.exponent > 0) {
- mp_add_integer(&t1, -2, &t1);
- mp_multiply(&t1, z, &t1);
- mpsin1(&t1, z, 0);
+ if (x_radians.exponent > 0) {
+ mp_add_integer(&x_radians, -2, &x_radians);
+ mp_multiply(&x_radians, z, &x_radians);
+ mpsin1(&x_radians, z, 0);
} else {
- mp_multiply(&t1, z, &t1);
- mpsin1(&t1, z, 1);
+ mp_multiply(&x_radians, z, &x_radians);
+ mpsin1(&x_radians, z, 1);
}
}
@@ -348,9 +258,7 @@ mp_sin(const MPNumber *x, MPAngleUnit unit, MPNumber *z)
void
mp_cos(const MPNumber *xi, MPAngleUnit unit, MPNumber *z)
{
- MPNumber t;
-
- /* COS(0) = 1 */
+ /* cos(0) = 1 */
if (xi->sign == 0) {
mp_set_from_integer(1, z);
return;
@@ -358,15 +266,14 @@ mp_cos(const MPNumber *xi, MPAngleUnit unit, MPNumber *z)
convert_to_radians(xi, unit, z);
- /* SEE IF ABS(X) <= 1 */
+ /* Use power series if -1 >= x >= 1 */
mp_abs(z, z);
if (mp_compare_mp_to_int(z, 1) <= 0) {
- /* HERE ABS(X) <= 1 SO USE POWER SERIES */
mpsin1(z, z, 0);
} else {
- /* HERE ABS(X) > 1 SO USE COS(X) = SIN(PI/2 - ABS(X)),
- * COMPUTING PI/2 WITH ONE GUARD DIGIT.
- */
+ MPNumber t;
+
+ /* cos(x) = sin(Ï?/2 - |x|) */
mp_get_pi(&t);
mp_divide_integer(&t, 2, &t);
mp_subtract(&t, z, z);
@@ -380,14 +287,16 @@ mp_tan(const MPNumber *x, MPAngleUnit unit, MPNumber *z)
{
MPNumber cos_x, sin_x;
- mp_sin(x, unit, &sin_x);
+ /* Check for undefined values */
mp_cos(x, unit, &cos_x);
- /* Check if COS(x) == 0 */
if (mp_is_zero(&cos_x)) {
/* Translators: Error displayed when tangent value is undefined */
mperr(_("Tangent is infinite"));
return;
}
+
+ /* tan(x) = sin(x) / cos(x) */
+ mp_sin(x, unit, &sin_x);
mp_divide(&sin_x, &cos_x, z);
}
@@ -404,6 +313,7 @@ mp_asin(const MPNumber *x, MPAngleUnit unit, MPNumber *z)
{
MPNumber t1, t2;
+ /* asin(0) = 0 */
if (x->sign == 0) {
mp_set_from_integer(0, z);
return;
@@ -452,32 +362,32 @@ mp_asin(const MPNumber *x, MPAngleUnit unit, MPNumber *z)
void
mp_acos(const MPNumber *x, MPAngleUnit unit, MPNumber *z)
{
- MPNumber MP1, MP2;
- MPNumber MPn1, MPpi, MPy;
+ MPNumber t1, t2;
+ MPNumber MPn1, pi, MPy;
- mp_get_pi(&MPpi);
- mp_set_from_integer(1, &MP1);
+ mp_get_pi(&pi);
+ mp_set_from_integer(1, &t1);
mp_set_from_integer(-1, &MPn1);
- if (mp_is_greater_than(x, &MP1) || mp_is_less_than(x, &MPn1)) {
+ if (mp_is_greater_than(x, &t1) || mp_is_less_than(x, &MPn1)) {
mperr("Error");
mp_set_from_integer(0, z);
} else if (x->sign == 0) {
- mp_divide_integer(&MPpi, 2, z);
- } else if (mp_is_equal(x, &MP1)) {
+ mp_divide_integer(&pi, 2, z);
+ } else if (mp_is_equal(x, &t1)) {
mp_set_from_integer(0, z);
} else if (mp_is_equal(x, &MPn1)) {
- mp_set_from_mp(&MPpi, z);
+ mp_set_from_mp(&pi, z);
} else {
- mp_multiply(x, x, &MP2);
- mp_subtract(&MP1, &MP2, &MP2);
- mp_sqrt(&MP2, &MP2);
- mp_divide(&MP2, x, &MP2);
- mp_atan(&MP2, MP_RADIANS, &MPy);
+ mp_multiply(x, x, &t2);
+ mp_subtract(&t1, &t2, &t2);
+ mp_sqrt(&t2, &t2);
+ mp_divide(&t2, x, &t2);
+ mp_atan(&t2, MP_RADIANS, &MPy);
if (x->sign > 0) {
mp_set_from_mp(&MPy, z);
} else {
- mp_add(&MPy, &MPpi, z);
+ mp_add(&MPy, &pi, z);
}
}
@@ -487,7 +397,7 @@ mp_acos(const MPNumber *x, MPAngleUnit unit, MPNumber *z)
/* RETURNS Z = ARCTAN(X) FOR MP X AND Z, USING AN O(T.M(T)) METHOD
* WHICH COULD EASILY BE MODIFIED TO AN O(SQRT(T)M(T))
- * METHOD (AS IN MPEXP1). Z IS IN THE RANGE -PI/2 TO +PI/2.
+ * METHOD (AS IN MPEXP). Z IS IN THE RANGE -PI/2 TO +PI/2.
* FOR AN ASYMPTOTICALLY FASTER METHOD, SEE - FAST MULTIPLE-
* PRECISION EVALUATION OF ELEMENTARY FUNCTIONS
* (BY R. P. BRENT), J. ACM 23 (1976), 242-251,
@@ -560,91 +470,59 @@ mp_atan(const MPNumber *x, MPAngleUnit unit, MPNumber *z)
}
-/* MP precision hyperbolic arc cosine.
- *
- * 1. If (x < 1) then report DOMAIN error and return 0.
- *
- * 2. acosh(x) = log(x + sqrt(x^2 - 1))
- */
void
-mp_acosh(const MPNumber *x, MPNumber *z)
+mp_sinh(const MPNumber *x, MPNumber *z)
{
- MPNumber MP1;
+ MPNumber abs_x;
- mp_set_from_integer(1, &MP1);
- if (mp_is_less_than(x, &MP1)) {
- mperr("Error");
+ /* sinh(0) = 0 */
+ if (x->sign == 0) {
mp_set_from_integer(0, z);
- } else {
- mp_multiply(x, x, &MP1);
- mp_add_integer(&MP1, -1, &MP1);
- mp_sqrt(&MP1, &MP1);
- mp_add(x, &MP1, &MP1);
- mp_ln(&MP1, z);
+ return;
}
-}
-
-
-/* MP precision hyperbolic arc sine.
- *
- * 1. asinh(x) = log(x + sqrt(x^2 + 1))
- */
-void
-mp_asinh(const MPNumber *x, MPNumber *z)
-{
- MPNumber MP1;
-
- mp_multiply(x, x, &MP1);
- mp_add_integer(&MP1, 1, &MP1);
- mp_sqrt(&MP1, &MP1);
- mp_add(x, &MP1, &MP1);
- mp_ln(&MP1, z);
-}
+ /* WORK WITH ABS(X) */
+ mp_abs(x, &abs_x);
-/* MP precision hyperbolic arc tangent.
- *
- * 1. If (x <= -1 or x >= 1) then report a DOMAIN error and return 0.
- *
- * 2. atanh(x) = 0.5 * log((1 + x) / (1 - x))
- */
-void
-mp_atanh(const MPNumber *x, MPNumber *z)
-{
- MPNumber MP1, MP2;
- MPNumber MP3, MPn1;
-
- mp_set_from_integer(1, &MP1);
- mp_set_from_integer(-1, &MPn1);
+ /* If |x| < 1 USE MPEXP TO AVOID CANCELLATION, otherwise IF TOO LARGE MP_EPOWY GIVES ERROR MESSAGE */
+ if (abs_x.exponent <= 0) {
+ MPNumber exp_x, a, b;
+
+ /* ((e^|x| + 1) * (e^|x| - 1)) / e^|x| */
+ // FIXME: Solves to e^|x| - e^-|x|, why not lower branch always? */
+ mp_epowy(&abs_x, &exp_x);
+ mp_add_integer(&exp_x, 1, &a);
+ mp_add_integer(&exp_x, -1, &b);
+ mp_multiply(&a, &b, z);
+ mp_divide(z, &exp_x, z);
+ }
+ else {
+ MPNumber exp_x;
- if (mp_is_greater_equal(x, &MP1) || mp_is_less_equal(x, &MPn1)) {
- mperr("Error");
- mp_set_from_integer(0, z);
- } else {
- mp_add(&MP1, x, &MP2);
- mp_subtract(&MP1, x, &MP3);
- mp_divide(&MP2, &MP3, &MP3);
- mp_ln(&MP3, &MP3);
- mp_set_from_string("0.5", &MP1);
- mp_multiply(&MP1, &MP3, z);
+ /* e^|x| - e^-|x| */
+ mp_epowy(&abs_x, &exp_x);
+ mp_reciprocal(&exp_x, z);
+ mp_subtract(&exp_x, z, z);
}
+
+ /* DIVIDE BY TWO AND RESTORE SIGN */
+ mp_divide_integer(z, 2, z);
+ mp_multiply_integer(z, x->sign, z);
}
-/* RETURNS Z = COSH(X) FOR MP NUMBERS X AND Z, X NOT TOO LARGE.
- * USES MP_EPOWY, DIMENSION OF R IN COMMON AT LEAST 5T+12
- */
void
mp_cosh(const MPNumber *x, MPNumber *z)
{
MPNumber t;
- /* COSH(0) == 1 */
+ /* cosh(0) = 1 */
if (x->sign == 0) {
- mp_set_from_integer(1, z);
- return;
+ mp_set_from_integer(1, z);
+ return;
}
+ /* cosh(x) = (e^x + e^-x) / 2 */
mp_abs(x, &t);
mp_epowy(&t, &t);
mp_reciprocal(&t, z);
@@ -653,91 +531,103 @@ mp_cosh(const MPNumber *x, MPNumber *z)
}
-/* RETURNS Z = SINH(X) FOR MP NUMBERS X AND Z, X NOT TOO LARGE.
- * METHOD IS TO USE MP_EPOWY OR MPEXP1, SPACE = 5T+12
- * SAVE SIGN OF X AND CHECK FOR ZERO, SINH(0) = 0
- */
-void
-mp_sinh(const MPNumber *x, MPNumber *z)
-{
- int xs;
- MPNumber t1, t2;
-
- xs = x->sign;
- if (xs == 0) {
- mp_set_from_integer(0, z);
- return;
- }
-
- /* WORK WITH ABS(X) */
- mp_abs(x, &t2);
-
- /* HERE ABS(X) < 1 SO USE MPEXP1 TO AVOID CANCELLATION */
- if (t2.exponent <= 0) {
- mpexp1(&t2, &t1);
- mp_add_integer(&t1, 2, &t2);
- mp_multiply(&t2, &t1, z);
- mp_add_integer(&t1, 1, &t2);
- mp_divide(z, &t2, z);
- }
- /* HERE ABS(X) >= 1, IF TOO LARGE MP_EPOWY GIVES ERROR MESSAGE
- * INCREASE M TO AVOID OVERFLOW IF SINH(X) REPRESENTABLE
- */
- else {
- mp_epowy(&t2, &t2);
- mp_reciprocal(&t2, z);
- mp_subtract(&t2, z, z);
- }
-
- /* DIVIDE BY TWO AND RESTORE SIGN */
- mp_divide_integer(z, xs << 1, z);
-}
-
-
-/* RETURNS Z = TANH(X) FOR MP NUMBERS X AND Z,
- * USING MP_EPOWY OR MPEXP1, SPACE = 5T+12
- */
void
mp_tanh(const MPNumber *x, MPNumber *z)
{
float r__1;
- int xs;
MPNumber t;
- /* TANH(0) = 0 */
+ /* tanh(0) = 0 */
if (x->sign == 0) {
mp_set_from_integer(0, z);
return;
}
- /* SAVE SIGN AND WORK WITH ABS(X) */
- xs = x->sign;
mp_abs(x, &t);
/* SEE IF ABS(X) SO LARGE THAT RESULT IS +-1 */
r__1 = (float) MP_T * 0.5 * log((float) MP_BASE);
mp_set_from_float(r__1, z);
if (mp_compare_mp_to_mp(&t, z) > 0) {
- /* HERE ABS(X) IS VERY LARGE */
- mp_set_from_integer(xs, z);
+ mp_set_from_integer(x->sign, z);
return;
}
- /* HERE ABS(X) NOT SO LARGE */
+ /* If |x| >= 1/2 use ?, otherwise use ? to avoid cancellation */
+ /* |tanh(x)| = (e^|2x| - 1) / (e^|2x| + 1) */
mp_multiply_integer(&t, 2, &t);
if (t.exponent > 0) {
- /* HERE ABS(X) >= 1/2 SO USE MP_EPOWY */
mp_epowy(&t, &t);
mp_add_integer(&t, -1, z);
mp_add_integer(&t, 1, &t);
mp_divide(z, &t, z);
} else {
- /* HERE ABS(X) < 1/2, SO USE MPEXP1 TO AVOID CANCELLATION */
- mpexp1(&t, &t);
- mp_add_integer(&t, 2, z);
+ mp_epowy(&t, &t);
+ mp_add_integer(&t, 1, z);
mp_divide(&t, z, z);
}
- /* RESTORE SIGN */
- z->sign = xs * z->sign;
+ /* Restore sign */
+ z->sign = x->sign * z->sign;
+}
+
+
+void
+mp_asinh(const MPNumber *x, MPNumber *z)
+{
+ MPNumber t;
+
+ /* asinh(x) = ln(x + sqrt(x^2 + 1)) */
+ mp_multiply(x, x, &t);
+ mp_add_integer(&t, 1, &t);
+ mp_sqrt(&t, &t);
+ mp_add(x, &t, &t);
+ mp_ln(&t, z);
+}
+
+
+void
+mp_acosh(const MPNumber *x, MPNumber *z)
+{
+ MPNumber t;
+
+ /* Check x >= 1 */
+ mp_set_from_integer(1, &t);
+ if (mp_is_less_than(x, &t)) {
+ mperr("Error");
+ mp_set_from_integer(0, z);
+ return;
+ }
+
+ /* acosh(x) = ln(x + sqrt(x^2 - 1)) */
+ mp_multiply(x, x, &t);
+ mp_add_integer(&t, -1, &t);
+ mp_sqrt(&t, &t);
+ mp_add(x, &t, &t);
+ mp_ln(&t, z);
+}
+
+
+void
+mp_atanh(const MPNumber *x, MPNumber *z)
+{
+ MPNumber one, minus_one, n, d;
+
+ /* Check -1 <= x <= 1 */
+ mp_set_from_integer(1, &one);
+ mp_set_from_integer(-1, &minus_one);
+ if (mp_is_greater_equal(x, &one) || mp_is_less_equal(x, &minus_one)) {
+ mperr("Error");
+ mp_set_from_integer(0, z);
+ return;
+ }
+
+ /* atanh(x) = 0.5 * ln((1 + x) / (1 - x)) */
+ mp_add_integer(x, 1, &n);
+ mp_set_from_mp(x, &d);
+ mp_invert_sign(&d, &d);
+ mp_add_integer(&d, 1, &d);
+ mp_divide(&n, &d, z);
+ mp_ln(z, z);
+ mp_divide_integer(z, 2, z);
}
diff --git a/src/mp.c b/src/mp.c
index 94e8e74..9400f3e 100644
--- a/src/mp.c
+++ b/src/mp.c
@@ -33,6 +33,22 @@
// FIXME: Re-add overflow and underflow detection
+/* THIS ROUTINE IS CALLED WHEN AN ERROR CONDITION IS ENCOUNTERED, AND
+ * AFTER A MESSAGE HAS BEEN WRITTEN TO STDERR.
+ */
+void
+mperr(const char *format, ...)
+{
+ char text[MAXLINE];
+ va_list args;
+
+ va_start(args, format);
+ vsnprintf(text, MAXLINE, format, args);
+ va_end(args);
+ doerr(text);
+}
+
+
/* ROUTINE CALLED BY MP_DIVIDE AND MP_SQRT TO ENSURE THAT
* RESULTS ARE REPRESENTED EXACTLY IN T-2 DIGITS IF THEY
* CAN BE. X IS AN MP NUMBER, I AND J ARE INTEGERS.
@@ -77,7 +93,6 @@ mp_get_eulers(MPNumber *z)
}
-/* SETS Y = ABS(X) FOR MP NUMBERS X AND Y */
void
mp_abs(const MPNumber *x, MPNumber *z)
{
@@ -294,14 +309,6 @@ mp_add(const MPNumber *x, const MPNumber *y, MPNumber *z)
}
-/* ADDS MULTIPLE-PRECISION X TO INTEGER Y
- * GIVING MULTIPLE-PRECISION Z.
- * DIMENSION OF R IN CALLING PROGRAM MUST BE
- * AT LEAST 2T+6 (BUT Z(1) MAY BE R(T+5)).
- * DIMENSION R(6) BECAUSE RALPH COMPILER ON UNIVAC 1100 COMPUTERS
- * OBJECTS TO DIMENSION R(1).
- * CHECK LEGALITY OF B, T, M AND MXR
- */
void
mp_add_integer(const MPNumber *x, int y, MPNumber *z)
{
@@ -311,10 +318,6 @@ mp_add_integer(const MPNumber *x, int y, MPNumber *z)
}
-/* ADDS THE RATIONAL NUMBER I/J TO MP NUMBER X, MP RESULT IN Y
- * DIMENSION OF R MUST BE AT LEAST 2T+6
- * CHECK LEGALITY OF B, T, M AND MXR
- */
void
mp_add_fraction(const MPNumber *x, int i, int j, MPNumber *y)
{
@@ -324,67 +327,61 @@ mp_add_fraction(const MPNumber *x, int i, int j, MPNumber *y)
}
-/* FOR MP X AND Y, RETURNS FRACTIONAL PART OF X IN Y,
- * I.E., Y = X - INTEGER PART OF X (TRUNCATED TOWARDS 0).
- */
void
mp_fractional_component(const MPNumber *x, MPNumber *z)
{
- /* RETURN 0 IF X = 0
- OR IF EXPONENT SO LARGE THAT NO FRACTIONAL PART */
- if (x->sign == 0 || x->exponent >= MP_T) {
+ int i, shift;
+
+ /* Fractional component of zero is 0 */
+ if (x->sign == 0) {
mp_set_from_integer(0, z);
return;
}
- /* IF EXPONENT NOT POSITIVE CAN RETURN X */
+ /* All fractional */
if (x->exponent <= 0) {
mp_set_from_mp(x, z);
return;
}
-
- /* MOVE FRACTIONAL PART OF X TO ACCUMULATOR */
+
+ /* Shift fractional component */
+ shift = x->exponent;
+ for (i = shift; i < MP_SIZE && x->fraction[i] == 0; i++)
+ shift++;
z->sign = x->sign;
- z->exponent = x->exponent;
- memset(z->fraction, 0, z->exponent*sizeof(int));
- if (x != z)
- memcpy(z->fraction + z->exponent, x->fraction + x->exponent,
- (MP_T - x->exponent)*sizeof(int));
- memset(z->fraction + MP_T, 0, 4*sizeof(int));
-
- /* NORMALIZE RESULT AND RETURN */
- mp_normalize(z, 1);
+ z->exponent = x->exponent - shift;
+ for (i = 0; i < MP_SIZE; i++) {
+ if (i + shift >= MP_SIZE)
+ z->fraction[i] = 0;
+ else
+ z->fraction[i] = x->fraction[i + shift];
+ }
+ if (z->fraction[0] == 0)
+ z->sign = 0;
}
-/* RETURNS Y = INTEGER PART OF X (TRUNCATED TOWARDS 0), FOR MP X AND Y.
- * USE IF Y TOO LARGE TO BE REPRESENTABLE AS A SINGLE-PRECISION INTEGER.
- * (ELSE COULD USE MPCMI).
- * CHECK LEGALITY OF B, T, M AND MXR
- */
+
void
mp_integer_component(const MPNumber *x, MPNumber *z)
{
int i;
-
- mp_set_from_mp(x, z);
- if (z->sign == 0)
- return;
-
- /* IF EXPONENT LARGE ENOUGH RETURN Y = X */
- if (z->exponent >= MP_T) {
+
+ /* Integer component of zero = 0 */
+ if (x->sign == 0) {
+ mp_set_from_mp(x, z);
return;
}
- /* IF EXPONENT SMALL ENOUGH RETURN ZERO */
- if (z->exponent <= 0) {
+ /* If all fractional then no integer component */
+ if (x->exponent <= 0) {
mp_set_from_integer(0, z);
return;
}
- /* SET FRACTION TO ZERO */
- for (i = z->exponent; i < MP_T; i++) {
+ /* Clear fraction */
+ mp_set_from_mp(x, z);
+ for (i = z->exponent; i < MP_SIZE; i++)
z->fraction[i] = 0;
- }
}
@@ -636,6 +633,7 @@ mp_is_integer(const MPNumber *x)
{
MPNumber MPtt, MP0, MP1;
+ /* 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
* limitation to the "fix" for Sun bugtraq bug #4006391 in the
* mp_integer_component() routine in mp.c, when the exponent is less than 1.
@@ -644,8 +642,26 @@ mp_is_integer(const MPNumber *x)
mp_multiply(x, &MPtt, &MP0);
mp_divide(&MP0, &MPtt, &MP0);
mp_integer_component(&MP0, &MP1);
-
return mp_is_equal(&MP0, &MP1);
+
+ /* Correct way to check for integer */
+ /*int i;
+
+ // Zero is an integer
+ if (x->sign == 0)
+ return 1;
+
+ // Fractional
+ if (x->exponent <= 0)
+ return 0;
+
+ // Look for fractional components
+ for (i = x->exponent; i < MP_SIZE; i++) {
+ if (x->fraction[i] != 0)
+ return 0;
+ }
+
+ return 1;*/
}
@@ -656,7 +672,6 @@ mp_is_natural(const MPNumber *x)
}
-/* RETURNS LOGICAL VALUE OF (X == Y) FOR MP X AND Y. */
int
mp_is_equal(const MPNumber *x, const MPNumber *y)
{
@@ -664,25 +679,85 @@ mp_is_equal(const MPNumber *x, const MPNumber *y)
}
-/* THIS ROUTINE IS CALLED WHEN AN ERROR CONDITION IS ENCOUNTERED, AND
- * AFTER A MESSAGE HAS BEEN WRITTEN TO STDERR.
+/* Return e^x for |x| < 1 USING AN O(SQRT(T).M(T)) ALGORITHM
+ * DESCRIBED IN - R. P. BRENT, THE COMPLEXITY OF MULTIPLE-
+ * PRECISION ARITHMETIC (IN COMPLEXITY OF COMPUTATIONAL PROBLEM
+ * SOLVING, UNIV. OF QUEENSLAND PRESS, BRISBANE, 1976, 126-165).
+ * ASYMPTOTICALLY FASTER METHODS EXIST, BUT ARE NOT USEFUL
+ * UNLESS T IS VERY LARGE. SEE COMMENTS TO MP_ATAN AND MPPIGL.
*/
-void
-mperr(const char *format, ...)
+static void
+mpexp(const MPNumber *x, MPNumber *z)
{
- char text[MAXLINE];
- va_list args;
+ int i, q, ib, ic;
+ float rlb;
+ MPNumber t1, t2;
+
+ /* e^0 = 1 */
+ if (x->sign == 0) {
+ mp_set_from_integer(1, z);
+ return;
+ }
+
+ /* CHECK THAT ABS(X) < 1 */
+ if (x->exponent > 0) {
+ mperr("*** ABS(X) NOT LESS THAN 1 IN CALL TO MPEXP ***");
+ mp_set_from_integer(0, z);
+ return;
+ }
+
+ mp_set_from_mp(x, &t1);
+ rlb = log((float)MP_BASE);
+
+ /* COMPUTE APPROXIMATELY OPTIMAL Q (AND DIVIDE X BY 2**Q) */
+ q = (int)(sqrt((float)MP_T * 0.48f * rlb) + (float) x->exponent * 1.44f * rlb);
+
+ /* HALVE Q TIMES */
+ if (q > 0) {
+ ib = MP_BASE << 2;
+ ic = 1;
+ for (i = 1; i <= q; ++i) {
+ ic <<= 1;
+ if (ic < ib && ic != MP_BASE && i < q)
+ continue;
+ mp_divide_integer(&t1, ic, &t1);
+ ic = 1;
+ }
+ }
+
+ if (t1.sign == 0) {
+ mp_set_from_integer(0, z);
+ return;
+ }
+ mp_set_from_mp(&t1, z);
+ mp_set_from_mp(&t1, &t2);
+
+ /* SUM SERIES, REDUCING T WHERE POSSIBLE */
+ for (i = 2; ; i++) {
+ if (MP_T + t2.exponent - z->exponent <= 0)
+ break;
+
+ mp_multiply(&t1, &t2, &t2);
+ mp_divide_integer(&t2, i, &t2);
+
+ mp_add(&t2, z, z);
+ if (t2.sign == 0)
+ break;
+ }
+
+ /* APPLY (X+1)**2 - 1 = X(2 + X) FOR Q ITERATIONS */
+ for (i = 1; i <= q; ++i) {
+ mp_add_integer(z, 2, &t1);
+ mp_multiply(&t1, z, z);
+ }
- va_start(args, format);
- vsnprintf(text, MAXLINE, format, args);
- va_end(args);
- doerr(text);
+ mp_add_integer(z, 1, z);
}
/* RETURNS Z = EXP(X) FOR MP X AND Z.
* EXP OF INTEGER AND FRACTIONAL PARTS OF X ARE COMPUTED
- * SEPARATELY. SEE ALSO COMMENTS IN MPEXP1.
+ * SEPARATELY. SEE ALSO COMMENTS IN MPEXP.
* TIME IS O(SQRT(T)M(T)).
* DIMENSION OF R MUST BE AT LEAST 4T+10 IN CALLING PROGRAM
* CHECK LEGALITY OF B, T, M AND MXR
@@ -695,17 +770,15 @@ mp_epowy(const MPNumber *x, MPNumber *z)
float rx, rz, rlb;
MPNumber t1, t2;
- /* CHECK FOR X == 0 */
+ /* x^0 = 1 */
if (x->sign == 0) {
mp_set_from_integer(1, z);
return;
}
- /* CHECK IF ABS(X) < 1 */
+ /* If |x| < 1 use mpexp */
if (x->exponent <= 0) {
- /* USE MPEXP1 HERE */
- mpexp1(x, z);
- mp_add_integer(z, 1, z);
+ mpexp(x, z);
return;
}
@@ -734,8 +807,7 @@ mp_epowy(const MPNumber *x, MPNumber *z)
/* ATTACH SIGN TO FRACTIONAL PART AND COMPUTE EXP OF IT */
t2.sign *= xs;
- mpexp1(&t2, z);
- mp_add_integer(z, 1, z);
+ mpexp(&t2, z);
/* COMPUTE E-2 OR 1/E USING TWO EXTRA DIGITS IN CASE ABS(X) LARGE
* (BUT ONLY ONE EXTRA DIGIT IF T < 4)
@@ -763,7 +835,7 @@ mp_epowy(const MPNumber *x, MPNumber *z)
/* RAISE E OR 1/E TO POWER IX */
if (xs > 0)
mp_add_integer(&t2, 2, &t2);
- mp_pwr_integer(&t2, ix, &t2);
+ mp_xpowy_integer(&t2, ix, &t2);
/* MULTIPLY EXPS OF INTEGER AND FRACTIONAL PARTS */
mp_multiply(z, &t2, z);
@@ -797,86 +869,6 @@ mp_epowy(const MPNumber *x, MPNumber *z)
}
-/* ASSUMES THAT X AND Y ARE MP NUMBERS, -1 < X < 1.
- * RETURNS Y = EXP(X) - 1 USING AN O(SQRT(T).M(T)) ALGORITHM
- * DESCRIBED IN - R. P. BRENT, THE COMPLEXITY OF MULTIPLE-
- * PRECISION ARITHMETIC (IN COMPLEXITY OF COMPUTATIONAL PROBLEM
- * SOLVING, UNIV. OF QUEENSLAND PRESS, BRISBANE, 1976, 126-165).
- * ASYMPTOTICALLY FASTER METHODS EXIST, BUT ARE NOT USEFUL
- * UNLESS T IS VERY LARGE. SEE COMMENTS TO MP_ATAN AND MPPIGL.
- * DIMENSION OF R IN CALLING PROGRAM MUST BE AT LEAST 3T+8
- * CHECK LEGALITY OF B, T, M AND MXR
- */
-void
-mpexp1(const MPNumber *x, MPNumber *z)
-{
- int i, q, ib, ic;
- float rlb;
- MPNumber t1, t2;
-
- /* CHECK FOR X = 0 */
- if (x->sign == 0) {
- mp_set_from_integer(0, z);
- return;
- }
-
- /* CHECK THAT ABS(X) < 1 */
- if (x->exponent > 0) {
- mperr("*** ABS(X) NOT LESS THAN 1 IN CALL TO MPEXP1 ***");
- mp_set_from_integer(0, z);
- return;
- }
-
- mp_set_from_mp(x, &t1);
- rlb = log((float)MP_BASE);
-
- /* COMPUTE APPROXIMATELY OPTIMAL Q (AND DIVIDE X BY 2**Q) */
- q = (int)(sqrt((float)MP_T * 0.48f * rlb) + (float) x->exponent * 1.44f * rlb);
-
- /* HALVE Q TIMES */
- if (q > 0) {
- ib = MP_BASE << 2;
- ic = 1;
- for (i = 1; i <= q; ++i) {
- ic <<= 1;
- if (ic < ib && ic != MP_BASE && i < q)
- continue;
- mp_divide_integer(&t1, ic, &t1);
- ic = 1;
- }
- }
-
- if (t1.sign == 0) {
- mp_set_from_integer(0, z);
- return;
- }
- mp_set_from_mp(&t1, z);
- mp_set_from_mp(&t1, &t2);
-
- /* SUM SERIES, REDUCING T WHERE POSSIBLE */
- for (i = 2; ; i++) {
- if (MP_T + t2.exponent - z->exponent <= 0)
- break;
-
- mp_multiply(&t1, &t2, &t2);
- mp_divide_integer(&t2, i, &t2);
-
- mp_add(&t2, z, z);
- if (t2.sign == 0)
- break;
- }
-
- if (q <= 0)
- return;
-
- /* APPLY (X+1)**2 - 1 = X(2 + X) FOR Q ITERATIONS */
- for (i = 1; i <= q; ++i) {
- mp_add_integer(z, 2, &t1);
- mp_multiply(&t1, z, z);
- }
-}
-
-
/* RETURNS K = K/GCD AND L = L/GCD, WHERE GCD IS THE
* GREATEST COMMON DIVISOR OF K AND L.
* SAVE INPUT PARAMETERS IN LOCAL VARIABLES
@@ -958,8 +950,7 @@ mp_is_less_equal(const MPNumber *x, const MPNumber *y)
* CONDITION ABS(X) < 1/B, ERROR OTHERWISE.
* USES NEWTONS METHOD TO SOLVE THE EQUATION
* EXP1(-Y) = X, THEN REVERSES SIGN OF Y.
- * (HERE EXP1(Y) = EXP(Y) - 1 IS COMPUTED USING MPEXP1).
- * TIME IS O(SQRT(T).M(T)) AS FOR MPEXP1, SPACE = 5T+12.
+ * TIME IS O(SQRT(T).M(T)) AS FOR MPEXP, SPACE = 5T+12.
* CHECK LEGALITY OF B, T, M AND MXR
*/
static void
@@ -968,14 +959,14 @@ mplns(const MPNumber *x, MPNumber *z)
int t, it0;
MPNumber t1, t2, t3;
- /* CHECK FOR X = 0 EXACTLY */
+ /* ln(1) = 0 */
if (x->sign == 0) {
mp_set_from_integer(0, z);
return;
}
/* CHECK THAT ABS(X) < 1/B */
- if (x->exponent + 1 > 0) {
+ if (x->exponent >= 0) {
mperr("*** ABS(X) >= 1/B IN CALL TO MPLNS ***");
mp_set_from_integer(0, z);
return;
@@ -1001,7 +992,8 @@ mplns(const MPNumber *x, MPNumber *z)
{
int ts2, ts3;
- mpexp1(z, &t3);
+ mp_epowy(z, &t3);
+ mp_add_integer(&t3, -1, &t3);
mp_multiply(&t2, &t3, &t1);
mp_add(&t3, &t1, &t3);
mp_add(&t2, &t3, &t3);
@@ -1039,7 +1031,7 @@ mplns(const MPNumber *x, MPNumber *z)
* FOR SMALL INTEGER X, MPLNI IS FASTER.
* ASYMPTOTICALLY FASTER METHODS EXIST (EG THE GAUSS-SALAMIN
* METHOD, SEE MPLNGS), BUT ARE NOT USEFUL UNLESS T IS LARGE.
- * SEE COMMENTS TO MP_ATAN, MPEXP1 AND MPPIGL.
+ * SEE COMMENTS TO MP_ATAN, MPEXP AND MPPIGL.
* DIMENSION OF R IN CALLING PROGRAM MUST BE AT LEAST 6T+14.
* CHECK LEGALITY OF B, T, M AND MXR
*/
@@ -1501,95 +1493,34 @@ mp_normalize(MPNumber *x, int trunc)
}
-/* RETURNS Y = X**N, FOR MP X AND Y, INTEGER N, WITH 0**0 = 1.
- * R MUST BE DIMENSIONED AT LEAST 4T+10 IN CALLING PROGRAM
- * (2T+6 IS ENOUGH IF N NONNEGATIVE).
- */
-void
-mp_pwr_integer(const MPNumber *x, int n, MPNumber *z)
-{
- int n2, ns;
- 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_PWR_INTEGER ***");
- mp_set_from_integer(0, z);
- return;
- }
- } else if (n2 == 0) {
- /* N == 0, RETURN Y = 1. */
- mp_set_from_integer(1, 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);
-
- /* SET PRODUCT TERM TO ONE */
- 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);
- }
-}
-
-
-/* RETURNS Z = X**Y FOR MP NUMBERS X, Y AND Z, WHERE X IS
- * POSITIVE (X == 0 ALLOWED IF Y > 0). SLOWER THAN
- * MP_PWR_INTEGER AND MPQPWR, SO USE THEM IF POSSIBLE.
- * DIMENSION OF R IN COMMON AT LEAST 7T+16
- * CHECK LEGALITY OF B, T, M AND MXR
- */
-void
+static void
mp_pwr(const MPNumber *x, const MPNumber *y, MPNumber *z)
{
MPNumber t;
+ /* (-x)^y imaginary */
if (x->sign < 0) {
mperr(_("Negative X and non-integer Y not supported"));
mp_set_from_integer(0, z);
+ return;
}
- else if (x->sign == 0)
- {
- /* HERE X IS ZERO, RETURN ZERO IF Y POSITIVE, OTHERWISE ERROR */
- if (y->sign <= 0) {
- mperr("*** X ZERO AND Y NONPOSITIVE IN CALL TO MP_PWR ***");
- }
+
+ /* 0^-y illegal */
+ if (x->sign == 0 && y->sign < 0) {
+ mperr("*** X ZERO AND Y NONPOSITIVE IN CALL TO MP_PWR ***");
mp_set_from_integer(0, z);
+ return;
}
- else {
- /* USUAL CASE HERE, X POSITIVE
- * USE MPLN AND MP_EPOWY TO COMPUTE POWER
- */
- mp_ln(x, &t);
- mp_multiply(y, &t, z);
-
- /* IF X**Y IS TOO LARGE, MP_EPOWY WILL PRINT ERROR MESSAGE */
- mp_epowy(z, z);
+
+ /* 0^0 = 1 */
+ if (x->sign == 0 && y->sign == 0) {
+ mp_set_from_integer(1, z);
+ return;
}
+
+ mp_ln(x, &t);
+ mp_multiply(y, &t, z);
+ mp_epowy(z, z);
}
@@ -1723,7 +1654,7 @@ mp_root(const MPNumber *x, int n, MPNumber *z)
return;
}
- if (x->sign < 0 && np % 2 == 0) {
+ if (x->sign < 0 && np % 2 == 0) {
mperr("*** X NEGATIVE AND N EVEN IN CALL TO MP_ROOT ***");
mp_set_from_integer(0, z);
return;
@@ -1766,7 +1697,7 @@ mp_root(const MPNumber *x, int n, MPNumber *z)
while(1) {
int ts2, ts3;
- mp_pwr_integer(&t1, np, &t2);
+ mp_xpowy_integer(&t1, np, &t2);
mp_multiply(x, &t2, &t2);
mp_add_integer(&t2, -1, &t2);
mp_multiply(&t1, &t2, &t2);
@@ -1801,7 +1732,7 @@ mp_root(const MPNumber *x, int n, MPNumber *z)
}
if (n >= 0) {
- mp_pwr_integer(&t1, n - 1, &t1);
+ mp_xpowy_integer(&t1, n - 1, &t1);
mp_multiply(x, &t1, z);
return;
}
@@ -1889,11 +1820,12 @@ mp_modulus_divide(const MPNumber *x, const MPNumber *y, MPNumber *z)
return 0;
}
+
void
mp_xpowy(const MPNumber *x, const MPNumber *y, MPNumber *z)
{
if (mp_is_integer(y)) {
- mp_pwr_integer(x, mp_cast_to_int(y), z);
+ mp_xpowy_integer(x, mp_cast_to_int(y), z);
} else {
MPNumber reciprocal;
mp_reciprocal(y, &reciprocal);
@@ -1903,3 +1835,55 @@ mp_xpowy(const MPNumber *x, const MPNumber *y, MPNumber *z)
mp_pwr(x, y, z);
}
}
+
+
+void
+mp_xpowy_integer(const MPNumber *x, int n, MPNumber *z)
+{
+ int n2, ns;
+ 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);
+ 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);
+
+ /* SET PRODUCT TERM TO ONE */
+ 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);
+ }
+}
diff --git a/src/mp.h b/src/mp.h
index a58a432..1cede21 100644
--- a/src/mp.h
+++ b/src/mp.h
@@ -2,6 +2,7 @@
/* $Header$
*
* Copyright (c) 1987-2008 Sun Microsystems, Inc. All Rights Reserved.
+ * Copyright (c) 2008-2009 Robert Ancell
*
* This program is free software; you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
@@ -43,28 +44,31 @@
/* Size of the multiple precision values */
#define MP_SIZE 1000
+/* Base for numbers */
+#define MP_BASE 10000
+
/* Object for a high precision floating point number representation
*
- * x = sign * (MP.b^(exponent-1) + MP.b^(exponent-2) + ...)
+ * x = sign * (MP_BASE^(exponent-1) + MP_BASE^(exponent-2) + ...)
*/
typedef struct
{
/* Sign (+1, -1) or 0 for the value zero */
int sign;
- /* Exponent (to base MP.b) */
+ /* Exponent (to base MP_BASE) */
int exponent;
/* Normalized fraction */
int fraction[MP_SIZE];
} MPNumber;
-typedef enum { MP_RADIANS, MP_DEGREES, MP_GRADIANS } MPAngleUnit;
-
-/* Initialise the MP state. Must be called only once and before any other MP function
- * 'accuracy' is the requested accuracy required.
- */
-void mp_init(int accuracy);
+typedef enum
+{
+ MP_RADIANS,
+ MP_DEGREES,
+ MP_GRADIANS
+} MPAngleUnit;
/* Returns:
* 0 if x == y
@@ -142,10 +146,10 @@ void mp_fractional_component(const MPNumber *x, MPNumber *z);
/* Sets z = integer part of x */
void mp_integer_component(const MPNumber *x, MPNumber *z);
-/* Sets z = ln(x) */
+/* Sets z = ln x */
void mp_ln(const MPNumber *x, MPNumber *z);
-/* Sets z = log_n(x) */
+/* Sets z = log_n x */
void mp_logarithm(int n, const MPNumber *x, MPNumber *z);
/* Sets z = Ï? */
@@ -154,12 +158,6 @@ void mp_get_pi(MPNumber *z);
/* Sets z = e */
void mp_get_eulers(MPNumber *z);
-/* Sets z = x^y */
-void mp_pwr(const MPNumber *x, const MPNumber *y, MPNumber *z);
-
-/* Sets z = x^y */
-void mp_pwr_integer(const MPNumber *x, int y, MPNumber *z);
-
/* Sets z = nâ??x */
void mp_root(const MPNumber *x, int n, MPNumber *z);
@@ -169,12 +167,15 @@ void mp_sqrt(const MPNumber *x, MPNumber *z);
/* Sets z = x! */
void mp_factorial(const MPNumber *x, MPNumber *z);
-/* Sets z = x modulo y */
+/* Sets z = x mod y */
int mp_modulus_divide(const MPNumber *x, const MPNumber *y, MPNumber *z);
/* Sets z = x^y */
void mp_xpowy(const MPNumber *x, const MPNumber *y, MPNumber *z);
+/* Sets z = x^y */
+void mp_xpowy_integer(const MPNumber *x, int y, MPNumber *z);
+
/* Sets z = e^x */
void mp_epowy(const MPNumber *x, MPNumber *z);
@@ -217,40 +218,40 @@ int mp_cast_to_int(const MPNumber *x);
*/
void mp_cast_to_string(const MPNumber *x, int base, int max_digits, int trim_zeroes, char *buffer, int buffer_length);
-/* Sets z = sin(x) */
+/* Sets z = sin x */
void mp_sin(const MPNumber *x, MPAngleUnit unit, MPNumber *z);
-/* Sets z = cos(x) */
+/* Sets z = cos x */
void mp_cos(const MPNumber *x, MPAngleUnit unit, MPNumber *z);
-/* Sets z = tan(x) */
+/* Sets z = tan x */
void mp_tan(const MPNumber *x, MPAngleUnit unit, MPNumber *z);
-/* Sets z = asin(x) */
+/* Sets z = sin�¹ x */
void mp_asin(const MPNumber *x, MPAngleUnit unit, MPNumber *z);
-/* Sets z = acos(x) */
+/* Sets z = cos�¹ x */
void mp_acos(const MPNumber *x, MPAngleUnit unit, MPNumber *z);
-/* Sets z = atan(x) */
+/* Sets z = tan�¹ x */
void mp_atan(const MPNumber *x, MPAngleUnit unit, MPNumber *z);
-/* Sets z = sinh(x) */
+/* Sets z = sinh x */
void mp_sinh(const MPNumber *x, MPNumber *z);
-/* Sets z = cosh(x) */
+/* Sets z = cosh x */
void mp_cosh(const MPNumber *x, MPNumber *z);
-/* Sets z = tanh(x) */
+/* Sets z = tanh x */
void mp_tanh(const MPNumber *x, MPNumber *z);
-/* Sets z = asinh(x) */
+/* Sets z = sinh�¹ x */
void mp_asinh(const MPNumber *x, MPNumber *z);
-/* Sets z = acosh(x) */
+/* Sets z = cosh�¹ x */
void mp_acosh(const MPNumber *x, MPNumber *z);
-/* Sets z = atanh(x) */
+/* Sets z = tanh�¹ x */
void mp_atanh(const MPNumber *x, MPNumber *z);
/* Returns true if x is cannot be represented in a binary word of length 'wordlen' */
[
Date Prev][
Date Next] [
Thread Prev][
Thread Next]
[
Thread Index]
[
Date Index]
[
Author Index]
