gcalctool r2279 - trunk/gcalctool
- From: kniederk svn gnome org
- To: svn-commits-list gnome org
- Subject: gcalctool r2279 - trunk/gcalctool
- Date: Sun, 26 Oct 2008 11:24:24 +0000 (UTC)
Author: kniederk
Date: Sun Oct 26 11:24:23 2008
New Revision: 2279
URL: http://svn.gnome.org/viewvc/gcalctool?rev=2279&view=rev
Log:
Put trigonometric functions in alphabetical order, and expand function
names mprec and mpart1
Modified:
trunk/gcalctool/mp-internal.h
trunk/gcalctool/mp-trigonometric.c
trunk/gcalctool/mp.c
trunk/gcalctool/mp.h
Modified: trunk/gcalctool/mp-internal.h
==============================================================================
--- trunk/gcalctool/mp-internal.h (original)
+++ trunk/gcalctool/mp-internal.h Sun Oct 26 11:24:23 2008
@@ -35,8 +35,8 @@
void mpnzr(int, int *, int *, int);
void mpexp1(const int *, int *);
void mpmulq(int *, int, int, int *);
-void mprec(const int *, int *);
+void mp_reciprocal(const int *, int *);
void mpadd2(const int *, const int *, int *, int, int);
-void mpart1(int, int *);
+void mp_atan1N(int n, int *z);
#endif /* MP_INTERNAL_H */
Modified: trunk/gcalctool/mp-trigonometric.c
==============================================================================
--- trunk/gcalctool/mp-trigonometric.c (original)
+++ trunk/gcalctool/mp-trigonometric.c Sun Oct 26 11:24:23 2008
@@ -31,6 +31,9 @@
// FIXME: Needed for doerr
#include "calctool.h"
+
+
+
/* COMPARES MP NUMBER X WITH INTEGER I, RETURNING
* +1 IF X > I,
* 0 IF X == I,
@@ -48,8 +51,9 @@
return mp_compare_mp_to_mp(x, &MP.r[MP.t + 4]);
}
+
/* COMPUTES Z = SIN(X) IF DO_SIN != 0, Z = COS(X) IF DO_SIN == 0,
- * USING TAYLOR SERIES. ASSUMES ABS(X) >= 1.
+ * 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
@@ -126,57 +130,21 @@
mp_add_integer(z, 1, z);
}
-/* RETURNS Z = COS(X) FOR MP X AND Z, USING MP_SIN AND MPSIN1.
- * DIMENSION OF R IN COMMON AT LEAST 5T+12.
- */
-void
-mp_cos(const int *x, int *z)
-{
- int i2;
-
- /* COS(0) = 1 */
- if (x[0] == 0) {
- mp_set_from_integer(1, z);
- return;
- }
-
- /* CHECK LEGALITY OF B, T, M AND MXR */
- mpchk(5, 12);
- i2 = MP.t * 3 + 12;
-
- /* SEE IF ABS(X) <= 1 */
- mp_abs(x, 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.
- */
- ++MP.t;
- mppi(&MP.r[i2 - 1]);
- mpdivi(&MP.r[i2 - 1], 2, &MP.r[i2 - 1]);
- --MP.t;
- mp_subtract(&MP.r[i2 - 1], z, z);
- mp_sin(z, z);
- }
-}
/* MP precision arc cosine.
*
- * 1. If (x < -1.0 or x > 1.0) then report DOMAIN error and return 0.0.
+ * 1. If (x < -1 or x > 1) then report DOMAIN error and return 0.
*
- * 2. If (x = 0.0) then acos(x) = PI/2.
+ * 2. If (x = 0) then acos(x) = PI/2.
*
- * 3. If (x = 1.0) then acos(x) = 0.0
+ * 3. If (x = 1) then acos(x) = 0
*
- * 4. If (x = -1.0) then acos(x) = PI.
+ * 4. If (x = -1) then acos(x) = PI.
*
- * 5. If (0.0 < x < 1.0) then acos(x) = atan(sqrt(1-(x**2)) / x)
+ * 5. If (0 < x < 1) then acos(x) = atan(sqrt(1-x^2) / x)
*
- * 6. If (-1.0 < x < 0.0) then acos(x) = atan(sqrt(1-(x**2)) / x) + PI
+ * 6. If (-1 < x < 0) then acos(x) = atan(sqrt(1-x^2) / x) + PI
*/
-
void
mp_acos(const int *x, int *z)
{
@@ -210,6 +178,241 @@
}
}
+
+/* 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 int *x, int *z)
+{
+ int MP1[MP_SIZE];
+
+ mp_set_from_integer(1, MP1);
+ if (mp_is_less_than(x, MP1)) {
+ doerr(_("Error"));
+ mp_set_from_integer(0, z);
+ } else {
+ mpmul(x, x, MP1);
+ mp_add_integer(MP1, -1, MP1);
+ mpsqrt(MP1, MP1);
+ mp_add(x, MP1, MP1);
+ mpln(MP1, z);
+ }
+}
+
+
+/* RETURNS Z = ARCSIN(X), ASSUMING ABS(X) <= 1,
+ * FOR MP NUMBERS X AND Z.
+ * Z IS IN THE RANGE -PI/2 TO +PI/2.
+ * METHOD IS TO USE MP_ATAN, SO TIME IS O(M(T)T).
+ * DIMENSION OF R MUST BE AT LEAST 5T+12
+ * CHECK LEGALITY OF B, T, M AND MXR
+ */
+void
+mp_asin(const int *x, int *z)
+{
+ int i2, i3;
+
+ mpchk(5, 12);
+ i3 = (MP.t << 2) + 11;
+ if (x[0] == 0) {
+ z[0] = 0;
+ return;
+ }
+
+ if (x[1] <= 0) {
+ /* HERE ABS(X) < 1, SO USE ARCTAN(X/SQRT(1 - X^2)) */
+ i2 = i3 - (MP.t + 2);
+ mp_set_from_integer(1, &MP.r[i2 - 1]);
+ mp_set_from_mp(&MP.r[i2 - 1], &MP.r[i3 - 1]);
+ mp_subtract(&MP.r[i2 - 1], x, &MP.r[i2 - 1]);
+ mp_add(&MP.r[i3 - 1], x, &MP.r[i3 - 1]);
+ mpmul(&MP.r[i2 - 1], &MP.r[i3 - 1], &MP.r[i3 - 1]);
+ mproot(&MP.r[i3 - 1], -2, &MP.r[i3 - 1]);
+ mpmul(x, &MP.r[i3 - 1], z);
+ mp_atan(z, z);
+ return;
+ }
+
+ /* HERE ABS(X) >= 1. SEE IF X == +-1 */
+ mp_set_from_integer(x[0], &MP.r[i3 - 1]);
+ if (! mp_is_equal(x, &MP.r[i3 - 1])) {
+ mperr("*** ABS(X) > 1 IN CALL TO MP_ASIN ***\n");
+ }
+
+ /* X == +-1 SO RETURN +-PI/2 */
+ mppi(z);
+ mpdivi(z, MP.r[i3 - 1] << 1, z);
+}
+
+
+/* MP precision hyperbolic arc sine.
+ *
+ * 1. asinh(x) = log(x + sqrt(x^2 + 1))
+ */
+void
+mp_asinh(const int *x, int *z)
+{
+ int MP1[MP_SIZE];
+
+ mpmul(x, x, MP1);
+ mp_add_integer(MP1, 1, MP1);
+ mpsqrt(MP1, MP1);
+ mp_add(x, MP1, MP1);
+ mpln(MP1, 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.
+ * FOR AN ASYMPTOTICALLY FASTER METHOD, SEE - FAST MULTIPLE-
+ * PRECISION EVALUATION OF ELEMENTARY FUNCTIONS
+ * (BY R. P. BRENT), J. ACM 23 (1976), 242-251,
+ * AND THE COMMENTS IN MPPIGL.
+ * DIMENSION OF R IN CALLING PROGRAM MUST BE AT LEAST 5T+12
+ * CHECK LEGALITY OF B, T, M AND MXR
+ */
+void
+mp_atan(const int *x, int *z)
+{
+ int i, q, i2, i3, ts;
+ float rx = 0.0, ry;
+
+
+ mpchk(5, 12);
+ i2 = MP.t * 3 + 9;
+ i3 = i2 + MP.t + 2;
+ if (x[0] == 0) {
+ z[0] = 0;
+ return;
+ }
+
+ mp_set_from_mp(x, &MP.r[i3 - 1]);
+ if (abs(x[1]) <= 2)
+ rx = mp_cast_to_float(x);
+
+ q = 1;
+
+ /* REDUCE ARGUMENT IF NECESSARY BEFORE USING SERIES */
+ while (MP.r[i3] >= 0)
+ {
+ if (MP.r[i3] == 0 && (MP.r[i3 + 1] + 1) << 1 <= MP.b)
+ break;
+
+ q <<= 1;
+ mpmul(&MP.r[i3 - 1], &MP.r[i3 - 1], z);
+ mp_add_integer(z, 1, z);
+ mpsqrt(z, z);
+ mp_add_integer(z, 1, z);
+ mpdiv(&MP.r[i3 - 1], z, &MP.r[i3 - 1]);
+ }
+
+ /* USE POWER SERIES NOW ARGUMENT IN (-0.5, 0.5) */
+ mp_set_from_mp(&MP.r[i3 - 1], z);
+ mpmul(&MP.r[i3 - 1], &MP.r[i3 - 1], &MP.r[i2 - 1]);
+ i = 1;
+ ts = MP.t;
+
+ /* SERIES LOOP. REDUCE T IF POSSIBLE. */
+ while ( (MP.t = ts + 2 + MP.r[i3]) > 1) {
+ MP.t = min(MP.t,ts);
+ mpmul(&MP.r[i3 - 1], &MP.r[i2 - 1], &MP.r[i3 - 1]);
+ mpmulq(&MP.r[i3 - 1], -i, i + 2, &MP.r[i3 - 1]);
+ i += 2;
+ MP.t = ts;
+ mp_add(z, &MP.r[i3 - 1], z);
+ if (MP.r[i3 - 1] == 0) break;
+ }
+
+ /* RESTORE T, CORRECT FOR ARGUMENT REDUCTION, AND EXIT */
+ MP.t = ts;
+ mpmuli(z, q, z);
+
+ /* CHECK THAT RELATIVE ERROR LESS THAN 0.01 UNLESS EXPONENT
+ * OF X IS LARGE (WHEN ATAN MIGHT NOT WORK)
+ */
+ if (abs(x[1]) > 2)
+ return;
+
+ ry = mp_cast_to_float(z);
+ if (fabs(ry - atan(rx)) < fabs(ry) * (float).01)
+ return;
+
+ /* THE FOLLOWING MESSAGE MAY INDICATE THAT B**(T-1) IS TOO SMALL. */
+ mperr("*** ERROR OCCURRED IN MP_ATAN, RESULT INCORRECT ***\n");
+}
+
+
+/* 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 int *x, int *z)
+{
+ int MP1[MP_SIZE], MP2[MP_SIZE];
+ int MP3[MP_SIZE], MPn1[MP_SIZE];
+
+ mp_set_from_integer(1, MP1);
+ mp_set_from_integer(-1, MPn1);
+
+ if (mp_is_greater_equal(x, MP1) || mp_is_less_equal(x, MPn1)) {
+ doerr(_("Error"));
+ z[0] = 0;
+ } else {
+ mp_add(MP1, x, MP2);
+ mp_subtract(MP1, x, MP3);
+ mpdiv(MP2, MP3, MP3);
+ mpln(MP3, MP3);
+ MPstr_to_num("0.5", 10, MP1);
+ mpmul(MP1, MP3, z);
+ }
+}
+
+
+/* RETURNS Z = COS(X) FOR MP X AND Z, USING MP_SIN AND MPSIN1.
+ * DIMENSION OF R IN COMMON AT LEAST 5T+12.
+ */
+void
+mp_cos(const int *x, int *z)
+{
+ int i2;
+
+ /* COS(0) = 1 */
+ if (x[0] == 0) {
+ mp_set_from_integer(1, z);
+ return;
+ }
+
+ /* CHECK LEGALITY OF B, T, M AND MXR */
+ mpchk(5, 12);
+ i2 = MP.t * 3 + 12;
+
+ /* SEE IF ABS(X) <= 1 */
+ mp_abs(x, 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.
+ */
+ ++MP.t;
+ mppi(&MP.r[i2 - 1]);
+ mpdivi(&MP.r[i2 - 1], 2, &MP.r[i2 - 1]);
+ --MP.t;
+ mp_subtract(&MP.r[i2 - 1], z, z);
+ mp_sin(z, z);
+ }
+}
+
+
/* RETURNS Z = COSH(X) FOR MP NUMBERS X AND Z, X NOT TOO LARGE.
* USES MPEXP, DIMENSION OF R IN COMMON AT LEAST 5T+12
*/
@@ -234,7 +437,7 @@
*/
MP.m += 2;
mpexp(&MP.r[i2 - 1], &MP.r[i2 - 1]);
- mprec(&MP.r[i2 - 1], z);
+ mp_reciprocal(&MP.r[i2 - 1], z);
mp_add(&MP.r[i2 - 1], z, z);
/* RESTORE M. IF RESULT OVERFLOWS OR UNDERFLOWS, MPDIVI WILL
@@ -244,29 +447,6 @@
mpdivi(z, 2, 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 int *x, int *z)
-{
- int MP1[MP_SIZE];
-
- mp_set_from_integer(1, MP1);
- if (mp_is_less_than(x, MP1)) {
- doerr(_("Error"));
- mp_set_from_integer(0, z);
- } else {
- mpmul(x, x, MP1);
- mp_add_integer(MP1, -1, MP1);
- mpsqrt(MP1, MP1);
- mp_add(x, MP1, MP1);
- mpln(MP1, z);
- }
-}
/* RETURNS Z = SIN(X) FOR MP X AND Z,
* METHOD IS TO REDUCE X TO (-1, 1) AND USE MPSIN1, SO
@@ -302,13 +482,13 @@
}
/* FIND ABS(X) MODULO 2PI (IT WOULD SAVE TIME IF PI WERE
* PRECOMPUTED AND SAVED IN COMMON).
- * FOR INCREASED ACCURACY COMPUTE PI/4 USING MPART1
+ * FOR INCREASED ACCURACY COMPUTE PI/4 USING MP_ATAN1N
*/
else {
i3 = (MP.t << 1) + 7;
- mpart1(5, &MP.r[i3 - 1]);
+ mp_atan1N(5, &MP.r[i3 - 1]);
mpmuli(&MP.r[i3 - 1], 4, &MP.r[i3 - 1]);
- mpart1(239, z);
+ mp_atan1N(239, z);
mp_subtract(&MP.r[i3 - 1], z, z);
mpdiv(&MP.r[i2 - 1], z, &MP.r[i2 - 1]);
mpdivi(&MP.r[i2 - 1], 8, &MP.r[i2 - 1]);
@@ -370,49 +550,6 @@
mperr("*** ERROR OCCURRED IN MPSIN, RESULT INCORRECT ***\n");
}
-/* RETURNS Z = ARCSIN(X), ASSUMING ABS(X) <= 1,
- * FOR MP NUMBERS X AND Z.
- * Z IS IN THE RANGE -PI/2 TO +PI/2.
- * METHOD IS TO USE MP_ATAN, SO TIME IS O(M(T)T).
- * DIMENSION OF R MUST BE AT LEAST 5T+12
- * CHECK LEGALITY OF B, T, M AND MXR
- */
-void
-mp_asin(const int *x, int *z)
-{
- int i2, i3;
-
- mpchk(5, 12);
- i3 = (MP.t << 2) + 11;
- if (x[0] == 0) {
- z[0] = 0;
- return;
- }
-
- if (x[1] <= 0) {
- /* HERE ABS(X) < 1, SO USE ARCTAN(X/SQRT(1 - X**2)) */
- i2 = i3 - (MP.t + 2);
- mp_set_from_integer(1, &MP.r[i2 - 1]);
- mp_set_from_mp(&MP.r[i2 - 1], &MP.r[i3 - 1]);
- mp_subtract(&MP.r[i2 - 1], x, &MP.r[i2 - 1]);
- mp_add(&MP.r[i3 - 1], x, &MP.r[i3 - 1]);
- mpmul(&MP.r[i2 - 1], &MP.r[i3 - 1], &MP.r[i3 - 1]);
- mproot(&MP.r[i3 - 1], -2, &MP.r[i3 - 1]);
- mpmul(x, &MP.r[i3 - 1], z);
- mp_atan(z, z);
- return;
- }
-
- /* HERE ABS(X) >= 1. SEE IF X == +-1 */
- mp_set_from_integer(x[0], &MP.r[i3 - 1]);
- if (! mp_is_equal(x, &MP.r[i3 - 1])) {
- mperr("*** ABS(X) > 1 IN CALL TO MP_ASIN ***\n");
- }
-
- /* X == +-1 SO RETURN +-PI/2 */
- mppi(z);
- mpdivi(z, MP.r[i3 - 1] << 1, z);
-}
/* RETURNS Z = SINH(X) FOR MP NUMBERS X AND Z, X NOT TOO LARGE.
* METHOD IS TO USE MPEXP OR MPEXP1, SPACE = 5T+12
@@ -451,7 +588,7 @@
else {
MP.m += 2;
mpexp(&MP.r[i3 - 1], &MP.r[i3 - 1]);
- mprec(&MP.r[i3 - 1], z);
+ mp_reciprocal(&MP.r[i3 - 1], z);
mp_subtract(&MP.r[i3 - 1], z, z);
/* RESTORE M. IF RESULT OVERFLOWS OR UNDERFLOWS, MPDIVI AT
@@ -464,22 +601,6 @@
mpdivi(z, xs << 1, z);
}
-/* MP precision hyperbolic arc sine.
- *
- * 1. asinh(x) = log(x + sqrt(x**2 + 1))
- */
-
-void
-mp_asinh(const int *x, int *z)
-{
- int MP1[MP_SIZE];
-
- mpmul(x, x, MP1);
- mp_add_integer(MP1, 1, MP1);
- mpsqrt(MP1, MP1);
- mp_add(x, MP1, MP1);
- mpln(MP1, z);
-}
void
mp_tan(const int x[MP_SIZE], int z[MP_SIZE])
@@ -496,85 +617,6 @@
mpdiv(MPsin, MPcos, 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.
- * FOR AN ASYMPTOTICALLY FASTER METHOD, SEE - FAST MULTIPLE-
- * PRECISION EVALUATION OF ELEMENTARY FUNCTIONS
- * (BY R. P. BRENT), J. ACM 23 (1976), 242-251,
- * AND THE COMMENTS IN MPPIGL.
- * DIMENSION OF R IN CALLING PROGRAM MUST BE AT LEAST 5T+12
- * CHECK LEGALITY OF B, T, M AND MXR
- */
-void
-mp_atan(const int *x, int *z)
-{
- int i, q, i2, i3, ts;
- float rx = 0.0, ry;
-
-
- mpchk(5, 12);
- i2 = MP.t * 3 + 9;
- i3 = i2 + MP.t + 2;
- if (x[0] == 0) {
- z[0] = 0;
- return;
- }
-
- mp_set_from_mp(x, &MP.r[i3 - 1]);
- if (abs(x[1]) <= 2)
- rx = mp_cast_to_float(x);
-
- q = 1;
-
- /* REDUCE ARGUMENT IF NECESSARY BEFORE USING SERIES */
- while (MP.r[i3] >= 0)
- {
- if (MP.r[i3] == 0 && (MP.r[i3 + 1] + 1) << 1 <= MP.b)
- break;
-
- q <<= 1;
- mpmul(&MP.r[i3 - 1], &MP.r[i3 - 1], z);
- mp_add_integer(z, 1, z);
- mpsqrt(z, z);
- mp_add_integer(z, 1, z);
- mpdiv(&MP.r[i3 - 1], z, &MP.r[i3 - 1]);
- }
-
- /* USE POWER SERIES NOW ARGUMENT IN (-0.5, 0.5) */
- mp_set_from_mp(&MP.r[i3 - 1], z);
- mpmul(&MP.r[i3 - 1], &MP.r[i3 - 1], &MP.r[i2 - 1]);
- i = 1;
- ts = MP.t;
-
- /* SERIES LOOP. REDUCE T IF POSSIBLE. */
- while ( (MP.t = ts + 2 + MP.r[i3]) > 1) {
- MP.t = min(MP.t,ts);
- mpmul(&MP.r[i3 - 1], &MP.r[i2 - 1], &MP.r[i3 - 1]);
- mpmulq(&MP.r[i3 - 1], -i, i + 2, &MP.r[i3 - 1]);
- i += 2;
- MP.t = ts;
- mp_add(z, &MP.r[i3 - 1], z);
- if (MP.r[i3 - 1] == 0) break;
- }
-
- /* RESTORE T, CORRECT FOR ARGUMENT REDUCTION, AND EXIT */
- MP.t = ts;
- mpmuli(z, q, z);
-
- /* CHECK THAT RELATIVE ERROR LESS THAN 0.01 UNLESS EXPONENT
- * OF X IS LARGE (WHEN ATAN MIGHT NOT WORK)
- */
- if (abs(x[1]) > 2)
- return;
-
- ry = mp_cast_to_float(z);
- if (fabs(ry - atan(rx)) < fabs(ry) * (float).01)
- return;
-
- /* THE FOLLOWING MESSAGE MAY INDICATE THAT B**(T-1) IS TOO SMALL. */
- mperr("*** ERROR OCCURRED IN MP_ATAN, RESULT INCORRECT ***\n");
-}
/* RETURNS Z = TANH(X) FOR MP NUMBERS X AND Z,
* USING MPEXP OR MPEXP1, SPACE = 5T+12
@@ -627,32 +669,3 @@
/* RESTORE SIGN */
z[0] = xs * z[0];
}
-
-/* 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 int *x, int *z)
-{
- int MP1[MP_SIZE], MP2[MP_SIZE];
- int MP3[MP_SIZE], MPn1[MP_SIZE];
-
- mp_set_from_integer(1, MP1);
- mp_set_from_integer(-1, MPn1);
-
- if (mp_is_greater_equal(x, MP1) || mp_is_less_equal(x, MPn1)) {
- doerr(_("Error"));
- z[0] = 0;
- } else {
- mp_add(MP1, x, MP2);
- mp_subtract(MP1, x, MP3);
- mpdiv(MP2, MP3, MP3);
- mpln(MP3, MP3);
- MPstr_to_num("0.5", 10, MP1);
- mpmul(MP1, MP3, z);
- }
-}
Modified: trunk/gcalctool/mp.c
==============================================================================
--- trunk/gcalctool/mp.c (original)
+++ trunk/gcalctool/mp.c Sun Oct 26 11:24:23 2008
@@ -299,21 +299,21 @@
}
-/* COMPUTES MP Y = ARCTAN(1/N), ASSUMING INTEGER N > 1.
+/* 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
*/
void
-mpart1(int n, int *y)
+mp_atan1N(int n, int *z)
{
int i, b2, i2, id, ts;
mpchk(2, 6);
if (n <= 1) {
- mperr("*** N .LE. 1 IN CALL TO MPART1 ***\n");
- y[0] = 0;
+ mperr("*** N <= 1 IN CALL TO MP_ATAN1N ***\n");
+ z[0] = 0;
return;
}
@@ -321,10 +321,10 @@
ts = MP.t;
/* SET SUM TO X = 1/N */
- mp_set_from_fraction(1, n, y);
+ mp_set_from_fraction(1, n, z);
/* SET ADDITIVE TERM TO X */
- mp_set_from_mp(y, &MP.r[i2 - 1]);
+ mp_set_from_mp(z, &MP.r[i2 - 1]);
i = 1;
id = 0;
@@ -334,7 +334,7 @@
id = b2 * 7 * b2 / (n * n);
/* MAIN LOOP. FIRST REDUCE T IF POSSIBLE */
- while ((MP.t = ts + 2 + MP.r[i2] - y[1]) > 1) {
+ while ((MP.t = ts + 2 + MP.r[i2] - z[1]) > 1) {
MP.t = min(MP.t,ts);
@@ -356,7 +356,7 @@
MP.t = ts;
/* ADD TO SUM */
- mp_add(&MP.r[i2 - 1], y, y);
+ mp_add(&MP.r[i2 - 1], z, z);
if (MP.r[i2 - 1] == 0) break;
}
MP.t = ts;
@@ -378,7 +378,7 @@
if (MP.t <= 1)
mperr("*** T = %d ILLEGAL IN CALL TO MPCHK.\nPERHAPS NOT SET BEFORE CALL TO AN MP ROUTINE ***\n", MP.t);
if (MP.m <= MP.t)
- mperr("*** M .LE. T IN CALL TO MPCHK.\nPERHAPS NOT SET BEFORE CALL TO AN MP ROUTINE ***\n");
+ mperr("*** M <= T IN CALL TO MPCHK.\nPERHAPS NOT SET BEFORE CALL TO AN MP ROUTINE ***\n");
/* 8*B*B-1 SHOULD BE REPRESENTABLE, IF NOT WILL OVERFLOW
* AND MAY BECOME NEGATIVE, SO CHECK FOR THIS
@@ -495,9 +495,9 @@
}
/* COMPARES MP NUMBER X WITH REAL NUMBER RI, RETURNING
- * +1 IF X .GT. RI,
- * 0 IF X .EQ. RI,
- * -1 IF X .LT. RI
+ * +1 IF X > RI,
+ * 0 IF X == RI,
+ * -1 IF X < RI
* DIMENSION OF R IN COMMON AT LEAST 2T+6
* CHECK LEGALITY OF B, T, M AND MXR
*/
@@ -573,14 +573,14 @@
return;
}
- /* SPACE USED BY MPREC IS 4T+10 WORDS, BUT CAN OVERLAP SLIGHTLY. */
+ /* SPACE USED BY MP_RECIPROCAL IS 4T+10 WORDS, BUT CAN OVERLAP SLIGHTLY. */
i2 = MP.t * 3 + 9;
- /* INCREASE M TO AVOID OVERFLOW IN MPREC */
+ /* INCREASE M TO AVOID OVERFLOW IN MP_RECIPROCAL */
MP.m += 2;
/* FORM RECIPROCAL OF Y */
- mprec(y, &MP.r[i2 - 1]);
+ mp_reciprocal(y, &MP.r[i2 - 1]);
/* SET EXPONENT OF R(I2) TO ZERO TO AVOID OVERFLOW IN MPMUL */
ie = MP.r[i2];
@@ -864,7 +864,7 @@
xs = x[0];
mp_abs(x, &MP.r[i3 - 1]);
- /* IF ABS(X) .GT. M POSSIBLE THAT INT(X) OVERFLOWS,
+ /* IF ABS(X) > M POSSIBLE THAT INT(X) OVERFLOWS,
* SO DIVIDE BY 32.
*/
if (fabs(rx) > (float) MP.m) {
@@ -881,7 +881,7 @@
mp_add_integer(y, 1, y);
/* COMPUTE E-2 OR 1/E USING TWO EXTRA DIGITS IN CASE ABS(X) LARGE
- * (BUT ONLY ONE EXTRA DIGIT IF T .LT. 4)
+ * (BUT ONLY ONE EXTRA DIGIT IF T < 4)
*/
tss = MP.t;
ts = MP.t + 2;
@@ -1216,7 +1216,7 @@
/* RETURNS MP Y = LN(1+X) IF X IS AN MP NUMBER SATISFYING THE
- * CONDITION ABS(X) .LT. 1/B, ERROR OTHERWISE.
+ * 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).
@@ -1241,7 +1241,7 @@
/* CHECK THAT ABS(X) < 1/B */
if (x[1] + 1 > 0) {
- mperr("*** ABS(X) .GE. 1/B IN CALL TO MPLNS ***\n");
+ mperr("*** ABS(X) >= 1/B IN CALL TO MPLNS ***\n");
y[0] = 0;
return;
}
@@ -1699,7 +1699,7 @@
b2 = MP.b / 2;
if (b2 << 1 != MP.b) {
round = 0;
- /* ODD BASE, ROUND IF R(T+1)... .GT. 1/2 */
+ /* ODD BASE, ROUND IF R(T+1)... > 1/2 */
for (i = 1; i <= 4; ++i) {
it = MP.t + i;
if ((i__1 = MP.r[it - 1] - b2) < 0)
@@ -1713,7 +1713,7 @@
}
}
else {
- /* B EVEN. ROUND IF R(T+1).GE.B2 UNLESS R(T) ODD AND ALL ZEROS
+ /* B EVEN. ROUND IF R(T+1) >= B2 UNLESS R(T) ODD AND ALL ZEROS
* AFTER R(T+2).
*/
round = 1;
@@ -1807,12 +1807,12 @@
mpchk(3, 8);
-/* ALLOW SPACE FOR MPART1 */
+/* ALLOW SPACE FOR MP_ATAN1N */
i2 = (MP.t << 1) + 7;
- mpart1(5, &MP.r[i2 - 1]);
+ mp_atan1N(5, &MP.r[i2 - 1]);
mpmuli(&MP.r[i2 - 1], 4, &MP.r[i2 - 1]);
- mpart1(239, x);
+ mp_atan1N(239, x);
mp_subtract(&MP.r[i2 - 1], x, x);
mpmuli(x, 4, x);
@@ -1862,9 +1862,9 @@
/* MOVE X */
mp_set_from_mp(x, y);
- /* IF N .LT. 0 FORM RECIPROCAL */
+ /* IF N < 0 FORM RECIPROCAL */
if (n < 0)
- mprec(y, y);
+ mp_reciprocal(y, y);
mp_set_from_mp(y, &MP.r[i2 - 1]);
/* SET PRODUCT TERM TO ONE */
@@ -1885,7 +1885,7 @@
/* RETURNS Z = X**Y FOR MP NUMBERS X, Y AND Z, WHERE X IS
- * POSITIVE (X .EQ. 0 ALLOWED IF Y .GT. 0). SLOWER THAN
+ * POSITIVE (X == 0 ALLOWED IF Y > 0). SLOWER THAN
* MPPWR AND MPQPWR, SO USE THEM IF POSSIBLE.
* DIMENSION OF R IN COMMON AT LEAST 7T+16
* CHECK LEGALITY OF B, T, M AND MXR
@@ -1924,15 +1924,15 @@
}
-/* RETURNS Y = 1/X, FOR MP X AND Y.
- * MPROOT (X, -1, Y) HAS THE SAME EFFECT.
+/* RETURNS Z = 1/X, FOR MP X AND Z.
+ * MPROOT (X, -1, Z) HAS THE SAME EFFECT.
* DIMENSION OF R MUST BE AT LEAST 4*T+10 IN CALLING PROGRAM
- * (BUT Y(1) MAY BE R(3T+9)).
+ * (BUT Z(1) MAY BE R(3T+9)).
* NEWTONS METHOD IS USED, SO FINAL ONE OR TWO DIGITS MAY
* NOT BE CORRECT.
*/
void
-mprec(const int *x, int *y)
+mp_reciprocal(const int *x, int *z)
{
/* Initialized data */
static int it[9] = { 0, 8, 6, 5, 4, 4, 4, 4, 4 };
@@ -1949,8 +1949,8 @@
i2 = (MP.t << 1) + 7;
i3 = i2 + MP.t + 2;
if (x[0] == 0) {
- mperr("*** ATTEMPTED DIVISION BY ZERO IN CALL TO MPREC ***\n");
- y[0] = 0;
+ mperr("*** ATTEMPTED DIVISION BY ZERO IN CALL TO MP_RECIPROCAL ***\n");
+ z[0] = 0;
return;
}
@@ -1974,7 +1974,7 @@
/* SAVE T (NUMBER OF DIGITS) */
ts = MP.t;
- /* START WITH SMALL T TO SAVE TIME. ENSURE THAT B**(T-1) .GE. 100 */
+ /* START WITH SMALL T TO SAVE TIME. ENSURE THAT B**(T-1) >= 100 */
MP.t = 3;
if (MP.b < 10)
MP.t = it[MP.b - 1];
@@ -2015,24 +2015,24 @@
/* THE FOLLOWING MESSAGE MAY INDICATE THAT B**(T-1) IS TOO SMALL,
* OR THAT THE STARTING APPROXIMATION IS NOT ACCURATE ENOUGH.
*/
- mperr("*** ERROR OCCURRED IN MPREC, NEWTON ITERATION NOT CONVERGING PROPERLY ***\n");
+ mperr("*** ERROR OCCURRED IN MP_RECIPROCAL, NEWTON ITERATION NOT CONVERGING PROPERLY ***\n");
}
}
/* MOVE RESULT TO Y AND RETURN AFTER RESTORING T */
MP.t = ts;
- mp_set_from_mp(&MP.r[i2 - 1], y);
+ mp_set_from_mp(&MP.r[i2 - 1], z);
/* RESTORE M AND CHECK FOR OVERFLOW (UNDERFLOW IMPOSSIBLE) */
MP.m += -2;
- if (y[1] <= MP.m)
+ if (z[1] <= MP.m)
return;
- mpovfl(y, "*** OVERFLOW OCCURRED IN MPREC ***\n");
+ mpovfl(z, "*** OVERFLOW OCCURRED IN MP_RECIPROCAL ***\n");
}
-/* RETURNS Y = X**(1/N) FOR INTEGER N, ABS(N) .LE. MAX (B, 64).
+/* RETURNS Y = X**(1/N) FOR INTEGER N, ABS(N) <= MAX (B, 64).
* AND MP NUMBERS X AND Y,
* USING NEWTONS METHOD WITHOUT DIVISIONS. SPACE = 4T+10
* (BUT Y(1) MAY BE R(3T+9))
@@ -2120,7 +2120,7 @@
/* START WITH SMALL T TO SAVE TIME */
MP.t = 3;
- /* ENSURE THAT B**(T-1) .GE. 100 */
+ /* ENSURE THAT B**(T-1) >= 100 */
if (MP.b < 10)
MP.t = it[MP.b - 1];
@@ -2189,12 +2189,12 @@
* ITMAX2 IS THE DIMENSION OF ARRAYS USED FOR MP NUMBERS,
* SO AN ERROR OCCURS IF THE COMPUTED T EXCEEDS ITMAX2 - 2.
* MPSET ALSO SETS
- * MXR = MAXDR (DIMENSION OF R IN COMMON, .GE. T+4), AND
+ * MXR = MAXDR (DIMENSION OF R IN COMMON, >= T+4), AND
* M = (W-1)/4 (MAXIMUM ALLOWABLE EXPONENT),
* WHERE W IS THE LARGEST INTEGER OF THE FORM 2**K-1 WHICH IS
- * REPRESENTABLE IN THE MACHINE, K .LE. 47
+ * REPRESENTABLE IN THE MACHINE, K <= 47
* THE COMPUTED B AND T SATISFY THE CONDITIONS
- * (T-1)*LN(B)/LN(10) .GE. IDECPL AND 8*B*B-1 .LE. W .
+ * (T-1)*LN(B)/LN(10) >= IDECPL AND 8*B*B-1 <= W .
* APPROXIMATELY MINIMAL T AND MAXIMAL B SATISFYING
* THESE CONDITIONS ARE CHOSEN.
* PARAMETERS IDECPL, ITMAX2 AND MAXDR ARE INTEGERS.
@@ -2216,13 +2216,13 @@
w = 0;
k = 0;
- /* ON CYBER 76 HAVE TO FIND K .LE. 47, SO ONLY LOOP
+ /* ON CYBER 76 HAVE TO FIND K <= 47, SO ONLY LOOP
* 47 TIMES AT MOST. IF GENUINE INTEGER WORDLENGTH
* EXCEEDS 47 BITS THIS RESTRICTION CAN BE RELAXED.
*/
for (i = 1; i <= 47; ++i) {
/* INTEGER OVERFLOW WILL EVENTUALLY OCCUR HERE
- * IF WORDLENGTH .LT. 48 BITS
+ * IF WORDLENGTH < 48 BITS
*/
w2 = w + w;
wn = w2 + 1;
@@ -2239,11 +2239,11 @@
/* SET MAXIMUM EXPONENT TO (W-1)/4 */
MP.m = w / 4;
if (idecpl <= 0) {
- mperr("*** IDECPL .LE. 0 IN CALL TO MPSET ***\n");
+ mperr("*** IDECPL <= 0 IN CALL TO MPSET ***\n");
return;
}
- /* B IS THE LARGEST POWER OF 2 SUCH THAT (8*B*B-1) .LE. W */
+ /* B IS THE LARGEST POWER OF 2 SUCH THAT (8*B*B-1) <= W */
MP.b = pow_ii(2, (k - 3) / 2);
/* 2E0 BELOW ENSURES AT LEAST ONE GUARD DIGIT */
@@ -2263,7 +2263,7 @@
mpchk(1, 4);
}
-/* RETURNS Y = SQRT(X), USING SUBROUTINE MPROOT IF X .GT. 0.
+/* RETURNS Y = SQRT(X), USING SUBROUTINE MPROOT IF X > 0.
* DIMENSION OF R IN CALLING PROGRAM MUST BE AT LEAST 4T+10
* (BUT Y(1) MAY BE R(3T+9)). X AND Y ARE MP NUMBERS.
* CHECK LEGALITY OF B, T, M AND MXR
Modified: trunk/gcalctool/mp.h
==============================================================================
--- trunk/gcalctool/mp.h (original)
+++ trunk/gcalctool/mp.h Sun Oct 26 11:24:23 2008
@@ -93,18 +93,18 @@
void mp_set_from_string(const char *number, int base, int t[MP_SIZE]);
/* mp-trigonometric.c */
+void mp_acos(const int *x, int *z);
+void mp_acosh(const int *x, int *z);
+void mp_asin(const int *x, int *z);
+void mp_asinh(const int *x, int *z);
+void mp_atan(const int *x, int *z);
+void mp_atanh(const int *x, int *z);
void mp_cos(const int *x, int *z);
void mp_cosh(const int *x, int *z);
void mp_sin(const int *x, int *z);
-void mp_asin(const int *x, int *z);
void mp_sinh(const int *x, int *z);
-void mp_asinh(const int *x, int *z);
-void mp_acosh(const int *x, int *z);
-void mp_acos(const int *x, int *z);
void mp_tan(const int *x, int *z);
-void mp_atan(const int *x, int *z);
void mp_tanh(const int *x, int *z);
-void mp_atanh(const int *x, int *z);
#endif /* MP_H */
[
Date Prev][
Date Next] [
Thread Prev][
Thread Next]
[
Thread Index]
[
Date Index]
[
Author Index]
