[gcalctool] Tidy up MP code, removing global variables etc
- From: Robert Ancell <rancell src gnome org>
- To: svn-commits-list gnome org
- Subject: [gcalctool] Tidy up MP code, removing global variables etc
- Date: Mon, 13 Jul 2009 02:52:15 +0000 (UTC)
commit 481b6e1dacaa0009d6352556c190bb0d331e52a3
Author: Robert Ancell <robert ancell gmail com>
Date: Mon Jul 13 12:51:26 2009 +1000
Tidy up MP code, removing global variables etc
src/calctool.c | 1 -
src/display.c | 5 +-
src/financial.c | 8 +-
src/gtk.c | 11 +-
src/mp-binary.c | 2 +-
src/mp-convert.c | 233 +++++++----
src/mp-equation-parser.y | 14 +-
src/mp-internal.h | 23 +-
src/mp-trigonometric.c | 630 ++++++++++----------------
src/mp.c | 1139 ++++++++++++++++------------------------------
src/mp.h | 67 ++--
11 files changed, 850 insertions(+), 1283 deletions(-)
---
diff --git a/src/calctool.c b/src/calctool.c
index 1516758..948c5c9 100644
--- a/src/calctool.c
+++ b/src/calctool.c
@@ -202,7 +202,6 @@ init_state(void)
int acc, i;
acc = MAX_DIGITS + 12; /* MP internal accuracy. */
- mp_init(acc);
v->error = FALSE; /* No calculator error initially. */
v->radix = get_radix(); /* Locale specific radix string. */
diff --git a/src/display.c b/src/display.c
index 2290cb6..583a949 100644
--- a/src/display.c
+++ b/src/display.c
@@ -711,9 +711,8 @@ make_eng_sci(GCDisplay *display, char *target, int target_len, const MPNumber *M
mp_set_from_mp(&MPval, &MPmant);
mp_set_from_integer(basevals[base], &MP1base);
- mp_pwr_integer(&MP1base, 3, &MP3base);
-
- mp_pwr_integer(&MP1base, 10, &MP10base);
+ mp_xpowy_integer(&MP1base, 3, &MP3base);
+ mp_xpowy_integer(&MP1base, 10, &MP10base);
mp_set_from_integer(1, &MP1);
mp_divide(&MP1, &MP10base, &MPatmp);
diff --git a/src/financial.c b/src/financial.c
index 17c6689..1058666 100644
--- a/src/financial.c
+++ b/src/financial.c
@@ -98,7 +98,7 @@ calc_fv(MPNumber *t, MPNumber *pmt, MPNumber *pint, MPNumber *n)
MPNumber MP1, MP2, MP3, MP4;
mp_add_integer(pint, 1, &MP1);
- mp_pwr(&MP1, n, &MP2);
+ mp_xpowy(&MP1, n, &MP2);
mp_add_integer(&MP2, -1, &MP3);
mp_multiply(pmt, &MP3, &MP4);
mp_divide(&MP4, pint, t);
@@ -138,7 +138,7 @@ calc_pmt(MPNumber *t, MPNumber *prin, MPNumber *pint, MPNumber *n)
mp_add_integer(pint, 1, &MP1);
mp_multiply_integer(n, -1, &MP2);
- mp_pwr(&MP1, &MP2, &MP3);
+ mp_xpowy(&MP1, &MP2, &MP3);
mp_multiply_integer(&MP3, -1, &MP4);
mp_add_integer(&MP4, 1, &MP1);
mp_divide(pint, &MP1, &MP2);
@@ -161,7 +161,7 @@ calc_pv(MPNumber *t, MPNumber *pmt, MPNumber *pint, MPNumber *n)
mp_add_integer(pint, 1, &MP1);
mp_multiply_integer(n, -1, &MP2);
- mp_pwr(&MP1, &MP2, &MP3);
+ mp_xpowy(&MP1, &MP2, &MP3);
mp_multiply_integer(&MP3, -1, &MP4);
mp_add_integer(&MP4, 1, &MP1);
mp_divide(&MP1, pint, &MP2);
@@ -185,7 +185,7 @@ calc_rate(MPNumber *t, MPNumber *fv, MPNumber *pv, MPNumber *n)
mp_divide(fv, pv, &MP1);
mp_set_from_integer(1, &MP2);
mp_divide(&MP2, n, &MP3);
- mp_pwr(&MP1, &MP3, &MP4);
+ mp_xpowy(&MP1, &MP3, &MP4);
mp_add_integer(&MP4, -1, t);
}
diff --git a/src/gtk.c b/src/gtk.c
index 107e505..49176e6 100644
--- a/src/gtk.c
+++ b/src/gtk.c
@@ -553,12 +553,21 @@ typedef enum {
POPUP_CENTERED /* Center popup within baseframe */
} PopupLocation;
+#include <sys/time.h>
+
static void load_ui(GtkBuilder *ui, const gchar *filename)
{
GError *error = NULL;
GtkWidget *dialog;
+ struct timeval start, stop;
+ gettimeofday(&start, NULL);
gtk_builder_add_from_file(ui, filename, &error);
+ gettimeofday(&stop, NULL);
+ if (start.tv_usec > stop.tv_usec)
+ printf("t=%lu.%-6lu\n", stop.tv_sec - start.tv_sec - 1, 1000000 + stop.tv_usec - start.tv_usec);
+ else
+ printf("t=%lu.%-6lu\n", stop.tv_sec - start.tv_sec, stop.tv_usec - start.tv_usec);
if (error == NULL)
return;
@@ -569,7 +578,7 @@ static void load_ui(GtkBuilder *ui, const gchar *filename)
N_("Error loading user interface"));
gtk_message_dialog_format_secondary_text(GTK_MESSAGE_DIALOG(dialog),
/* Translators: Description in UI error dialog when unable to load the UI files. %s is replaced with the error message provided by GTK+ */
- N_("A required file is missing, please check your installation.\n\n%s"), error->message);
+ N_("A required file is missing or damaged, please check your installation.\n\n%s"), error->message);
gtk_dialog_add_buttons(GTK_DIALOG(dialog), GTK_STOCK_QUIT, GTK_RESPONSE_ACCEPT, NULL);
gtk_dialog_run(GTK_DIALOG(dialog));
diff --git a/src/mp-binary.c b/src/mp-binary.c
index e9b10b7..b97fc98 100644
--- a/src/mp-binary.c
+++ b/src/mp-binary.c
@@ -66,7 +66,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 4a3512e..d09473c 100644
--- a/src/mp-convert.c
+++ b/src/mp-convert.c
@@ -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
@@ -31,17 +32,10 @@
* SEE IF X AND Y HAVE THE SAME ADDRESS (THEY OFTEN DO)
*/
void
-mp_set_from_mp(const MPNumber *x, MPNumber *y)
+mp_set_from_mp(const MPNumber *x, MPNumber *z)
{
- /* HERE X AND Y MUST HAVE THE SAME ADDRESS */
- if (x == y)
- return;
-
- /* NO NEED TO COPY X[1],X[2],... IF X[0] == 0 */
- if (x->sign == 0)
- y->sign = 0;
- else
- memcpy (y, x, (MP.t + 2) * sizeof(int));
+ if (z != x)
+ memcpy(z, x, sizeof(MPNumber));
}
/* CONVERTS SINGLE-PRECISION NUMBER RX TO MULTIPLE-PRECISION Z.
@@ -55,8 +49,9 @@ mp_set_from_float(float rx, MPNumber *z)
int i, k, i2, ib, ie, tp;
float rb, rj;
- mpchk();
- i2 = MP.t + 4;
+ i2 = MP_T + 4;
+
+ memset(z, 0, sizeof(MPNumber));
/* CHECK SIGN */
if (rx < (float) 0.0) {
@@ -67,7 +62,7 @@ mp_set_from_float(float rx, MPNumber *z)
rj = rx;
} else {
/* IF RX = 0E0 RETURN 0 */
- z->sign = 0;
+ mp_set_from_integer(0, z);
return;
}
@@ -88,7 +83,7 @@ mp_set_from_float(float rx, MPNumber *z)
* SET EXPONENT TO 0
*/
z->exponent = 0;
- rb = (float) MP.b;
+ rb = (float) MP_BASE;
/* CONVERSION LOOP (ASSUME SINGLE-PRECISION OPS. EXACT) */
for (i = 0; i < i2; i++) {
@@ -101,7 +96,7 @@ mp_set_from_float(float rx, MPNumber *z)
mp_normalize(z, 0);
/* Computing MAX */
- ib = max(MP.b * 7 * MP.b, 32767) / 16;
+ ib = max(MP_BASE * 7 * MP_BASE, 32767) / 16;
tp = 1;
/* NOW MULTIPLY BY 16**IE */
@@ -109,7 +104,7 @@ mp_set_from_float(float rx, MPNumber *z)
k = -ie;
for (i = 1; i <= k; ++i) {
tp <<= 4;
- if (tp <= ib && tp != MP.b && i < k)
+ if (tp <= ib && tp != MP_BASE && i < k)
continue;
mp_divide_integer(z, tp, z);
tp = 1;
@@ -117,7 +112,7 @@ mp_set_from_float(float rx, MPNumber *z)
} else if (ie > 0) {
for (i = 1; i <= ie; ++i) {
tp <<= 4;
- if (tp <= ib && tp != MP.b && i < ie)
+ if (tp <= ib && tp != MP_BASE && i < ie)
continue;
mp_multiply_integer(z, tp, z);
tp = 1;
@@ -146,8 +141,9 @@ mp_set_from_double(double dx, MPNumber *z)
int i, k, i2, ib, ie, tp;
double db, dj;
- mpchk();
- i2 = MP.t + 4;
+ i2 = MP_T + 4;
+
+ memset(z, 0, sizeof(MPNumber));
/* CHECK SIGN */
if (dx < 0.) {
@@ -157,7 +153,7 @@ mp_set_from_double(double dx, MPNumber *z)
z->sign = 1;
dj = dx;
} else {
- z->sign = 0;
+ mp_set_from_integer(0, z);
return;
}
@@ -173,7 +169,7 @@ mp_set_from_double(double dx, MPNumber *z)
*/
z->exponent = 0;
- db = (double) MP.b;
+ db = (double) MP_BASE;
/* CONVERSION LOOP (ASSUME DOUBLE-PRECISION OPS. EXACT) */
for (i = 0; i < i2; i++) {
@@ -186,7 +182,7 @@ mp_set_from_double(double dx, MPNumber *z)
mp_normalize(z, 0);
/* Computing MAX */
- ib = max(MP.b * 7 * MP.b, 32767) / 16;
+ ib = max(MP_BASE * 7 * MP_BASE, 32767) / 16;
tp = 1;
/* NOW MULTIPLY BY 16**IE */
@@ -194,7 +190,7 @@ mp_set_from_double(double dx, MPNumber *z)
k = -ie;
for (i = 1; i <= k; ++i) {
tp <<= 4;
- if (tp <= ib && tp != MP.b && i < k)
+ if (tp <= ib && tp != MP_BASE && i < k)
continue;
mp_divide_integer(z, tp, z);
tp = 1;
@@ -202,7 +198,7 @@ mp_set_from_double(double dx, MPNumber *z)
} else if (ie > 0) {
for (i = 1; i <= ie; ++i) {
tp <<= 4;
- if (tp <= ib && tp != MP.b && i < ie)
+ if (tp <= ib && tp != MP_BASE && i < ie)
continue;
mp_multiply_integer(z, tp, z);
tp = 1;
@@ -219,10 +215,10 @@ mp_set_from_double(double dx, MPNumber *z)
void
mp_set_from_integer(int ix, MPNumber *z)
{
- mpchk();
+ memset(z, 0, sizeof(MPNumber));
if (ix == 0) {
- z->sign = 0;
+ z->exponent = 1;
return;
}
@@ -234,13 +230,10 @@ mp_set_from_integer(int ix, MPNumber *z)
z->sign = 1;
/* SET EXPONENT TO T */
- z->exponent = MP.t;
-
- /* CLEAR FRACTION */
- memset(z->fraction, 0, (MP.t-1)*sizeof(int));
+ z->exponent = MP_T;
/* INSERT IX */
- z->fraction[MP.t - 1] = ix;
+ z->fraction[MP_T - 1] = ix;
/* NORMALIZE BY CALLING MPMUL2 */
mpmul2(z, 1, z, 1);
@@ -248,23 +241,24 @@ mp_set_from_integer(int ix, MPNumber *z)
/* CONVERTS THE RATIONAL NUMBER I/J TO MULTIPLE PRECISION Q. */
void
-mp_set_from_fraction(int i, int j, MPNumber *q)
+mp_set_from_fraction(int i, int j, MPNumber *z)
{
mpgcd(&i, &j);
if (j == 0) {
- mperr("*** J == 0 IN CALL TO MP_SET_FROM_FRACTION ***\n");
- q->sign = 0;
- return;
+ mperr("*** J == 0 IN CALL TO MP_SET_FROM_FRACTION ***\n");
+ mp_set_from_integer(0, z);
+ return;
}
if (j < 0) {
- i = -i;
- j = -j;
+ i = -i;
+ j = -j;
}
- mp_set_from_integer(i, q);
- if (j != 1) mp_divide_integer(q, j, q);
+ mp_set_from_integer(i, z);
+ if (j != 1)
+ mp_divide_integer(z, j, z);
}
/* CONVERTS MULTIPLE-PRECISION X TO INTEGER, AND
@@ -291,8 +285,8 @@ mp_cast_to_int(const MPNumber *x)
for (i = 0; i < x2; i++) {
int izs;
izs = ret_val;
- ret_val = MP.b * ret_val;
- if (i < MP.t)
+ ret_val = MP_BASE * ret_val;
+ if (i < MP_T)
ret_val += x->fraction[i];
/* CHECK FOR SIGNS OF INTEGER OVERFLOW */
@@ -307,11 +301,11 @@ mp_cast_to_int(const MPNumber *x)
for (i = x2 - 1; i >= 0; i--) {
int j1, kx;
- j1 = j / MP.b;
+ j1 = j / MP_BASE;
kx = 0;
- if (i < MP.t)
+ if (i < MP_T)
kx = x->fraction[i];
- if (kx != j - MP.b * j1)
+ if (kx != j - MP_BASE * j1)
return 0;
j = j1;
}
@@ -366,12 +360,11 @@ mp_cast_to_float(const MPNumber *x)
int i;
float rb, rz = 0.0;
- mpchk();
if (x->sign == 0)
return 0.0;
- rb = (float) MP.b;
- for (i = 0; i < MP.t; i++) {
+ rb = (float) MP_BASE;
+ for (i = 0; i < MP_T; i++) {
rz = rb * rz + (float)x->fraction[i];
/* CHECK IF FULL SINGLE-PRECISION ACCURACY ATTAINED */
@@ -385,7 +378,7 @@ mp_cast_to_float(const MPNumber *x)
/* 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.b) + (float).5)) > (float).6) {
+ 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.
*/
@@ -433,12 +426,11 @@ mp_cast_to_double(const MPNumber *x)
int i, tm = 0;
double d__1, db, dz2, ret_val = 0.0;
- mpchk();
if (x->sign == 0)
return 0.0;
- db = (double) MP.b;
- for (i = 0; i < MP.t; i++) {
+ db = (double) MP_BASE;
+ for (i = 0; i < MP_T; i++) {
ret_val = db * ret_val + (double) x->fraction[i];
tm = i;
@@ -459,7 +451,7 @@ mp_cast_to_double(const MPNumber *x)
/* 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.b)) + .5), abs(d__1)) > .6)) {
+ ((float) MP_BASE)) + .5), abs(d__1)) > .6)) {
/* FOLLOWING MESSAGE INDICATES THAT X IS TOO LARGE OR SMALL -
* TRY USING MPCMDE INSTEAD.
*/
@@ -495,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);
@@ -556,20 +548,69 @@ mp_cast_to_string(const MPNumber *x, int base, int accuracy, int trim_zeroes, ch
static int
-char_val(char chr, int base)
+char_val(char **c, int base)
{
- int value;
- if (chr >= '0' && chr <= '9') {
- value = chr - '0';
- } else if (chr >= 'a' && chr <= 'f') {
- value = chr - 'a' + 10;
- } else if (chr >= 'A' && chr <= 'F') {
- value = chr - 'A' + 10;
+ int i, j, value, offset;
+ const char *digits[][10] = {{"Ù ", "Ù¡", "Ù¢", "Ù£", "Ù¤", "Ù¥", "Ù¦", "Ù§", "Ù¨", "Ù©"},
+ {"Û°", "Û±", "Û²", "Û³", "Û´", "Ûµ", "Û¶", "Û·", "Û¸", "Û¹"},
+ {"ß?", "ß?", "ß?", "ß?", "ß?", "ß?", "ß?", "ß?", "ß?", "ß?"},
+ {"०", "१", "२", "३", "४", "५", "६", "à¥", "८", "९"},
+ {"০", "১", "২", "৩", "৪", "৫", "৬", "à§", "৮", "৯"},
+ {"੦", "੧", "੨", "à©©", "੪", "à©«", "੬", "à©", "à©®", "੯"},
+ {"૦", "૧", "૨", "à«©", "૪", "à««", "૬", "à«", "à«®", "૯"},
+ {"à¦", "à§", "à¨", "à©", "àª", "à«", "à¬", "à", "à®", "à¯"},
+ {"௦", "௧", "௨", "௩", "௪", "௫", "௬", "à¯", "௮", "௯"},
+ {"౦", "౧", "౨", "౩", "౪", "౫", "౬", "à±", "à±®", "౯"},
+ {"೦", "೧", "೨", "೩", "೪", "೫", "೬", "à³", "à³®", "೯"},
+ {"൦", "൧", "൨", "൩", "൪", "൫", "൬", "àµ", "൮", "൯"},
+ {"�", "�", "�", "�", "�", "�", "�", "�", "�", "�"},
+ {"à»?", "à»?", "à»?", "à»?", "à»?", "à»?", "à»?", "à»?", "à»?", "à»?"},
+ {"༠", "༡", "༢", "༣", "༤", "༥", "༦", "༧", "༨", "༩"},
+ {"á??", "á??", "á??", "á??", "á??", "á??", "á??", "á??", "á??", "á??"},
+ {"á??", "á??", "á??", "á??", "á??", "á??", "á??", "á??", "á??", "á??"},
+ {"á? ", "á?¡", "á?¢", "á?£", "á?¤", "á?¥", "á?¦", "á?§", "á?¨", "á?©"},
+ {"á ?", "á ?", "á ?", "á ?", "á ?", "á ?", "á ?", "á ?", "á ?", "á ?"},
+ {"�", "�", "�", "�", "�", "�", "�", "�", "�", "�"},
+ {"�", "�", "�", "�", "�", "�", "�", "�", "�", "�"},
+ {"á?", "á?", "á?", "á?", "á?", "á?", "á?", "á?", "á?", "á?"},
+ {"᮰", "᮱", "᮲", "᮳", "᮴", "᮵", "᮶", "᮷", "᮸", "᮹"},
+ {"á±?", "á±?", "á±?", "á±?", "á±?", "á±?", "á±?", "á±?", "á±?", "á±?"},
+ {"á±?", "á±?", "á±?", "á±?", "á±?", "á±?", "á±?", "á±?", "á±?", "á±?"},
+ {"ê? ", "ê?¡", "ê?¢", "ê?£", "ê?¤", "ê?¥", "ê?¦", "ê?§", "ê?¨", "ê?©"},
+ {"�", "�", "�", "�", "�", "�", "�", "�", "�", "�"},
+ {"�", "�", "�", "�", "�", "�", "�", "�", "�", "�"},
+ {"ê©?", "ê©?", "ê©?", "ê©?", "ê©?", "ê©?", "ê©?", "ê©?", "ê©?", "ê©?"},
+ {"ð?? ", "ð??¡", "ð??¢", "ð??£", "ð??¤", "ð??¥", "ð??¦", "ð??§", "ð??¨", "ð??©"},
+ {NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL}};
+
+ if (**c >= '0' && **c <= '9') {
+ value = **c - '0';
+ offset = 1;
+ } else if (**c >= 'a' && **c <= 'f') {
+ value = **c - 'a' + 10;
+ offset = 1;
+ } else if (**c >= 'A' && **c <= 'F') {
+ value = **c - 'A' + 10;
+ offset = 1;
} else {
- return -1;
+ for (i = 0; digits[i][0]; i++) {
+ for (j = 0; j < 10; j++) {
+ if (strncmp(*c, digits[i][j], strlen(digits[i][j])) == 0)
+ break;
+ }
+ if (j != 10)
+ break;
+ }
+ if (digits[i][0] == NULL)
+ return -1;
+ value = j;
+ offset = strlen(digits[i][j]);
}
if (value >= base)
return -1;
+
+ *c += offset;
+
return value;
}
@@ -578,19 +619,37 @@ char_val(char chr, int base)
void
mp_set_from_string(const char *str, int base, MPNumber *z)
{
- int i, negate = 0;
+ int i, negate = 0, multiplier = 0;
const char *c, *end;
- const char *fractions[] = {"½", "â??", "â??", "¼", "¾", "â??", "â??", "â??", "â??", "â??", "â??", "â??", "â??", "â??", "â??", NULL};
- int numerators[] = { 1, 1, 2, 1, 3, 1, 2, 3, 4, 1, 5, 1, 3, 5, 7};
- int denominators[] = { 2, 3, 3, 4, 4, 5, 5, 5, 5, 6, 6, 8, 8, 8, 8};
+ gboolean has_fraction = FALSE;
+
+ const char *base_suffixes[] = {"â??", "â??", "â??â??", NULL};
+ int base_values[] = {2, 8, 16, 10};
+ const char *fractions[] = {"½", "â??", "â??", "¼", "¾", "â??", "â??", "â??", "â??", "â??", "â??", "â??", "â??", "â??", "â??", NULL};
+ int numerators[] = { 1, 1, 2, 1, 3, 1, 2, 3, 4, 1, 5, 1, 3, 5, 7};
+ int denominators[] = { 2, 3, 3, 4, 4, 5, 5, 5, 5, 6, 6, 8, 8, 8, 8};
+ const char *si_suffixes[] = {"T", "G", "M", "k", "d", "c", "m", "u", "µ", "n", "p", "f", NULL};
+ int si_multipliers[] = { 12, 9, 6, 3, -1, -2, -3, -6, -6, -9, -12, -15};
+ /* Find the base */
end = str;
while (*end != '\0')
end++;
+ if (base < 0) {
+ for (i = 0; base_suffixes[i] != NULL; i++) {
+ if (end - strlen(base_suffixes[i]) < str)
+ continue;
+ if (strcmp(end - strlen(base_suffixes[i]), base_suffixes[i]) == 0)
+ break;
+ }
+ base = base_values[i];
+ }
- /* Check if this is a negative number */
+ /* Check if this has a sign */
c = str;
- if (*c == '-') {
+ if (*c == '+') {
+ c++;
+ } else if (*c == '-') {
negate = 1;
c++;
} else if (strncmp(c, "â??", strlen("â??")) == 0) {
@@ -600,13 +659,12 @@ mp_set_from_string(const char *str, int base, MPNumber *z)
/* Convert integer part */
mp_set_from_integer(0, z);
- while ((i = char_val(*c, base)) >= 0) {
+ while ((i = char_val((char **)&c, base)) >= 0) {
mp_multiply_integer(z, base, z);
mp_add_integer(z, i, z);
- c++;
}
- /* Look for fraction characters */
+ /* Look for fraction characters, e.g. â?? */
for (i = 0; fractions[i] != NULL; i++) {
if (end - strlen(fractions[i]) < str)
continue;
@@ -618,20 +676,32 @@ mp_set_from_string(const char *str, int base, MPNumber *z)
mp_set_from_fraction(numerators[i], denominators[i], &fraction);
mp_add(z, &fraction, z);
}
+
+ if (*c == '.') {
+ has_fraction = TRUE;
+ c++;
+ } else {
+ for (i = 0; si_suffixes[i] != NULL; i++) {
+ if (strncmp(c, si_suffixes[i], strlen(si_suffixes[i])) == 0)
+ break;
+ }
+ if (si_suffixes[i] != NULL) {
+ has_fraction = TRUE;
+ multiplier = si_multipliers[i];
+ c += strlen(si_suffixes[i]);
+ }
+ }
/* Convert fractional part */
- if (*c == '.') {
+ if (has_fraction) {
MPNumber numerator, denominator;
-
- c++;
mp_set_from_integer(0, &numerator);
mp_set_from_integer(1, &denominator);
- while ((i = char_val(*c, base)) >= 0) {
+ while ((i = char_val((char **)&c, base)) >= 0) {
mp_multiply_integer(&denominator, base, &denominator);
mp_multiply_integer(&numerator, base, &numerator);
mp_add_integer(&numerator, i, &numerator);
- c++;
}
mp_divide(&numerator, &denominator, &numerator);
mp_add(z, &numerator, z);
@@ -657,7 +727,7 @@ mp_set_from_string(const char *str, int base, MPNumber *z)
/* Get magnitude */
mp_set_from_integer(0, &MPexponent);
- while ((i = char_val(*c, base)) >= 0) {
+ while ((i = char_val((char **)&c, base)) >= 0) {
mp_multiply_integer(&MPexponent, base, &MPexponent);
mp_add_integer(&MPexponent, i, &MPexponent);
c++;
@@ -667,13 +737,20 @@ mp_set_from_string(const char *str, int base, MPNumber *z)
}
mp_set_from_integer(base, &MPbase);
- mp_pwr(&MPbase, &MPexponent, &temp);
+ mp_xpowy(&MPbase, &MPexponent, &temp);
mp_multiply(z, &temp, z);
}
if (c != end) {
// FIXME: Error decoding
}
+
+ if (multiplier != 0) {
+ MPNumber t;
+ mp_set_from_integer(10, &t);
+ mp_xpowy_integer(&t, multiplier, &t);
+ mp_multiply(z, &t, z);
+ }
if (negate == 1)
mp_invert_sign(z, z);
diff --git a/src/mp-equation-parser.y b/src/mp-equation-parser.y
index e14b8c7..7fe28bd 100644
--- a/src/mp-equation-parser.y
+++ b/src/mp-equation-parser.y
@@ -202,7 +202,7 @@ term:
| tABS exp tABS {mp_abs(&$2, &$$);}
| 'e' '^' term {mp_epowy(&$3, &$$);}
| term '!' {mp_factorial(&$1, &$$);}
-| term tSUPNUM {mp_pwr_integer(&$1, $2, &$$);}
+| term tSUPNUM {mp_xpowy_integer(&$1, $2, &$$);}
| term '%' {mp_divide_integer(&$1, 100, &$$);}
| tNOT term %prec LNEG {
if (!mp_is_natural(&$2)) {
@@ -233,18 +233,18 @@ func:
| tSIN term %prec HIGH {mp_sin(&$2, _mp_equation_get_extra(yyscanner)->angle_units, &$$);}
| tCOS term %prec HIGH {mp_cos(&$2, _mp_equation_get_extra(yyscanner)->angle_units, &$$);}
| tTAN term %prec HIGH {mp_tan(&$2, _mp_equation_get_extra(yyscanner)->angle_units, &$$);}
-| tSIN tSUPNUM term %prec HIGH {MPNumber t; mp_sin(&$3, _mp_equation_get_extra(yyscanner)->angle_units, &t); mp_pwr_integer(&t, $2, &$$);}
-| tCOS tSUPNUM term %prec HIGH {MPNumber t; mp_cos(&$3, _mp_equation_get_extra(yyscanner)->angle_units, &t); mp_pwr_integer(&t, $2, &$$);}
-| tTAN tSUPNUM term %prec HIGH {MPNumber t; mp_tan(&$3, _mp_equation_get_extra(yyscanner)->angle_units, &t); mp_pwr_integer(&t, $2, &$$);}
+| tSIN tSUPNUM term %prec HIGH {MPNumber t; mp_sin(&$3, _mp_equation_get_extra(yyscanner)->angle_units, &t); mp_xpowy_integer(&t, $2, &$$);}
+| tCOS tSUPNUM term %prec HIGH {MPNumber t; mp_cos(&$3, _mp_equation_get_extra(yyscanner)->angle_units, &t); mp_xpowy_integer(&t, $2, &$$);}
+| tTAN tSUPNUM term %prec HIGH {MPNumber t; mp_tan(&$3, _mp_equation_get_extra(yyscanner)->angle_units, &t); mp_xpowy_integer(&t, $2, &$$);}
| tASIN term %prec HIGH {mp_asin(&$2, _mp_equation_get_extra(yyscanner)->angle_units, &$$);}
| tACOS term %prec HIGH {mp_acos(&$2, _mp_equation_get_extra(yyscanner)->angle_units, &$$);}
| tATAN term %prec HIGH {mp_atan(&$2, _mp_equation_get_extra(yyscanner)->angle_units, &$$);}
| tSINH term %prec HIGH {mp_sinh(&$2, &$$);}
| tCOSH term %prec HIGH {mp_cosh(&$2, &$$);}
| tTANH term %prec HIGH {mp_tanh(&$2, &$$);}
-| tSINH tSUPNUM term %prec HIGH {MPNumber t; mp_sinh(&$3, &t); mp_pwr_integer(&t, $2, &$$);}
-| tCOSH tSUPNUM term %prec HIGH {MPNumber t; mp_cosh(&$3, &t); mp_pwr_integer(&t, $2, &$$);}
-| tTANH tSUPNUM term %prec HIGH {MPNumber t; mp_tanh(&$3, &t); mp_pwr_integer(&t, $2, &$$);}
+| tSINH tSUPNUM term %prec HIGH {MPNumber t; mp_sinh(&$3, &t); mp_xpowy_integer(&t, $2, &$$);}
+| tCOSH tSUPNUM term %prec HIGH {MPNumber t; mp_cosh(&$3, &t); mp_xpowy_integer(&t, $2, &$$);}
+| tTANH tSUPNUM term %prec HIGH {MPNumber t; mp_tanh(&$3, &t); mp_xpowy_integer(&t, $2, &$$);}
| tASINH term %prec HIGH {mp_asinh(&$2, &$$);}
| tACOSH term %prec HIGH {mp_acosh(&$2, &$$);}
| tATANH term %prec HIGH {mp_atanh(&$2, &$$);}
diff --git a/src/mp-internal.h b/src/mp-internal.h
index 2fc231a..d47695b 100644
--- a/src/mp-internal.h
+++ b/src/mp-internal.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
@@ -32,25 +33,17 @@
#define min(a, b) ((a) <= (b) ? (a) : (b))
#define max(a, b) ((a) >= (b) ? (a) : (b))
-/* Evil global variables that must be removed */
-struct {
- /* Base */
- int b;
-
- /* Number of digits */
- // This number is temporarily changed across calls to add/subtract/multiply/divide - it should be passed to those calls
- int t;
-
- /* Min/max exponent value */
- int m;
-} MP;
+//2E0 BELOW ENSURES AT LEAST ONE GUARD DIGIT
+//MP.t = (int) ((float) (accuracy) * log((float)10.) / log((float) MP_BASE) + (float) 2.0);
+//if (MP.t > MP_SIZE) {
+// mperr("MP_SIZE TOO SMALL IN CALL TO MPSET, INCREASE MP_SIZE AND DIMENSIONS OF MP ARRAYS TO AT LEAST %d ***", MP.t);
+// MP.t = MP_SIZE;
+//}
+#define MP_T 55
void mperr(const char *format, ...) __attribute__((format(printf, 1, 2)));
-void mpchk();
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 ddfaaf7..ac46df7 100644
--- a/src/mp-trigonometric.c
+++ b/src/mp-trigonometric.c
@@ -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
@@ -31,161 +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;
-
- mpchk();
-
- /* 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;
-
- mpchk();
-
- /* SIN(0) = 0, COS(0) = 1 */
- if (x->sign == 0) {
- z->sign = 0;
- if (do_sin == 0)
- mp_set_from_integer(1, z);
- return;
- }
-
- b2 = max(MP.b, 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);
-
- z->sign = 0;
- i = 1;
- if (do_sin != 0) {
- mp_set_from_mp(&t1, z);
- i = 2;
- }
-
- /* POWER SERIES LOOP. REDUCE T IF POSSIBLE */
- for (; ; i+= 2) {
- int t, ts;
-
- t = MP.t + t1.exponent + 2;
- if (t <= 2)
- break;
- t = min(t, MP.t);
-
- /* IF I*(I+1) IS NOT REPRESENTABLE AS AN INTEGER, THE FOLLOWING
- * DIVISION BY I*(I+1) HAS TO BE SPLIT UP.
- */
- ts = MP.t;
- MP.t = t;
- 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.t = ts;
- 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;
-
- mpchk();
- if (n <= 1) {
- mperr("*** N <= 1 IN CALL TO MP_ATAN1N ***");
- z->sign = 0;
- 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.b, 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) {
- int t, ts;
-
- t = MP.t + 2 + t2.exponent - z->exponent;
- if (t <= 1)
- break;
- t = min(t, MP.t);
-
- /* IF (I+2)*N**2 IS NOT REPRESENTABLE AS AN INTEGER THE DIVISION
- * FOLLOWING HAS TO BE PERFORMED IN SEVERAL STEPS.
- */
- ts = MP.t;
- MP.t = t;
- 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);
- }
- MP.t = ts;
-
- /* ADD TO SUM */
- mp_add(&t2, z, z);
- if (t2.sign == 0)
- break;
- }
-}
-
/* Convert x to radians */
static void
@@ -213,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", 10, z);
+}
+
+
static void
convert_from_radians(const MPNumber *x, MPAngleUnit unit, MPNumber *z)
{
@@ -238,30 +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_multiply(x, x, &t2);
+ if (mp_compare_mp_to_int(&t2, 1) > 0) {
+ mperr("*** ABS(X) > 1 IN CALL TO MPSIN1 ***");
+ }
- mpchk();
+ 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;
+ }
- 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);
+ /* 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);
}
@@ -276,75 +172,69 @@ mp_sin(const MPNumber *x, MPAngleUnit unit, MPNumber *z)
{
int ie, xs;
float rx = 0.0;
- MPNumber t1, t2;
+ MPNumber x_radians;
- mpchk();
-
- convert_to_radians(x, unit, &t1);
-
- if (t1.sign == 0) {
- z->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) {
- z->sign = 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) {
- z->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);
}
}
@@ -368,28 +258,22 @@ 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;
}
- /* CHECK LEGALITY OF B, T, M AND MXR */
- mpchk();
-
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);
@@ -403,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);
}
@@ -427,9 +313,9 @@ mp_asin(const MPNumber *x, MPAngleUnit unit, MPNumber *z)
{
MPNumber t1, t2;
- mpchk();
+ /* asin(0) = 0 */
if (x->sign == 0) {
- z->sign = 0;
+ mp_set_from_integer(0, z);
return;
}
@@ -476,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");
- z->sign = 0;
+ mp_set_from_integer(0, z);
} else if (x->sign == 0) {
- mp_divide_integer(&MPpi, 2, z);
- } else if (mp_is_equal(x, &MP1)) {
- z->sign = 0;
+ 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);
}
}
@@ -511,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,
@@ -526,9 +412,8 @@ mp_atan(const MPNumber *x, MPAngleUnit unit, MPNumber *z)
float rx = 0.0;
MPNumber t1, t2;
- mpchk();
if (x->sign == 0) {
- z->sign = 0;
+ mp_set_from_integer(0, z);
return;
}
@@ -540,7 +425,7 @@ 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.b)
+ if (t2.exponent == 0 && (t2.fraction[0] + 1) << 1 <= MP_BASE)
break;
q <<= 1;
@@ -557,18 +442,11 @@ mp_atan(const MPNumber *x, MPAngleUnit unit, MPNumber *z)
/* SERIES LOOP. REDUCE T IF POSSIBLE. */
for (i = 1; ; i += 2) {
- int t, ts;
-
- t = MP.t + 2 + t2.exponent;
- if (t <= 1)
+ if (MP_T + 2 + t2.exponent <= 1)
break;
- t = min(t, MP.t);
- ts = MP.t;
- MP.t = t;
mp_multiply(&t2, &t1, &t2);
mp_multiply_fraction(&t2, -i, i + 2, &t2);
- MP.t = ts;
mp_add(z, &t2, z);
if (t2.sign == 0)
@@ -592,208 +470,164 @@ 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");
- z->sign = 0;
- } 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", 10, &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;
}
- /* CHECK LEGALITY OF B, T, M AND MXR */
- mpchk();
+ /* cosh(x) = (e^x + e^-x) / 2 */
mp_abs(x, &t);
-
- /* IF ABS(X) TOO LARGE MP_EPOWY WILL PRINT ERROR MESSAGE
- * INCREASE M TO AVOID OVERFLOW WHEN COSH(X) REPRESENTABLE
- */
- MP.m += 2;
mp_epowy(&t, &t);
mp_reciprocal(&t, z);
mp_add(&t, z, z);
-
- /* RESTORE M. IF RESULT OVERFLOWS OR UNDERFLOWS, mp_divide_integer WILL
- * ACT ACCORDINGLY.
- */
- MP.m += -2;
mp_divide_integer(z, 2, 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) {
- z->sign = 0;
- return;
- }
-
- /* CHECK LEGALITY OF B, T, M AND MXR */
- mpchk();
-
- /* 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.m += 2;
- mp_epowy(&t2, &t2);
- mp_reciprocal(&t2, z);
- mp_subtract(&t2, z, z);
-
- /* RESTORE M. IF RESULT OVERFLOWS OR UNDERFLOWS, mp_divide_integer AT
- * STATEMENT 30 WILL ACT ACCORDINGLY.
- */
- MP.m += -2;
- }
-
- /* 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) {
- z->sign = 0;
+ mp_set_from_integer(0, z);
return;
}
- /* CHECK LEGALITY OF B, T, M AND MXR */
- mpchk();
-
- /* 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.b);
+ 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 5689ddb..9400f3e 100644
--- a/src/mp.c
+++ b/src/mp.c
@@ -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
@@ -28,90 +29,23 @@
#include "mp-internal.h"
#include "calctool.h" // FIXME: Required for doerr() and MAXLINE
-// FIXME: MP.t and MP.m modified inside functions, needs to be fixed to be thread safe
-/* SETS X TO THE LARGEST POSSIBLE POSITIVE MP NUMBER
- * CHECK LEGALITY OF B, T, M AND MXR
- */
-static void
-mpmaxr(MPNumber *x)
-{
- int i, it;
-
- mpchk();
- it = MP.b - 1;
-
- /* SET FRACTION DIGITS TO B-1 */
- for (i = 0; i < MP.t; i++)
- x->fraction[i] = it;
-
- /* SET SIGN AND EXPONENT */
- x->sign = 1;
- x->exponent = MP.m;
-}
-
-
-/* CALLED ON MULTIPLE-PRECISION OVERFLOW, IE WHEN THE
- * EXPONENT OF MP NUMBER X WOULD EXCEED M.
- * AT PRESENT EXECUTION IS TERMINATED WITH AN ERROR MESSAGE
- * AFTER CALLING MPMAXR(X), BUT IT WOULD BE POSSIBLE TO RETURN,
- * POSSIBLY UPDATING A COUNTER AND TERMINATING EXECUTION AFTER
- * A PRESET NUMBER OF OVERFLOWS. ACTION COULD EASILY BE DETERMINED
- * BY A FLAG IN LABELLED COMMON.
- * M MAY HAVE BEEN OVERWRITTEN, SO CHECK B, T, M ETC.
- */
-static void
-mpovfl(MPNumber *x, const char *error)
-{
- fprintf(stderr, "%s\n", error);
-
- mpchk();
-
- /* SET X TO LARGEST POSSIBLE POSITIVE NUMBER */
- mpmaxr(x);
-
- /* TERMINATE EXECUTION BY CALLING MPERR */
- mperr("*** CALL TO MPOVFL, MP OVERFLOW OCCURRED ***");
-}
+// FIXME: Re-add overflow and underflow detection
-/* CALLED ON MULTIPLE-PRECISION UNDERFLOW, IE WHEN THE
- * EXPONENT OF MP NUMBER X WOULD BE LESS THAN -M.
- * SINCE M MAY HAVE BEEN OVERWRITTEN, CHECK B, T, M ETC.
+/* THIS ROUTINE IS CALLED WHEN AN ERROR CONDITION IS ENCOUNTERED, AND
+ * AFTER A MESSAGE HAS BEEN WRITTEN TO STDERR.
*/
-static void
-mpunfl(MPNumber *x)
-{
- mpchk();
-
- /* THE UNDERFLOWING NUMBER IS SET TO ZERO
- * AN ALTERNATIVE WOULD BE TO CALL MPMINR (X) AND RETURN,
- * POSSIBLY UPDATING A COUNTER AND TERMINATING EXECUTION
- * AFTER A PRESET NUMBER OF UNDERFLOWS. ACTION COULD EASILY
- * BE DETERMINED BY A FLAG IN LABELLED COMMON.
- */
- x->sign = 0;
-}
-
-
-static int
-pow_ii(int x, int n)
+void
+mperr(const char *format, ...)
{
- int p = 1;
+ char text[MAXLINE];
+ va_list args;
- if (n <= 0)
- return 1;
-
- for (;;) {
- if (n & 01)
- p *= x;
- if (n >>= 1)
- x *= x;
- else
- break;
- }
-
- return p;
+ va_start(args, format);
+ vsnprintf(text, MAXLINE, format, args);
+ va_end(args);
+ doerr(text);
}
@@ -124,132 +58,41 @@ mpext(int i, int j, MPNumber *x)
{
int q, s;
- if (x->sign == 0 || MP.t <= 2 || i == 0)
+ if (x->sign == 0 || MP_T <= 2 || i == 0)
return;
/* COMPUTE MAXIMUM POSSIBLE ERROR IN THE LAST PLACE */
q = (j + 1) / i + 1;
- s = MP.b * x->fraction[MP.t - 2] + x->fraction[MP.t - 1];
+ s = MP_BASE * x->fraction[MP_T - 2] + x->fraction[MP_T - 1];
/* SET LAST TWO DIGITS TO ZERO */
if (s <= q) {
- x->fraction[MP.t - 2] = 0;
- x->fraction[MP.t - 1] = 0;
+ x->fraction[MP_T - 2] = 0;
+ x->fraction[MP_T - 1] = 0;
return;
}
- if (s + q < MP.b * MP.b)
+ if (s + q < MP_BASE * MP_BASE)
return;
/* ROUND UP HERE */
- x->fraction[MP.t - 2] = MP.b - 1;
- x->fraction[MP.t - 1] = MP.b;
+ x->fraction[MP_T - 2] = MP_BASE - 1;
+ x->fraction[MP_T - 1] = MP_BASE;
/* NORMALIZE X (LAST DIGIT B IS OK IN MP_MULTIPLY_INTEGER) */
mp_multiply_integer(x, 1, x);
}
-/* SETS BASE (B) AND NUMBER OF DIGITS (T) TO GIVE THE
- * EQUIVALENT OF AT LEAST IDECPL DECIMAL DIGITS.
- * IDECPL SHOULD BE POSITIVE.
- * 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, >= 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 <= 47
- * THE COMPUTED B AND T SATISFY THE CONDITIONS
- * (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.
- * BEWARE - MPSET WILL CAUSE AN INTEGER OVERFLOW TO OCCUR
- * ****** IF WORDLENGTH IS LESS THAN 48 BITS.
- * IF THIS IS NOT ALLOWABLE, CHANGE THE DETERMINATION
- * OF W (DO 30 ... TO 30 W = WN) OR SET B, T, M,
- * AND MXR WITHOUT CALLING MPSET.
- * FIRST SET MXR
- */
void
-mp_init(int accuracy)
+mp_get_eulers(MPNumber *z)
{
- int i, k, w, b;
-
- /* DETERMINE LARGE REPRESENTABLE INTEGER W OF FORM 2**K - 1 */
- /* 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.
- */
- w = 0;
- k = 0;
- for (i = 1; i <= 47; ++i) {
- int w2, wn;
-
- /* INTEGER OVERFLOW WILL EVENTUALLY OCCUR HERE
- * IF WORDLENGTH < 48 BITS
- */
- w2 = w + w;
- wn = w2 + 1;
-
- /* APPARENTLY REDUNDANT TESTS MAY BE NECESSARY ON SOME
- * MACHINES, DEPENDING ON HOW INTEGER OVERFLOW IS HANDLED
- */
- if (wn <= w || wn <= w2 || wn <= 0)
- break;
- k = i;
- w = wn;
- }
-
- /* SET MAXIMUM EXPONENT TO (W-1)/4 */
- MP.m = w / 4;
- if (accuracy <= 0) {
- mperr("*** ACCURACY <= 0 IN CALL TO MPSET ***");
- return;
- }
-
- /* B IS THE LARGEST POWER OF 2 SUCH THAT (8*B*B-1) <= W */
- MP.b = pow_ii(2, (k - 3) / 2);
-
- /* Make a multiple of 10 so fractions can be represented exactly */
- b = 1;
- while (MP.b % (10 * b) != MP.b)
- b *= 10;
- MP.b = b;
-
- /* 2E0 BELOW ENSURES AT LEAST ONE GUARD DIGIT */
- MP.t = (int) ((float) (accuracy) * log((float)10.) / log((float) MP.b) +
- (float) 2.0);
-
- /* SEE IF T TOO LARGE FOR DIMENSION STATEMENTS */
- if (MP.t > MP_SIZE) {
- mperr("MP_SIZE TOO SMALL IN CALL TO MPSET, INCREASE MP_SIZE AND DIMENSIONS OF MP ARRAYS TO AT LEAST %d ***", MP.t);
- MP.t = MP_SIZE;
- }
-
- /* CHECK LEGALITY OF B, T, M AND MXR (AT LEAST T+4) */
- mpchk();
-}
-
-
-/* CHECKS LEGALITY OF B, T, M AND MXR.
- * THE CONDITION ON MXR (THE DIMENSION OF R IN COMMON) IS THAT
- * MXR >= (I*T + J)
- */
-// FIXME: MP.t is changed around calls to add/subtract/multiply/divide - it should not be global
-void
-mpchk()
-{
- /* CHECK LEGALITY OF T AND M */
- if (MP.t <= 1)
- mperr("*** T = %d ILLEGAL IN CALL TO MPCHK, PERHAPS NOT SET BEFORE CALL TO AN MP ROUTINE ***", MP.t);
- if (MP.m <= MP.t)
- mperr("*** M <= T IN CALL TO MPCHK, PERHAPS NOT SET BEFORE CALL TO AN MP ROUTINE ***");
+ MPNumber t;
+ mp_set_from_integer(1, &t);
+ mp_epowy(&t, z);
}
-/* SETS Y = ABS(X) FOR MP NUMBERS X AND Y */
void
mp_abs(const MPNumber *x, MPNumber *z)
{
@@ -268,24 +111,24 @@ mp_add3(const MPNumber *x, const MPNumber *y, int *r, int s, int med)
/* CLEAR GUARD DIGITS TO RIGHT OF X DIGITS */
for(i = 3; i >= med; i--)
- r[MP.t + i] = 0;
+ r[MP_T + i] = 0;
if (s >= 0) {
/* HERE DO ADDITION, EXPONENT(Y) >= EXPONENT(X) */
- for (i = MP.t + 3; i >= MP.t; i--)
+ for (i = MP_T + 3; i >= MP_T; i--)
r[i] = x->fraction[i - med];
c = 0;
for (; i >= med; i--) {
c = y->fraction[i] + x->fraction[i - med] + c;
- if (c < MP.b) {
+ if (c < MP_BASE) {
/* NO CARRY GENERATED HERE */
r[i] = c;
c = 0;
} else {
/* CARRY GENERATED HERE */
- r[i] = c - MP.b;
+ r[i] = c - MP_BASE;
c = 1;
}
}
@@ -293,7 +136,7 @@ mp_add3(const MPNumber *x, const MPNumber *y, int *r, int s, int med)
for (; i >= 0; i--)
{
c = y->fraction[i] + c;
- if (c < MP.b) {
+ if (c < MP_BASE) {
r[i] = c;
i--;
@@ -310,7 +153,7 @@ mp_add3(const MPNumber *x, const MPNumber *y, int *r, int s, int med)
/* MUST SHIFT RIGHT HERE AS CARRY OFF END */
if (c != 0) {
- for (i = MP.t + 3; i > 0; i--)
+ for (i = MP_T + 3; i > 0; i--)
r[i] = r[i - 1];
r[0] = 1;
return 1;
@@ -320,7 +163,7 @@ mp_add3(const MPNumber *x, const MPNumber *y, int *r, int s, int med)
}
c = 0;
- for (i = MP.t + med - 1; i >= MP.t; i--) {
+ for (i = MP_T + med - 1; i >= MP_T; i--) {
/* HERE DO SUBTRACTION, ABS(Y) > ABS(X) */
r[i] = c - x->fraction[i - med];
c = 0;
@@ -328,7 +171,7 @@ mp_add3(const MPNumber *x, const MPNumber *y, int *r, int s, int med)
/* BORROW GENERATED HERE */
if (r[i] < 0) {
c = -1;
- r[i] += MP.b;
+ r[i] += MP_BASE;
}
}
@@ -340,7 +183,7 @@ mp_add3(const MPNumber *x, const MPNumber *y, int *r, int s, int med)
c = 0;
} else {
/* BORROW GENERATED HERE */
- r[i] = c + MP.b;
+ r[i] = c + MP_BASE;
c = -1;
}
}
@@ -359,7 +202,7 @@ mp_add3(const MPNumber *x, const MPNumber *y, int *r, int s, int med)
return 0;
}
- r[i] = c + MP.b;
+ r[i] = c + MP_BASE;
c = -1;
}
@@ -375,6 +218,8 @@ mp_add2(const MPNumber *x, const MPNumber *y, int y_sign, MPNumber *z)
int exp_diff, med;
int x_largest = 0;
MPNumber zt; // Use stack variable because of mp_normalize brokeness
+
+ memset(&zt, 0, sizeof(MPNumber));
/* X = 0 OR NEGLIGIBLE, SO RESULT = +-Y */
if (x->sign == 0) {
@@ -392,9 +237,8 @@ mp_add2(const MPNumber *x, const MPNumber *y, int y_sign, MPNumber *z)
/* COMPARE SIGNS */
sign_prod = y_sign * x->sign;
if (abs(sign_prod) > 1) {
- mpchk();
mperr("*** SIGN NOT 0, +1 OR -1 IN MP_ADD2 CALL, POSSIBLE OVERWRITING PROBLEM ***");
- z->sign = 0;
+ mp_set_from_integer(0, z);
return;
}
@@ -403,7 +247,7 @@ mp_add2(const MPNumber *x, const MPNumber *y, int y_sign, MPNumber *z)
med = abs(exp_diff);
if (exp_diff < 0) {
/* HERE EXPONENT(Y) > EXPONENT(X) */
- if (med > MP.t) {
+ if (med > MP_T) {
/* 'y' so much larger than 'x' that 'x+-y = y'. Warning: still true with rounding?? */
mp_set_from_mp(y, z);
z->sign = y_sign;
@@ -412,7 +256,7 @@ mp_add2(const MPNumber *x, const MPNumber *y, int y_sign, MPNumber *z)
x_largest = 0;
} else if (exp_diff > 0) {
/* HERE EXPONENT(X) > EXPONENT(Y) */
- if (med > MP.t) {
+ if (med > MP_T) {
/* 'x' so much larger than 'y' that 'x+-y = x'. Warning: still true with rounding?? */
mp_set_from_mp(x, z);
return;
@@ -423,7 +267,7 @@ mp_add2(const MPNumber *x, const MPNumber *y, int y_sign, MPNumber *z)
if (sign_prod < 0) {
/* Signs are not equal. find out which mantissa is larger. */
int j;
- for (j = 0; j < MP.t; j++) {
+ for (j = 0; j < MP_T; j++) {
int i = x->fraction[j] - y->fraction[j];
if (i == 0)
continue;
@@ -435,8 +279,8 @@ mp_add2(const MPNumber *x, const MPNumber *y, int y_sign, MPNumber *z)
}
/* Both mantissas equal, so result is zero. */
- if (j >= MP.t) {
- z->sign = 0;
+ if (j >= MP_T) {
+ mp_set_from_integer(0, z);
return;
}
}
@@ -451,7 +295,7 @@ mp_add2(const MPNumber *x, const MPNumber *y, int y_sign, MPNumber *z)
zt.exponent = y->exponent + mp_add3(x, y, zt.fraction, sign_prod, med);
}
mp_normalize(&zt, 0);
- mp_set_from_mp(&zt, z);
+ mp_set_from_mp(&zt, z);
}
@@ -465,119 +309,81 @@ 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)
{
MPNumber t;
- mpchk();
mp_set_from_integer(y, &t);
mp_add(x, &t, 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)
{
MPNumber t;
- mpchk();
mp_set_from_fraction(i, j, &t);
mp_add(x, &t, 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 *y)
+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) {
- y->sign = 0;
+ 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, y);
+ mp_set_from_mp(x, z);
return;
}
-
- /* MOVE FRACTIONAL PART OF X TO ACCUMULATOR */
- y->sign = x->sign;
- y->exponent = x->exponent;
- memset(y->fraction, 0, y->exponent*sizeof(int));
- if (x != y)
- memcpy(y->fraction + y->exponent, x->fraction + x->exponent,
- (MP.t - x->exponent)*sizeof(int));
- memset(y->fraction + MP.t, 0, 4*sizeof(int));
-
- /* NORMALIZE RESULT AND RETURN */
- mp_normalize(y, 1);
+
+ /* 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 - 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 *y)
+mp_integer_component(const MPNumber *x, MPNumber *z)
{
int i;
-
- mpchk();
- mp_set_from_mp(x, y);
- if (y->sign == 0)
- return;
-
- /* IF EXPONENT LARGE ENOUGH RETURN Y = X */
- if (y->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 (y->exponent <= 0) {
- y->sign = 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 = y->exponent; i < MP.t; i++) {
- y->fraction[i] = 0;
- }
+ /* Clear fraction */
+ mp_set_from_mp(x, z);
+ for (i = z->exponent; i < MP_SIZE; i++)
+ z->fraction[i] = 0;
}
-/* COMPARES MP NUMBER X WITH REAL NUMBER RI, RETURNING
- * +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
- */
-static int
-mp_compare_mp_to_float(const MPNumber *x, float ri)
-{
- MPNumber t;
- mpchk();
-
- /* CONVERT RI TO MULTIPLE-PRECISION AND COMPARE */
- mp_set_from_float(ri, &t);
- return mp_compare_mp_to_mp(x, &t);
-}
/* COMPARES THE MULTIPLE-PRECISION NUMBERS X AND Y,
* RETURNING +1 IF X > Y,
@@ -608,7 +414,7 @@ mp_compare_mp_to_mp(const MPNumber *x, const MPNumber *y)
/* ABS(X) > ABS(Y) */
return x->sign;
}
- for (i = 0; i < MP.t; i++) {
+ for (i = 0; i < MP_T; i++) {
i2 = x->fraction[i] - y->fraction[i];
if (i2 < 0) {
/* ABS(X) < ABS(Y) */
@@ -635,24 +441,19 @@ mp_divide(const MPNumber *x, const MPNumber *y, MPNumber *z)
int i, ie, iz3;
MPNumber t;
- mpchk();
-
- /* CHECK FOR DIVISION BY ZERO */
+ /* x/0 */
if (y->sign == 0) {
mperr("*** ATTEMPTED DIVISION BY ZERO IN CALL TO MP_DIVIDE ***");
- z->sign = 0;
+ mp_set_from_integer(0, z);
return;
}
- /* CHECK FOR X = 0 */
+ /* 0/y = 0 */
if (x->sign == 0) {
- z->sign = 0;
+ mp_set_from_integer(0, z);
return;
}
- /* INCREASE M TO AVOID OVERFLOW IN MP_RECIPROCAL */
- MP.m += 2;
-
/* FORM RECIPROCAL OF Y */
mp_reciprocal(y, &t);
@@ -666,15 +467,7 @@ mp_divide(const MPNumber *x, const MPNumber *y, MPNumber *z)
iz3 = z->fraction[0];
mpext(i, iz3, z);
- /* RESTORE M, CORRECT EXPONENT AND RETURN */
- MP.m += -2;
z->exponent += ie;
-
- /* Check for overflow or underflow */
- if (z->exponent < -MP.m)
- mpunfl(z);
- else if (z->exponent > MP.m)
- mpovfl(z, "*** OVERFLOW OCCURRED IN MP_DIVIDE ***");
}
@@ -690,7 +483,7 @@ mp_divide_integer(const MPNumber *x, int iy, MPNumber *z)
if (iy == 0) {
mperr("*** ATTEMPTED DIVISION BY ZERO IN CALL TO MP_DIVIDE_INTEGER ***");
- z->sign = 0;
+ mp_set_from_integer(0, z);
return;
}
@@ -704,12 +497,7 @@ mp_divide_integer(const MPNumber *x, int iy, MPNumber *z)
z->exponent = x->exponent;
/* CHECK FOR DIVISION BY B */
- if (iy == MP.b) {
- if (x->exponent <= -MP.m)
- {
- mpunfl(z);
- return;
- }
+ if (iy == MP_BASE) {
mp_set_from_mp(x, z);
z->exponent--;
return;
@@ -724,7 +512,7 @@ mp_divide_integer(const MPNumber *x, int iy, MPNumber *z)
}
c = 0;
- i2 = MP.t + 4;
+ i2 = MP_T + 4;
i = 0;
/* IF IY*B NOT REPRESENTABLE AS AN INTEGER HAVE TO SIMULATE
@@ -732,12 +520,12 @@ mp_divide_integer(const MPNumber *x, int iy, MPNumber *z)
*/
/* Computing MAX */
- b2 = max(MP.b << 3, 32767 / MP.b);
+ b2 = max(MP_BASE << 3, 32767 / MP_BASE);
if (iy < b2) {
/* LOOK FOR FIRST NONZERO DIGIT IN QUOTIENT */
do {
- c = MP.b * c;
- if (i < MP.t)
+ c = MP_BASE * c;
+ if (i < MP_T)
c += x->fraction[i];
i++;
r1 = c / iy;
@@ -748,14 +536,14 @@ mp_divide_integer(const MPNumber *x, int iy, MPNumber *z)
/* ADJUST EXPONENT AND GET T+4 DIGITS IN QUOTIENT */
z->exponent += 1 - i;
z->fraction[0] = r1;
- c = MP.b * (c - iy * r1);
+ c = MP_BASE * (c - iy * r1);
kh = 1;
- if (i < MP.t) {
- kh = MP.t + 1 - i;
+ if (i < MP_T) {
+ kh = MP_T + 1 - i;
for (k = 1; k < kh; k++) {
c += x->fraction[i];
z->fraction[k] = c / iy;
- c = MP.b * (c - iy * z->fraction[k]);
+ c = MP_BASE * (c - iy * z->fraction[k]);
i++;
}
if (c < 0)
@@ -764,7 +552,7 @@ mp_divide_integer(const MPNumber *x, int iy, MPNumber *z)
for (k = kh; k < i2; k++) {
z->fraction[k] = c / iy;
- c = MP.b * (c - iy * z->fraction[k]);
+ c = MP_BASE * (c - iy * z->fraction[k]);
}
if (c < 0)
goto L210;
@@ -776,15 +564,15 @@ mp_divide_integer(const MPNumber *x, int iy, MPNumber *z)
}
/* HERE NEED SIMULATED DOUBLE-PRECISION DIVISION */
- j1 = iy / MP.b;
+ j1 = iy / MP_BASE;
/* LOOK FOR FIRST NONZERO DIGIT */
c2 = 0;
- j2 = iy - j1 * MP.b;
+ j2 = iy - j1 * MP_BASE;
do {
- c = MP.b * c + c2;
+ c = MP_BASE * c + c2;
i__1 = c - j1;
- c2 = i < MP.t ? x->fraction[i] : 0;
+ c2 = i < MP_T ? x->fraction[i] : 0;
i++;
} while (i__1 < 0 || (i__1 == 0 && c2 < j2));
@@ -807,14 +595,14 @@ mp_divide_integer(const MPNumber *x, int iy, MPNumber *z)
iq -= j1;
}
- iq = iq * MP.b - ir * j2;
+ iq = iq * MP_BASE - ir * j2;
if (iq < 0) {
/* HERE IQ NEGATIVE SO IR WAS TOO LARGE */
ir--;
iq += iy;
}
- if (i < MP.t)
+ if (i < MP_T)
iq += x->fraction[i];
i++;
iqj = iq / iy;
@@ -835,9 +623,8 @@ mp_divide_integer(const MPNumber *x, int iy, MPNumber *z)
L210:
/* CARRY NEGATIVE SO OVERFLOW MUST HAVE OCCURRED */
- mpchk();
mperr("*** INTEGER OVERFLOW IN MP_DIVIDE_INTEGER, B TOO LARGE ***");
- z->sign = 0;
+ mp_set_from_integer(0, z);
}
@@ -846,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.
@@ -854,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;*/
}
@@ -866,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)
{
@@ -874,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
@@ -901,41 +766,26 @@ void
mp_epowy(const MPNumber *x, MPNumber *z)
{
float r__1;
- int i, ie, ix, ts, xs, tss;
+ int i, ix, xs, tss;
float rx, rz, rlb;
MPNumber t1, t2;
- mpchk();
-
- /* 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;
}
/* SEE IF ABS(X) SO LARGE THAT EXP(X) WILL CERTAINLY OVERFLOW
* OR UNDERFLOW. 1.01 IS TO ALLOW FOR ERRORS IN ALOG.
*/
- rlb = log((float) MP.b) * (float)1.01;
- if (mp_compare_mp_to_float(x, -(double)((float) (MP.m + 1)) * rlb) < 0) {
- /* UNDERFLOW SO CALL MPUNFL AND RETURN */
- mpunfl(z);
- return;
- }
-
- if (mp_compare_mp_to_float(x, (float) MP.m * rlb) > 0) {
- /* OVERFLOW HERE */
- mpovfl(z, "*** OVERFLOW IN SUBROUTINE MP_EPOWY ***");
- return;
- }
+ rlb = log((float)MP_BASE) * 1.01f;
/* NOW SAFE TO CONVERT X TO REAL */
rx = mp_cast_to_float(x);
@@ -947,9 +797,9 @@ mp_epowy(const MPNumber *x, MPNumber *z)
/* IF ABS(X) > M POSSIBLE THAT INT(X) OVERFLOWS,
* SO DIVIDE BY 32.
*/
- if (fabs(rx) > (float) MP.m) {
+ /*if (fabs(rx) > (float)MP.m) {
mp_divide_integer(&t2, 32, &t2);
- }
+ }*/
/* GET FRACTIONAL AND INTEGER PARTS OF ABS(X) */
ix = mp_cast_to_int(&t2);
@@ -957,85 +807,58 @@ 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)
*/
- if (MP.t < 4)
- tss = MP.t + 1;
+ if (MP_T < 4)
+ tss = MP_T + 1;
else
- tss = MP.t + 2;
+ tss = MP_T + 2;
/* LOOP FOR E COMPUTATION. DECREASE T IF POSSIBLE. */
/* Computing MIN */
- ts = MP.t;
- MP.t = tss;
mp_set_from_integer(xs, &t1);
- MP.t = ts;
t2.sign = 0;
for (i = 2 ; ; i++) {
- int t;
-
- t = min(tss, tss + 2 + t1.exponent);
- if (t <= 2)
+ if (min(tss, tss + 2 + t1.exponent) <= 2)
break;
- ts = MP.t;
- MP.t = t;
mp_divide_integer(&t1, i * xs, &t1);
- MP.t = ts;
-
- ts = MP.t;
- MP.t = tss;
mp_add(&t2, &t1, &t2);
- MP.t = ts;
if (t1.sign == 0)
break;
}
/* RAISE E OR 1/E TO POWER IX */
- ts = MP.t;
- MP.t = tss;
if (xs > 0)
mp_add_integer(&t2, 2, &t2);
- mp_pwr_integer(&t2, ix, &t2);
- MP.t = ts;
+ mp_xpowy_integer(&t2, ix, &t2);
/* MULTIPLY EXPS OF INTEGER AND FRACTIONAL PARTS */
mp_multiply(z, &t2, z);
/* MUST CORRECT RESULT IF DIVIDED BY 32 ABOVE. */
- if (fabs(rx) > (float) MP.m && z->sign != 0) {
+ /*if (fabs(rx) > (float)MP.m && z->sign != 0) {
for (i = 1; i <= 5; ++i) {
- /* SAVE EXPONENT TO AVOID OVERFLOW IN MP_MULTIPLY */
+ // SAVE EXPONENT TO AVOID OVERFLOW IN MP_MULTIPLY
ie = z->exponent;
z->exponent = 0;
mp_multiply(z, z, z);
z->exponent += ie << 1;
-
- /* CHECK FOR UNDERFLOW AND OVERFLOW */
- if (z->exponent < -MP.m) {
- mpunfl(z);
- return;
- }
- if (z->exponent > MP.m) {
- mpovfl(z, "*** OVERFLOW IN SUBROUTINE MP_EPOWY ***");
- return;
- }
}
- }
+ }*/
/* CHECK THAT RELATIVE ERROR LESS THAN 0.01 UNLESS ABS(X) LARGE
* (WHEN EXP MIGHT OVERFLOW OR UNDERFLOW)
*/
- if (fabs(rx) > (float)10.0)
+ if (fabs(rx) > 10.0f)
return;
rz = mp_cast_to_float(z);
- if ((r__1 = rz - exp(rx), fabs(r__1)) < rz * (float)0.01)
+ if ((r__1 = rz - exp(rx), fabs(r__1)) < rz * 0.01f)
return;
/* THE FOLLOWING MESSAGE MAY INDICATE THAT
@@ -1046,96 +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 *y)
-{
- int i, q, ib, ic;
- float rlb;
- MPNumber t1, t2;
-
- mpchk();
-
- /* CHECK FOR X = 0 */
- if (x->sign == 0) {
- y->sign = 0;
- return;
- }
-
- /* CHECK THAT ABS(X) < 1 */
- if (x->exponent > 0) {
- mperr("*** ABS(X) NOT LESS THAN 1 IN CALL TO MPEXP1 ***");
- y->sign = 0;
- return;
- }
-
- mp_set_from_mp(x, &t1);
- rlb = log((float) MP.b);
-
- /* COMPUTE APPROXIMATELY OPTIMAL Q (AND DIVIDE X BY 2**Q) */
- q = (int) (sqrt((float) MP.t * (float).48 * rlb) + (float) x->exponent *
- (float)1.44 * rlb);
-
- /* HALVE Q TIMES */
- if (q > 0) {
- ib = MP.b << 2;
- ic = 1;
- for (i = 1; i <= q; ++i) {
- ic <<= 1;
- if (ic < ib && ic != MP.b && i < q)
- continue;
- mp_divide_integer(&t1, ic, &t1);
- ic = 1;
- }
- }
-
- if (t1.sign == 0) {
- y->sign = 0;
- return;
- }
- mp_set_from_mp(&t1, y);
- mp_set_from_mp(&t1, &t2);
-
- /* SUM SERIES, REDUCING T WHERE POSSIBLE */
- for (i = 2; ; i++) {
- int t, ts;
-
- t = MP.t + 2 + t2.exponent - y->exponent;
- if (t <= 2)
- break;
- t = min(t, MP.t);
-
- ts = MP.t;
- MP.t = t;
- mp_multiply(&t1, &t2, &t2);
- mp_divide_integer(&t2, i, &t2);
- MP.t = ts;
-
- mp_add(&t2, y, y);
- 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(y, 2, &t1);
- mp_multiply(&t1, y, y);
- }
-}
-
-
/* 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
@@ -1217,28 +950,25 @@ 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
-mplns(const MPNumber *x, MPNumber *y)
+mplns(const MPNumber *x, MPNumber *z)
{
int t, it0;
MPNumber t1, t2, t3;
- mpchk();
-
- /* CHECK FOR X = 0 EXACTLY */
+ /* ln(1) = 0 */
if (x->sign == 0) {
- y->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 ***");
- y->sign = 0;
+ mp_set_from_integer(0, z);
return;
}
@@ -1250,27 +980,25 @@ mplns(const MPNumber *x, MPNumber *y)
mp_add_fraction(&t1, 1, 2, &t1);
mp_multiply(x, &t1, &t1);
mp_add_integer(&t1, -1, &t1);
- mp_multiply(x, &t1, y);
+ mp_multiply(x, &t1, z);
/* START NEWTON ITERATION USING SMALL T, LATER INCREASE */
- t = max(5, 13 - (MP.b << 1));
- if (t <= MP.t)
+ t = max(5, 13 - (MP_BASE << 1));
+ if (t <= MP_T)
{
it0 = (t + 5) / 2;
while(1)
{
- int ts, ts2, ts3;
+ int ts2, ts3;
- ts = MP.t;
- MP.t = t;
- mpexp1(y, &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);
- mp_subtract(y, &t3, y);
- MP.t = ts;
- if (t >= MP.t)
+ mp_subtract(z, &t3, z);
+ if (t >= MP_T)
break;
/* FOLLOWING LOOP COMPUTES NEXT VALUE OF T TO USE.
@@ -1278,7 +1006,7 @@ mplns(const MPNumber *x, MPNumber *y)
* WE CAN ALMOST DOUBLE T EACH TIME.
*/
ts3 = t;
- t = MP.t;
+ t = MP_T;
do {
ts2 = t;
t = (t + it0) / 2;
@@ -1287,13 +1015,13 @@ mplns(const MPNumber *x, MPNumber *y)
}
/* CHECK THAT NEWTON ITERATION WAS CONVERGING AS EXPECTED */
- if (t3.sign != 0 && t3.exponent << 1 > it0 - MP.t) {
+ if (t3.sign != 0 && t3.exponent << 1 > it0 - MP_T) {
mperr("*** ERROR OCCURRED IN MPLNS, NEWTON ITERATION NOT CONVERGING PROPERLY ***");
}
}
/* REVERSE SIGN OF Y AND RETURN */
- y->sign = -y->sign;
+ z->sign = -z->sign;
}
@@ -1303,7 +1031,7 @@ mplns(const MPNumber *x, MPNumber *y)
* 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
*/
@@ -1314,18 +1042,16 @@ mp_ln(const MPNumber *x, MPNumber *z)
float rx, rlx;
MPNumber t1, t2;
- mpchk();
-
/* CHECK THAT X IS POSITIVE */
if (x->sign <= 0) {
mperr("*** X NONPOSITIVE IN CALL TO MP_LN ***");
- z->sign = 0;
+ mp_set_from_integer(0, z);
return;
}
/* MOVE X TO LOCAL STORAGE */
mp_set_from_mp(x, &t1);
- z->sign = 0;
+ mp_set_from_integer(0, z);
k = 0;
/* LOOP TO GET APPROXIMATE LN(X) USING SINGLE-PRECISION */
@@ -1348,7 +1074,7 @@ mp_ln(const MPNumber *x, MPNumber *z)
/* RESTORE EXPONENT AND COMPUTE SINGLE-PRECISION LOG */
t1.exponent = e;
- rlx = log(rx) + (float) e * log((float) MP.b);
+ rlx = log(rx) + (float)e * log((float)MP_BASE);
mp_set_from_float(-(double)rlx, &t2);
/* UPDATE Z AND COMPUTE ACCURATE EXP OF APPROXIMATE LOG */
@@ -1373,7 +1099,7 @@ mp_ln(const MPNumber *x, MPNumber *z)
* 1. log10(x) = log10(e) * log(x)
*/
void
-mp_logarithm(int n, MPNumber *x, MPNumber *z)
+mp_logarithm(int n, const MPNumber *x, MPNumber *z)
{
MPNumber t1, t2;
@@ -1411,9 +1137,10 @@ mp_multiply(const MPNumber *x, const MPNumber *y, MPNumber *z)
int c, i, j, i2, xi;
MPNumber r;
- mpchk();
- i2 = MP.t + 4;
-
+ i2 = MP_T + 4;
+
+ memset(&r, 0, sizeof(MPNumber));
+
/* FORM SIGN OF PRODUCT */
z->sign = x->sign * y->sign;
if (z->sign == 0)
@@ -1422,13 +1149,9 @@ mp_multiply(const MPNumber *x, const MPNumber *y, MPNumber *z)
/* FORM EXPONENT OF PRODUCT */
z->exponent = x->exponent + y->exponent;
- /* CLEAR ACCUMULATOR */
- for (i = 0; i < i2; i++)
- r.fraction[i] = 0;
-
/* PERFORM MULTIPLICATION */
c = 8;
- for (i = 0; i < MP.t; i++) {
+ for (i = 0; i < MP_T; i++) {
xi = x->fraction[i];
/* FOR SPEED, PUT THE NUMBER WITH MANY ZEROS FIRST */
@@ -1436,16 +1159,16 @@ mp_multiply(const MPNumber *x, const MPNumber *y, MPNumber *z)
continue;
/* Computing MIN */
- for (j = 0; j < min(MP.t, i2 - i - 1); j++)
+ for (j = 0; j < min(MP_T, i2 - i - 1); j++)
r.fraction[i+j+1] += xi * y->fraction[j];
c--;
if (c > 0)
continue;
/* CHECK FOR LEGAL BASE B DIGIT */
- if (xi < 0 || xi >= MP.b) {
+ if (xi < 0 || xi >= MP_BASE) {
mperr("*** ILLEGAL BASE B DIGIT IN CALL TO MP_MULTIPLY, POSSIBLE OVERWRITING PROBLEM ***");
- z->sign = 0;
+ mp_set_from_integer(0, z);
return;
}
@@ -1456,24 +1179,24 @@ mp_multiply(const MPNumber *x, const MPNumber *y, MPNumber *z)
int ri = r.fraction[j] + c;
if (ri < 0) {
mperr("*** INTEGER OVERFLOW IN MP_MULTIPLY, B TOO LARGE ***");
- z->sign = 0;
+ mp_set_from_integer(0, z);
return;
}
- c = ri / MP.b;
- r.fraction[j] = ri - MP.b * c;
+ c = ri / MP_BASE;
+ r.fraction[j] = ri - MP_BASE * c;
}
if (c != 0) {
mperr("*** ILLEGAL BASE B DIGIT IN CALL TO MP_MULTIPLY, POSSIBLE OVERWRITING PROBLEM ***");
- z->sign = 0;
+ mp_set_from_integer(0, z);
return;
}
c = 8;
}
if (c != 8) {
- if (xi < 0 || xi >= MP.b) {
+ if (xi < 0 || xi >= MP_BASE) {
mperr("*** ILLEGAL BASE B DIGIT IN CALL TO MP_MULTIPLY, POSSIBLE OVERWRITING PROBLEM ***");
- z->sign = 0;
+ mp_set_from_integer(0, z);
return;
}
@@ -1482,23 +1205,23 @@ mp_multiply(const MPNumber *x, const MPNumber *y, MPNumber *z)
int ri = r.fraction[j] + c;
if (ri < 0) {
mperr("*** INTEGER OVERFLOW IN MP_MULTIPLY, B TOO LARGE ***");
- z->sign = 0;
+ mp_set_from_integer(0, z);
return;
}
- c = ri / MP.b;
- r.fraction[j] = ri - MP.b * c;
+ c = ri / MP_BASE;
+ r.fraction[j] = ri - MP_BASE * c;
}
if (c != 0) {
mperr("*** ILLEGAL BASE B DIGIT IN CALL TO MP_MULTIPLY, POSSIBLE OVERWRITING PROBLEM ***");
- z->sign = 0;
+ mp_set_from_integer(0, z);
return;
}
}
/* NORMALIZE AND ROUND RESULT */
// FIXME: Use stack variable because of mp_normalize brokeness
- for (i = 0; i < i2; i++)
+ for (i = 0; i < MP_SIZE; i++)
z->fraction[i] = r.fraction[i];
mp_normalize(z, 0);
}
@@ -1517,7 +1240,7 @@ mpmul2(const MPNumber *x, int iy, MPNumber *z, int trunc)
z->sign = x->sign;
if (z->sign == 0 || iy == 0) {
- z->sign = 0;
+ mp_set_from_integer(0, z);
return;
}
@@ -1526,16 +1249,10 @@ mpmul2(const MPNumber *x, int iy, MPNumber *z, int trunc)
z->sign = -z->sign;
/* CHECK FOR MULTIPLICATION BY B */
- if (iy == MP.b) {
- if (x->exponent < MP.m) {
- mp_set_from_mp(x, z);
- z->sign = z->sign;
- z->exponent = x->exponent + 1;
- }
- else {
- mpchk();
- mpovfl(z, "*** OVERFLOW OCCURRED IN MPMUL2 ***");
- }
+ if (iy == MP_BASE) {
+ mp_set_from_mp(x, z);
+ z->sign = z->sign;
+ z->exponent = x->exponent + 1;
return;
}
}
@@ -1545,48 +1262,47 @@ mpmul2(const MPNumber *x, int iy, MPNumber *z, int trunc)
/* FORM PRODUCT IN ACCUMULATOR */
c = 0;
- t1 = MP.t + 1;
- t3 = MP.t + 3;
- t4 = MP.t + 4;
+ t1 = MP_T + 1;
+ t3 = MP_T + 3;
+ t4 = MP_T + 4;
/* IF IY*B NOT REPRESENTABLE AS AN INTEGER WE HAVE TO SIMULATE
* DOUBLE-PRECISION MULTIPLICATION.
*/
/* Computing MAX */
- if (iy >= max(MP.b << 3, 32767 / MP.b)) {
+ if (iy >= max(MP_BASE << 3, 32767 / MP_BASE)) {
/* HERE J IS TOO LARGE FOR SINGLE-PRECISION MULTIPLICATION */
- j1 = iy / MP.b;
- j2 = iy - j1 * MP.b;
+ j1 = iy / MP_BASE;
+ j2 = iy - j1 * MP_BASE;
/* FORM PRODUCT */
for (ij = 1; ij <= t4; ++ij) {
- c1 = c / MP.b;
- c2 = c - MP.b * c1;
+ c1 = c / MP_BASE;
+ c2 = c - MP_BASE * c1;
i = t1 - ij;
ix = 0;
if (i > 0)
ix = x->fraction[i - 1];
ri = j2 * ix + c2;
- is = ri / MP.b;
+ is = ri / MP_BASE;
c = j1 * ix + c1 + is;
- z->fraction[i + 3] = ri - MP.b * is;
+ z->fraction[i + 3] = ri - MP_BASE * is;
}
}
else
{
- for (ij = 1; ij <= MP.t; ++ij) {
+ for (ij = 1; ij <= MP_T; ++ij) {
i = t1 - ij;
ri = iy * x->fraction[i - 1] + c;
- c = ri / MP.b;
- z->fraction[i + 3] = ri - MP.b * c;
+ c = ri / MP_BASE;
+ z->fraction[i + 3] = ri - MP_BASE * c;
}
/* CHECK FOR INTEGER OVERFLOW */
if (ri < 0) {
- mpchk();
mperr("*** INTEGER OVERFLOW IN MPMUL2, B TOO LARGE ***");
- z->sign = 0;
+ mp_set_from_integer(0, z);
return;
}
@@ -1594,8 +1310,8 @@ mpmul2(const MPNumber *x, int iy, MPNumber *z, int trunc)
for (ij = 1; ij <= 4; ++ij) {
i = 5 - ij;
ri = c;
- c = ri / MP.b;
- z->fraction[i - 1] = ri - MP.b * c;
+ c = ri / MP_BASE;
+ z->fraction[i - 1] = ri - MP_BASE * c;
}
}
@@ -1609,9 +1325,8 @@ mpmul2(const MPNumber *x, int iy, MPNumber *z, int trunc)
}
if (c < 0) {
- mpchk();
mperr("*** INTEGER OVERFLOW IN MPMUL2, B TOO LARGE ***");
- z->sign = 0;
+ mp_set_from_integer(0, z);
return;
}
@@ -1620,8 +1335,8 @@ mpmul2(const MPNumber *x, int iy, MPNumber *z, int trunc)
z->fraction[i] = z->fraction[i - 1];
}
ri = c;
- c = ri / MP.b;
- z->fraction[0] = ri - MP.b * c;
+ c = ri / MP_BASE;
+ z->fraction[0] = ri - MP_BASE * c;
z->exponent++;
}
}
@@ -1646,14 +1361,13 @@ mp_multiply_fraction(const MPNumber *x, int numerator, int denominator, MPNumber
int is, js;
if (denominator == 0) {
- mpchk();
mperr("*** ATTEMPTED DIVISION BY ZERO IN MP_MULTIPLY_FRACTION ***");
- z->sign = 0;
+ mp_set_from_integer(0, z);
return;
}
if (numerator == 0) {
- z->sign = 0;
+ mp_set_from_integer(0, z);
return;
}
@@ -1683,7 +1397,7 @@ mp_invert_sign(const MPNumber *x, MPNumber *y)
* AND RETURNS MP RESULT IN Z. INTEGER ARGUMENTS REG_EXP IS
* NOT PRESERVED. R*-ROUNDING IS USED IF TRUNC == 0
*/
-// FIXME: Is r->fraction large enough? It seems to be in practise but it may be MP.t+4 instead of MP.t
+// FIXME: Is r->fraction large enough? It seems to be in practise but it may be MP_T+4 instead of MP_T
// FIXME: There is some sort of stack corruption/use of unitialised variables here. Some functions are
// using stack variables as x otherwise there are corruption errors. e.g. "Cos(45) - 1/Sqrt(2) = -0"
// (try in scientific mode)
@@ -1699,18 +1413,18 @@ mp_normalize(MPNumber *x, int trunc)
/* CHECK THAT SIGN = +-1 */
if (abs(x->sign) > 1) {
mperr("*** SIGN NOT 0, +1 OR -1 IN CALL TO MP_NORMALIZE, POSSIBLE OVERWRITING PROBLEM ***");
- x->sign = 0;
+ mp_set_from_integer(0, x);
return;
}
- i2 = MP.t + 4;
+ i2 = MP_T + 4;
/* Normalize by shifting the fraction to the left */
if (x->fraction[0] == 0) {
/* Find how many places to shift and detect 0 fraction */
for (i = 1; i < i2 && x->fraction[i] == 0; i++);
if (i == i2) {
- x->sign = 0;
+ mp_set_from_integer(0, x);
return;
}
@@ -1727,12 +1441,12 @@ mp_normalize(MPNumber *x, int trunc)
/* SEE IF ROUNDING NECESSARY
* TREAT EVEN AND ODD BASES DIFFERENTLY
*/
- b2 = MP.b / 2;
- if (b2 << 1 != MP.b) {
+ b2 = MP_BASE / 2;
+ if (b2 << 1 != MP_BASE) {
round = 0;
/* ODD BASE, ROUND IF R(T+1)... > 1/2 */
for (i = 0; i < 4; i++) {
- i__1 = x->fraction[MP.t + i] - b2;
+ i__1 = x->fraction[MP_T + i] - b2;
if (i__1 < 0)
break;
else if (i__1 == 0)
@@ -1748,12 +1462,12 @@ mp_normalize(MPNumber *x, int trunc)
* AFTER R(T+2).
*/
round = 1;
- i__1 = x->fraction[MP.t] - b2;
+ i__1 = x->fraction[MP_T] - b2;
if (i__1 < 0)
round = 0;
else if (i__1 == 0) {
- if (x->fraction[MP.t - 1] % 2 != 0) {
- if (x->fraction[MP.t + 1] + x->fraction[MP.t + 2] + x->fraction[MP.t + 3] == 0) {
+ if (x->fraction[MP_T - 1] % 2 != 0) {
+ if (x->fraction[MP_T + 1] + x->fraction[MP_T + 2] + x->fraction[MP_T + 3] == 0) {
round = 0;
}
}
@@ -1762,9 +1476,9 @@ mp_normalize(MPNumber *x, int trunc)
/* ROUND */
if (round) {
- for (j = MP.t - 1; j >= 0; j--) {
+ for (j = MP_T - 1; j >= 0; j--) {
++x->fraction[j];
- if (x->fraction[j] < MP.b)
+ if (x->fraction[j] < MP_BASE)
break;
x->fraction[j] = 0;
}
@@ -1776,112 +1490,37 @@ mp_normalize(MPNumber *x, int trunc)
}
}
}
-
- /* Check for over and underflow */
- if (x->exponent > MP.m) {
- mpovfl(x, "*** OVERFLOW OCCURRED IN MP_NORMALIZE ***");
- return;
- }
- if (x->exponent < -MP.m) {
- mpunfl(x);
- return;
- }
}
-/* 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 *y)
-{
- int n2, ns;
- MPNumber t;
-
- n2 = n;
- if (n2 < 0) {
- /* N < 0 */
- mpchk();
- n2 = -n2;
- if (x->sign == 0) {
- mperr("*** ATTEMPT TO RAISE ZERO TO NEGATIVE POWER IN CALL TO SUBROUTINE MP_PWR_INTEGER ***");
- y->sign = 0;
- return;
- }
- } else if (n2 == 0) {
- /* N == 0, RETURN Y = 1. */
- mp_set_from_integer(1, y);
- return;
- } else {
- /* N > 0 */
- mpchk();
- if (x->sign == 0) {
- y->sign = 0;
- return;
- }
- }
-
- /* MOVE X */
- mp_set_from_mp(x, y);
-
- /* IF N < 0 FORM RECIPROCAL */
- if (n < 0)
- mp_reciprocal(y, y);
- mp_set_from_mp(y, &t);
-
- /* SET PRODUCT TERM TO ONE */
- mp_set_from_integer(1, y);
-
- /* MAIN LOOP, LOOK AT BITS OF N2 FROM RIGHT */
- while(1) {
- ns = n2;
- n2 /= 2;
- if (n2 << 1 != ns)
- mp_multiply(y, &t, y);
- 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;
- mpchk();
-
+ /* (-x)^y imaginary */
if (x->sign < 0) {
mperr(_("Negative X and non-integer Y not supported"));
- z->sign = 0;
- }
- 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 ***");
- }
- z->sign = 0;
+ 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^-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;
+ }
+
+ /* 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);
}
@@ -1903,21 +1542,15 @@ mp_reciprocal(const MPNumber *x, MPNumber *z)
int ex, it0, t;
float rx;
- /* CHECK LEGALITY OF B, T, M AND MXR */
- mpchk();
-
/* MP_ADD_INTEGER REQUIRES 2T+6 WORDS. */
if (x->sign == 0) {
mperr("*** ATTEMPTED DIVISION BY ZERO IN CALL TO MP_RECIPROCAL ***");
- z->sign = 0;
+ mp_set_from_integer(0, z);
return;
}
ex = x->exponent;
- /* TEMPORARILY INCREASE M TO AVOID OVERFLOW */
- MP.m += 2;
-
/* SET EXPONENT TO ZERO SO RX NOT TOO LARGE OR SMALL. */
/* work-around to avoid touching X */
mp_set_from_mp(x, &tmp_x);
@@ -1925,56 +1558,44 @@ mp_reciprocal(const MPNumber *x, MPNumber *z)
rx = mp_cast_to_float(&tmp_x);
/* USE SINGLE-PRECISION RECIPROCAL AS FIRST APPROXIMATION */
- mp_set_from_float((float)1. / rx, &t1);
+ mp_set_from_float(1.0f / rx, &t1);
/* CORRECT EXPONENT OF FIRST APPROXIMATION */
t1.exponent -= ex;
/* START WITH SMALL T TO SAVE TIME. ENSURE THAT B**(T-1) >= 100 */
- if (MP.b < 10)
- t = it[MP.b - 1];
+ if (MP_BASE < 10)
+ t = it[MP_BASE - 1];
else
t = 3;
it0 = (t + 4) / 2;
/* MAIN ITERATION LOOP */
- if (t <= MP.t) {
+ if (t <= MP_T) {
while(1) {
- int ts, ts2, ts3;
+ int ts2, ts3;
- ts = MP.t;
- MP.t = t;
mp_multiply(x, &t1, &t2);
mp_add_integer(&t2, -1, &t2);
- MP.t = ts;
-
- /* TEMPORARILY REDUCE T */
- ts = MP.t;
- MP.t = (t + it0) / 2;
mp_multiply(&t1, &t2, &t2);
- MP.t = ts;
-
- ts = MP.t;
- MP.t = t;
mp_subtract(&t1, &t2, &t1);
- MP.t = ts;
- if (t >= MP.t)
+ if (t >= MP_T)
break;
/* FOLLOWING LOOP ALMOST DOUBLES T (POSSIBLE
* BECAUSE NEWTONS METHOD HAS 2ND ORDER CONVERGENCE).
*/
ts3 = t;
- t = MP.t;
+ t = MP_T;
do {
ts2 = t;
t = (t + it0) / 2;
} while (t > ts3);
- t = min(ts2, MP.t);
+ t = min(ts2, MP_T);
}
/* RETURN IF NEWTON ITERATION WAS CONVERGING */
- if (t2.sign != 0 && (t1.exponent - t2.exponent) << 1 < MP.t - it0) {
+ if (t2.sign != 0 && (t1.exponent - t2.exponent) << 1 < MP_T - it0) {
/* THE FOLLOWING MESSAGE MAY INDICATE THAT B**(T-1) IS TOO SMALL,
* OR THAT THE STARTING APPROXIMATION IS NOT ACCURATE ENOUGH.
*/
@@ -1984,11 +1605,6 @@ mp_reciprocal(const MPNumber *x, MPNumber *z)
/* MOVE RESULT TO Y AND RETURN AFTER RESTORING T */
mp_set_from_mp(&t1, z);
-
- /* RESTORE M AND CHECK FOR OVERFLOW (UNDERFLOW IMPOSSIBLE) */
- MP.m += -2;
- if (z->exponent > MP.m)
- mpovfl(z, "*** OVERFLOW OCCURRED IN MP_RECIPROCAL ***");
}
@@ -2009,9 +1625,6 @@ mp_root(const MPNumber *x, int n, MPNumber *z)
float rx;
MPNumber t1, t2;
- /* CHECK LEGALITY OF B, T, M AND MXR */
- mpchk();
-
if (n == 1) {
mp_set_from_mp(x, z);
return;
@@ -2019,34 +1632,31 @@ mp_root(const MPNumber *x, int n, MPNumber *z)
if (n == 0) {
mperr("*** N == 0 IN CALL TO MP_ROOT ***");
- z->sign = 0;
+ mp_set_from_integer(0, z);
return;
}
np = abs(n);
/* LOSS OF ACCURACY IF NP LARGE, SO ONLY ALLOW NP <= MAX (B, 64) */
- if (np > max(MP.b, 64)) {
+ if (np > max(MP_BASE, 64)) {
mperr("*** ABS(N) TOO LARGE IN CALL TO MP_ROOT ***");
- z->sign = 0;
+ mp_set_from_integer(0, z);
return;
}
/* LOOK AT SIGN OF X */
if (x->sign == 0) {
/* X == 0 HERE, RETURN 0 IF N POSITIVE, ERROR IF NEGATIVE */
- z->sign = 0;
- if (n > 0)
- return;
-
- mperr("*** X == 0 AND N NEGATIVE IN CALL TO MP_ROOT ***");
- z->sign = 0;
+ mp_set_from_integer(0, z);
+ if (n <= 0)
+ mperr("*** X == 0 AND N NEGATIVE IN CALL TO MP_ROOT ***");
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 ***");
- z->sign = 0;
+ mp_set_from_integer(0, z);
return;
}
@@ -2062,8 +1672,8 @@ mp_root(const MPNumber *x, int n, MPNumber *z)
}
/* USE SINGLE-PRECISION ROOT FOR FIRST APPROXIMATION */
- r__1 = exp(((float) (np * ex - x->exponent) * log((float) MP.b) -
- log((fabs(rx)))) / (float) np);
+ r__1 = exp(((float)(np * ex - x->exponent) * log((float)MP_BASE) -
+ log((fabs(rx)))) / (float)np);
mp_set_from_float(r__1, &t1);
/* SIGN OF APPROXIMATION SAME AS THAT OF X */
@@ -2074,57 +1684,45 @@ mp_root(const MPNumber *x, int n, MPNumber *z)
/* START WITH SMALL T TO SAVE TIME */
/* ENSURE THAT B**(T-1) >= 100 */
- if (MP.b < 10)
- t = it[MP.b - 1];
+ if (MP_BASE < 10)
+ t = it[MP_BASE - 1];
else
t = 3;
- if (t <= MP.t) {
+ if (t <= MP_T) {
/* IT0 IS A NECESSARY SAFETY FACTOR */
it0 = (t + 4) / 2;
/* MAIN ITERATION LOOP */
while(1) {
- int ts, ts2, ts3;
+ int ts2, ts3;
- ts = MP.t;
- MP.t = t;
- mp_pwr_integer(&t1, np, &t2);
+ mp_xpowy_integer(&t1, np, &t2);
mp_multiply(x, &t2, &t2);
mp_add_integer(&t2, -1, &t2);
- MP.t = ts;
-
- /* TEMPORARILY REDUCE T */
- ts = MP.t;
- MP.t = (t + it0) / 2;
mp_multiply(&t1, &t2, &t2);
mp_divide_integer(&t2, np, &t2);
- MP.t = ts;
-
- ts = MP.t;
- MP.t = t;
mp_subtract(&t1, &t2, &t1);
- MP.t = ts;
/* FOLLOWING LOOP ALMOST DOUBLES T (POSSIBLE BECAUSE
* NEWTONS METHOD HAS 2ND ORDER CONVERGENCE).
*/
- if (t >= MP.t)
+ if (t >= MP_T)
break;
ts3 = t;
- t = MP.t;
+ t = MP_T;
do {
ts2 = t;
t = (t + it0) / 2;
} while (t > ts3);
- t = min(ts2, MP.t);
+ t = min(ts2, MP_T);
}
/* NOW R(I2) IS X**(-1/NP)
* CHECK THAT NEWTON ITERATION WAS CONVERGING
*/
- if (t2.sign != 0 && (t1.exponent - t2.exponent) << 1 < MP.t - it0) {
+ if (t2.sign != 0 && (t1.exponent - t2.exponent) << 1 < MP_T - it0) {
/* THE FOLLOWING MESSAGE MAY INDICATE THAT B**(T-1) IS TOO SMALL,
* OR THAT THE INITIAL APPROXIMATION OBTAINED USING ALOG AND EXP
* IS NOT ACCURATE ENOUGH.
@@ -2134,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;
}
@@ -2154,14 +1752,12 @@ mp_sqrt(const MPNumber *x, MPNumber *z)
int i, i2, iy3;
MPNumber t;
- mpchk();
-
/* MP_ROOT NEEDS 4T+10 WORDS, BUT CAN OVERLAP SLIGHTLY. */
- i2 = MP.t * 3 + 9;
+ i2 = MP_T * 3 + 9;
if (x->sign < 0) {
mperr("*** X NEGATIVE IN CALL TO SUBROUTINE MP_SQRT ***");
} else if (x->sign == 0) {
- z->sign = 0;
+ mp_set_from_integer(0, z);
} else {
mp_root(x, -2, &t);
i = t.fraction[0];
@@ -2224,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);
@@ -2238,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 36a14c4..2e46492 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,25 +44,31 @@
/* Size of the multiple precision values */
#define MP_SIZE 1000
-/* Object for a high precision floating point number representation */
+/* Base for numbers */
+#define MP_BASE 10000
+
+/* Object for a high precision floating point number representation
+ *
+ * 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]; // Size MP.t?
+ 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
@@ -139,20 +146,17 @@ 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) */
-void mp_logarithm(int n, MPNumber *x, MPNumber *z);
+/* Sets z = log_n x */
+void mp_logarithm(int n, const MPNumber *x, MPNumber *z);
/* Sets z = Ï? */
void mp_get_pi(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 = e */
+void mp_get_eulers(MPNumber *z);
/* Sets z = nâ??x */
void mp_root(const MPNumber *x, int n, MPNumber *z);
@@ -163,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);
@@ -211,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]
