summaryrefslogtreecommitdiff
path: root/src
diff options
context:
space:
mode:
Diffstat (limited to 'src')
-rw-r--r--src/Makefile.am16
-rw-r--r--src/currency-manager.c13
-rw-r--r--src/currency.c1
-rw-r--r--src/financial.c68
-rw-r--r--src/lexer.c5
-rw-r--r--src/mate-calc-cmd.c3
-rw-r--r--src/mate-calc.c5
-rw-r--r--src/math-buttons.c20
-rw-r--r--src/math-converter.c10
-rw-r--r--src/math-equation.c33
-rw-r--r--src/math-variable-popup.c6
-rw-r--r--src/math-variables.c5
-rw-r--r--src/mp-binary.c31
-rw-r--r--src/mp-convert.c436
-rw-r--r--src/mp-private.h39
-rw-r--r--src/mp-serializer.c88
-rw-r--r--src/mp-trigonometric.c578
-rw-r--r--src/mp.c1962
-rw-r--r--src/mp.h63
-rw-r--r--src/parser.c7
-rw-r--r--src/parserfunc.c327
-rw-r--r--src/preferences.ui2
-rw-r--r--src/test-mp-equation.c5
-rw-r--r--src/test-mp.c27
-rw-r--r--src/unit-category.c11
-rw-r--r--src/unittest.c45
26 files changed, 877 insertions, 2929 deletions
diff --git a/src/Makefile.am b/src/Makefile.am
index 0555049..e9e9490 100644
--- a/src/Makefile.am
+++ b/src/Makefile.am
@@ -50,7 +50,6 @@ mate_calc_SOURCES = \
mp-enums.h \
mp-equation.c \
mp-equation.h \
- mp-private.h \
mp-serializer.c \
mp-serializer.h \
mp-trigonometric.c \
@@ -107,7 +106,10 @@ mate_calc_cmd_SOURCES = \
mate_calc_cmd_LDADD = \
$(MATE_CALC_CMD_LIBS) \
- -lm
+ -lm \
+ -lgmp \
+ -lmpfr \
+ -lmpc
test_mp_SOURCES = \
test-mp.c \
@@ -122,7 +124,10 @@ test_mp_SOURCES = \
test_mp_LDADD = \
$(MATE_CALC_CMD_LIBS) \
- -lm
+ -lm \
+ -lgmp \
+ -lmpfr \
+ -lmpc
test_mp_equation_SOURCES = \
test-mp-equation.c \
@@ -156,7 +161,10 @@ test_mp_equation_SOURCES = \
test_mp_equation_LDADD = \
$(MATE_CALC_CMD_LIBS) \
- -lm
+ -lm \
+ -lgmp \
+ -lmpfr \
+ -lmpc
CLEANFILES = \
mp-enums.c \
diff --git a/src/currency-manager.c b/src/currency-manager.c
index 40e87a7..1299522 100644
--- a/src/currency-manager.c
+++ b/src/currency-manager.c
@@ -394,7 +394,7 @@ load_imf_rates(CurrencyManager *manager)
}
if (name_map[name_index].name) {
Currency *c = currency_manager_get_currency(manager, name_map[name_index].symbol);
- MPNumber value;
+ MPNumber value = mp_new();
if (!c) {
g_debug ("Using IMF rate of %s for %s", tokens[value_index], name_map[name_index].symbol);
@@ -403,6 +403,7 @@ load_imf_rates(CurrencyManager *manager)
mp_set_from_string(tokens[value_index], 10, &value);
mp_reciprocal(&value, &value);
currency_set_value(c, &value);
+ mp_clear(&value);
}
else
g_warning("Unknown currency '%s'", tokens[0]);
@@ -435,7 +436,8 @@ set_ecb_rate(CurrencyManager *manager, xmlNodePtr node, Currency *eur_rate)
/* Use data if value and no rate currently defined */
if (name && value && !currency_manager_get_currency(manager, name)) {
Currency *c;
- MPNumber r, v;
+ MPNumber r = mp_new();
+ MPNumber v = mp_new();
g_debug ("Using ECB rate of %s for %s", value, name);
c = add_currency(manager, name);
@@ -443,6 +445,8 @@ set_ecb_rate(CurrencyManager *manager, xmlNodePtr node, Currency *eur_rate)
mp_set_from_mp(currency_get_value(eur_rate), &v);
mp_multiply(&v, &r, &v);
currency_set_value(c, &v);
+ mp_clear(&r);
+ mp_clear(&v);
}
if (name)
@@ -456,7 +460,8 @@ static void
set_ecb_fixed_rate(CurrencyManager *manager, const gchar *name, const gchar *value, Currency *eur_rate)
{
Currency *c;
- MPNumber r, v;
+ MPNumber r = mp_new();
+ MPNumber v = mp_new();
g_debug ("Using ECB fixed rate of %s for %s", value, name);
c = add_currency(manager, name);
@@ -464,6 +469,8 @@ set_ecb_fixed_rate(CurrencyManager *manager, const gchar *name, const gchar *val
mp_set_from_mp(currency_get_value(eur_rate), &v);
mp_divide(&v, &r, &v);
currency_set_value(c, &v);
+ mp_clear(&r);
+ mp_clear(&v);
}
diff --git a/src/currency.c b/src/currency.c
index 312e661..a6295ee 100644
--- a/src/currency.c
+++ b/src/currency.c
@@ -36,6 +36,7 @@ currency_new(const gchar *name,
currency->priv->name = g_strdup(name);
currency->priv->display_name = g_strdup(display_name);
currency->priv->symbol = g_strdup(symbol);
+ currency->priv->value = mp_new();
return currency;
}
diff --git a/src/financial.c b/src/financial.c
index fa0e33d..9552cc1 100644
--- a/src/financial.c
+++ b/src/financial.c
@@ -24,13 +24,18 @@ calc_ctrm(MathEquation *equation, MPNumber *t, MPNumber *pint, MPNumber *fv, MPN
*
* RESULT = log(fv / pv) / log(1 + pint)
*/
- MPNumber MP1, MP2, MP3, MP4;
+ MPNumber MP1 = mp_new();
+ MPNumber MP2 = mp_new();
+ MPNumber MP3 = mp_new();
+ MPNumber MP4 = mp_new();
mp_divide(fv, pv, &MP1);
mp_ln(&MP1, &MP2);
mp_add_integer(pint, 1, &MP3);
mp_ln(&MP3, &MP4);
mp_divide(&MP2, &MP4, t);
+
+ mp_clear(&MP1); mp_clear(&MP2); mp_clear(&MP3); mp_clear(&MP4);
}
@@ -54,10 +59,12 @@ calc_ddb(MathEquation *equation, MPNumber *t, MPNumber *cost, MPNumber *life, MP
int i;
int len;
- MPNumber MPbv, MP1, MP2;
+ MPNumber MPbv = mp_new();
+ MPNumber MP1 = mp_new();
+ MPNumber MP2 = mp_new();
mp_set_from_integer(0, &MPbv);
- len = mp_cast_to_int(period);
+ len = mp_to_integer(period);
for (i = 0; i < len; i++) {
mp_subtract(cost, &MPbv, &MP1);
mp_multiply_integer(&MP1, 2, &MP2);
@@ -70,6 +77,8 @@ calc_ddb(MathEquation *equation, MPNumber *t, MPNumber *cost, MPNumber *life, MP
math_equation_set_status (equation, _("Error: the number of periods must be positive"));
mp_set_from_integer(0, t);
}
+
+ mp_clear(&MPbv); mp_clear(&MP1); mp_clear(&MP2);
}
@@ -84,13 +93,18 @@ calc_fv(MathEquation *equation, MPNumber *t, MPNumber *pmt, MPNumber *pint, MPNu
* RESULT = pmt * (pow(1 + pint, n) - 1) / pint
*/
- MPNumber MP1, MP2, MP3, MP4;
+ MPNumber MP1 = mp_new();
+ MPNumber MP2 = mp_new();
+ MPNumber MP3 = mp_new();
+ MPNumber MP4 = mp_new();
mp_add_integer(pint, 1, &MP1);
mp_xpowy(&MP1, n, &MP2);
mp_add_integer(&MP2, -1, &MP3);
mp_multiply(pmt, &MP3, &MP4);
mp_divide(&MP4, pint, t);
+
+ mp_clear(&MP1); mp_clear(&MP2); mp_clear(&MP3); mp_clear(&MP4);
}
@@ -104,11 +118,14 @@ calc_gpm(MathEquation *equation, MPNumber *t, MPNumber *cost, MPNumber *margin)
* RESULT = cost / (1 - margin)
*/
- MPNumber MP1, MP2;
+ MPNumber MP1 = mp_new();
+ MPNumber MP2 = mp_new();
mp_set_from_integer(1, &MP1);
mp_subtract(&MP1, margin, &MP2);
mp_divide(cost, &MP2, t);
+
+ mp_clear(&MP1); mp_clear(&MP2);
}
@@ -123,7 +140,10 @@ calc_pmt(MathEquation *equation, MPNumber *t, MPNumber *prin, MPNumber *pint, MP
* RESULT = prin * (pint / (1 - pow(pint + 1, -1 * n)))
*/
- MPNumber MP1, MP2, MP3, MP4;
+ MPNumber MP1 = mp_new();
+ MPNumber MP2 = mp_new();
+ MPNumber MP3 = mp_new();
+ MPNumber MP4 = mp_new();
mp_add_integer(pint, 1, &MP1);
mp_multiply_integer(n, -1, &MP2);
@@ -132,6 +152,8 @@ calc_pmt(MathEquation *equation, MPNumber *t, MPNumber *prin, MPNumber *pint, MP
mp_add_integer(&MP4, 1, &MP1);
mp_divide(pint, &MP1, &MP2);
mp_multiply(prin, &MP2, t);
+
+ mp_clear(&MP1); mp_clear(&MP2); mp_clear(&MP3); mp_clear(&MP4);
}
@@ -146,7 +168,10 @@ calc_pv(MathEquation *equation, MPNumber *t, MPNumber *pmt, MPNumber *pint, MPNu
* RESULT = pmt * (1 - pow(1 + pint, -1 * n)) / pint
*/
- MPNumber MP1, MP2, MP3, MP4;
+ MPNumber MP1 = mp_new();
+ MPNumber MP2 = mp_new();
+ MPNumber MP3 = mp_new();
+ MPNumber MP4 = mp_new();
mp_add_integer(pint, 1, &MP1);
mp_multiply_integer(n, -1, &MP2);
@@ -155,6 +180,8 @@ calc_pv(MathEquation *equation, MPNumber *t, MPNumber *pmt, MPNumber *pint, MPNu
mp_add_integer(&MP4, 1, &MP1);
mp_divide(&MP1, pint, &MP2);
mp_multiply(pmt, &MP2, t);
+
+ mp_clear(&MP1); mp_clear(&MP2); mp_clear(&MP3); mp_clear(&MP4);
}
@@ -169,13 +196,18 @@ calc_rate(MathEquation *equation, MPNumber *t, MPNumber *fv, MPNumber *pv, MPNum
* RESULT = pow(fv / pv, 1 / n) - 1
*/
- MPNumber MP1, MP2, MP3, MP4;
+ MPNumber MP1 = mp_new();
+ MPNumber MP2 = mp_new();
+ MPNumber MP3 = mp_new();
+ MPNumber MP4 = mp_new();
mp_divide(fv, pv, &MP1);
mp_set_from_integer(1, &MP2);
mp_divide(&MP2, n, &MP3);
mp_xpowy(&MP1, &MP3, &MP4);
mp_add_integer(&MP4, -1, t);
+
+ mp_clear(&MP1); mp_clear(&MP2); mp_clear(&MP3); mp_clear(&MP4);
}
@@ -190,9 +222,10 @@ calc_sln(MathEquation *equation, MPNumber *t, MPNumber *cost, MPNumber *salvage,
* RESULT = (cost - salvage) / life
*/
- MPNumber MP1;
+ MPNumber MP1 = mp_new();
mp_subtract(cost, salvage, &MP1);
mp_divide(&MP1, life, t);
+ mp_clear(&MP1);
}
@@ -209,7 +242,10 @@ calc_syd(MathEquation *equation, MPNumber *t, MPNumber *cost, MPNumber *salvage,
* (life * (life + 1)) / 2
*/
- MPNumber MP1, MP2, MP3, MP4;
+ MPNumber MP1 = mp_new();
+ MPNumber MP2 = mp_new();
+ MPNumber MP3 = mp_new();
+ MPNumber MP4 = mp_new();
mp_subtract(life, period, &MP2);
mp_add_integer(&MP2, 1, &MP3);
@@ -220,6 +256,8 @@ calc_syd(MathEquation *equation, MPNumber *t, MPNumber *cost, MPNumber *salvage,
mp_divide(&MP3, &MP1, &MP2);
mp_subtract(cost, salvage, &MP1);
mp_multiply(&MP1, &MP2, t);
+
+ mp_clear(&MP1); mp_clear(&MP2); mp_clear(&MP3); mp_clear(&MP4);
}
@@ -234,7 +272,10 @@ calc_term(MathEquation *equation, MPNumber *t, MPNumber *pmt, MPNumber *fv, MPNu
* RESULT = log(1 + (fv * pint / pmt)) / log(1 + pint)
*/
- MPNumber MP1, MP2, MP3, MP4;
+ MPNumber MP1 = mp_new();
+ MPNumber MP2 = mp_new();
+ MPNumber MP3 = mp_new();
+ MPNumber MP4 = mp_new();
mp_add_integer(pint, 1, &MP1);
mp_ln(&MP1, &MP2);
@@ -243,13 +284,15 @@ calc_term(MathEquation *equation, MPNumber *t, MPNumber *pmt, MPNumber *fv, MPNu
mp_add_integer(&MP3, 1, &MP4);
mp_ln(&MP4, &MP1);
mp_divide(&MP1, &MP2, t);
+
+ mp_clear(&MP1); mp_clear(&MP2); mp_clear(&MP3); mp_clear(&MP4);
}
void
do_finc_expression(MathEquation *equation, int function, MPNumber *arg1, MPNumber *arg2, MPNumber *arg3, MPNumber *arg4)
{
- MPNumber result;
+ MPNumber result = mp_new();
switch (function) {
case FINC_CTRM_DIALOG:
calc_ctrm(equation, &result, arg1, arg2, arg3);
@@ -283,4 +326,5 @@ do_finc_expression(MathEquation *equation, int function, MPNumber *arg1, MPNumbe
break;
}
math_equation_set_number(equation, &result);
+ mp_clear(&result);
}
diff --git a/src/lexer.c b/src/lexer.c
index 176c773..10bfe36 100644
--- a/src/lexer.c
+++ b/src/lexer.c
@@ -30,12 +30,13 @@ l_check_if_function(LexerState* state)
static gboolean
l_check_if_number(LexerState* state)
{
- MPNumber tmp;
+ MPNumber tmp = mp_new();
int count = 0;
gchar* text = pl_get_marked_substring(state->prelexer);
if(mp_set_from_string(text, state->parent->options->base, &tmp) == 0)
{
free(text);
+ mp_clear(&tmp);
return TRUE;
}
else
@@ -46,6 +47,7 @@ l_check_if_number(LexerState* state)
if(mp_set_from_string(text, state->parent->options->base, &tmp) == 0)
{
free(text);
+ mp_clear(&tmp);
return TRUE;
}
free(text);
@@ -57,6 +59,7 @@ l_check_if_number(LexerState* state)
while(count--)
pl_get_next_token (state->prelexer);
free(text);
+ mp_clear(&tmp);
return FALSE;
}
}
diff --git a/src/mate-calc-cmd.c b/src/mate-calc-cmd.c
index b4fb301..cf8db3d 100644
--- a/src/mate-calc-cmd.c
+++ b/src/mate-calc-cmd.c
@@ -27,7 +27,7 @@ solve(const char *equation)
{
int ret;
MPEquationOptions options;
- MPNumber z;
+ MPNumber z = mp_new();
gchar *result_str = NULL;
memset(&options, 0, sizeof(options));
@@ -46,6 +46,7 @@ solve(const char *equation)
printf("%s\n", result_str);
}
g_free(result_str);
+ mp_clear(&z);
}
diff --git a/src/mate-calc.c b/src/mate-calc.c
index 0c2645c..7dc20d8 100644
--- a/src/mate-calc.c
+++ b/src/mate-calc.c
@@ -45,7 +45,7 @@ solve(const char *equation)
{
MPEquationOptions options;
MPErrorCode error;
- MPNumber result;
+ MPNumber result = mp_new();
char *result_str;
memset(&options, 0, sizeof(options));
@@ -57,15 +57,18 @@ solve(const char *equation)
error = mp_equation_parse(equation, &options, &result, NULL);
if(error == PARSER_ERR_MP) {
fprintf(stderr, "Error: %s\n", mp_get_error());
+ mp_clear(&result);
exit(1);
}
else if(error != 0) {
fprintf(stderr, "Error: %s\n", mp_error_code_to_string(error));
+ mp_clear(&result);
exit(1);
}
else {
result_str = mp_serializer_to_string(mp_serializer_new(MP_DISPLAY_FORMAT_AUTOMATIC, 10, 9), &result);
printf("%s\n", result_str);
+ mp_clear(&result);
exit(0);
}
}
diff --git a/src/math-buttons.c b/src/math-buttons.c
index 4e6e6da..69e993a 100644
--- a/src/math-buttons.c
+++ b/src/math-buttons.c
@@ -360,7 +360,7 @@ load_finc_dialogs(MathButtons *buttons)
static void
update_bit_panel(MathButtons *buttons)
{
- MPNumber x;
+ MPNumber x = mp_new();
gboolean enabled;
guint64 bits;
int i;
@@ -373,14 +373,17 @@ update_bit_panel(MathButtons *buttons)
enabled = math_equation_get_number(buttons->priv->equation, &x);
if (enabled) {
- MPNumber max, fraction;
+ MPNumber max = mp_new();
+ MPNumber fraction = mp_new();
mp_set_from_unsigned_integer(G_MAXUINT64, &max);
mp_fractional_part(&x, &fraction);
if (mp_is_negative(&x) || mp_is_greater_than(&x, &max) || !mp_is_zero(&fraction))
enabled = FALSE;
else
- bits = mp_cast_to_unsigned_int(&x);
+ bits = mp_to_unsigned_integer(&x);
+ mp_clear(&max);
+ mp_clear(&fraction);
}
gtk_widget_set_sensitive(buttons->priv->bit_panel, enabled);
@@ -422,6 +425,7 @@ update_bit_panel(MathButtons *buttons)
gtk_label_set_text(GTK_LABEL(buttons->priv->base_label), label->str);
g_string_free(label, TRUE);
+ mp_clear(&x);
}
@@ -1177,6 +1181,7 @@ finc_response_cb(GtkWidget *widget, gint response_id, MathButtons *buttons)
dialog = GPOINTER_TO_INT(g_object_get_data(G_OBJECT(widget), "finc_dialog"));
for (i = 0; i < 4; i++) {
+ arg[i] = mp_new();
if (finc_dialog_fields[dialog][i] == NULL) {
continue;
}
@@ -1187,6 +1192,10 @@ finc_response_cb(GtkWidget *widget, gint response_id, MathButtons *buttons)
gtk_widget_grab_focus(GET_WIDGET(buttons->priv->financial_ui, finc_dialog_fields[dialog][0]));
do_finc_expression(buttons->priv->equation, dialog, &arg[0], &arg[1], &arg[2], &arg[3]);
+
+ for (i = 0; i < 4; i++) {
+ mp_clear(&arg[i]);
+ }
}
@@ -1200,7 +1209,7 @@ character_code_dialog_response_cb(GtkWidget *dialog, gint response_id, MathButto
text = gtk_entry_get_text(GTK_ENTRY(buttons->priv->character_code_entry));
if (response_id == GTK_RESPONSE_OK) {
- MPNumber x;
+ MPNumber x = mp_new();
int i = 0;
mp_set_from_integer(0, &x);
@@ -1213,10 +1222,9 @@ character_code_dialog_response_cb(GtkWidget *dialog, gint response_id, MathButto
else
break;
}
-
math_equation_insert_number(buttons->priv->equation, &x);
+ mp_clear(&x);
}
-
gtk_widget_hide(dialog);
}
diff --git a/src/math-converter.c b/src/math-converter.c
index fa06a09..1d3e274 100644
--- a/src/math-converter.c
+++ b/src/math-converter.c
@@ -81,7 +81,8 @@ convert_equation(MathConverter *converter, const MPNumber *x, MPNumber *z)
static void
update_result_label(MathConverter *converter)
{
- MPNumber x, z;
+ MPNumber x = mp_new();
+ MPNumber z = mp_new();
gboolean enabled;
if (!converter->priv->result_label)
@@ -111,6 +112,8 @@ update_result_label(MathConverter *converter)
g_object_unref(source_unit);
g_object_unref(target_unit);
}
+ mp_clear(&x);
+ mp_clear(&z);
}
@@ -363,7 +366,8 @@ static void
swap_button_clicked_cb(GtkButton *button, MathConverter *converter)
{
Unit *from_unit, *to_unit;
- MPNumber x, z;
+ MPNumber x = mp_new();
+ MPNumber z = mp_new();
if (math_equation_get_number(converter->priv->equation, &x) &&
convert_equation(converter, &x, &z))
@@ -377,6 +381,8 @@ swap_button_clicked_cb(GtkButton *button, MathConverter *converter)
g_object_unref(from_unit);
g_object_unref(to_unit);
+ mp_clear(&x);
+ mp_clear(&z);
}
static void
diff --git a/src/math-equation.c b/src/math-equation.c
index 1d349e2..45d61cc 100644
--- a/src/math-equation.c
+++ b/src/math-equation.c
@@ -300,6 +300,7 @@ get_current_state(MathEquation *equation)
gint ans_start = -1, ans_end = -1;
state = g_malloc0(sizeof(MathEquationState));
+ state->ans = mp_new();
if (equation->priv->ans_start)
{
@@ -939,7 +940,7 @@ math_equation_get_answer(MathEquation *equation)
void
math_equation_store(MathEquation *equation, const gchar *name)
{
- MPNumber t;
+ MPNumber t = mp_new();
g_return_if_fail(equation != NULL);
g_return_if_fail(name != NULL);
@@ -948,6 +949,7 @@ math_equation_store(MathEquation *equation, const gchar *name)
math_equation_set_status(equation, _("No sane value to store"));
else
math_variables_set(equation->priv->variables, name, &t);
+ mp_clear(&t);
}
@@ -1148,7 +1150,6 @@ get_variable(const char *name, MPNumber *z, void *data)
else
result = 0;
}
-
free(lower_name);
return result;
@@ -1203,7 +1204,7 @@ math_equation_solve_real(gpointer data)
gint n_brackets = 0, result;
gchar *c, *text, *error_token;
GString *equation_text;
- MPNumber z;
+ MPNumber z = mp_new();
text = math_equation_get_equation(equation);
equation_text = g_string_new(text);
@@ -1227,6 +1228,7 @@ math_equation_solve_real(gpointer data)
switch (result) {
case PARSER_ERR_NONE:
solvedata->number_result = g_slice_new(MPNumber);
+ *solvedata->number_result = mp_new();
mp_set_from_mp(&z, solvedata->number_result);
break;
@@ -1267,6 +1269,7 @@ math_equation_solve_real(gpointer data)
break;
}
g_async_queue_push(equation->priv->queue, solvedata);
+ mp_clear(&z);
return NULL;
}
@@ -1302,6 +1305,7 @@ math_equation_look_for_answer(gpointer data)
}
else if (result->number_result != NULL) {
math_equation_set_number(equation, result->number_result);
+ mp_clear(result->number_result);
g_slice_free(MPNumber, result->number_result);
}
else if (result->text_result != NULL) {
@@ -1348,7 +1352,7 @@ math_equation_factorize_real(gpointer data)
{
GString *text;
GList *factors, *factor;
- MPNumber x;
+ MPNumber x = mp_new();
MathEquation *equation = MATH_EQUATION(data);
SolveData *result = g_slice_new0(SolveData);
@@ -1366,6 +1370,7 @@ math_equation_factorize_real(gpointer data)
g_string_append(text, temp);
if (factor->next)
g_string_append(text, "×");
+ mp_clear(n);
g_slice_free(MPNumber, n);
g_free(temp);
}
@@ -1374,6 +1379,7 @@ math_equation_factorize_real(gpointer data)
result->text_result = g_strndup(text->str, text->len);
g_async_queue_push(equation->priv->queue, result);
g_string_free(text, TRUE);
+ mp_clear(&x);
return NULL;
}
@@ -1382,7 +1388,7 @@ math_equation_factorize_real(gpointer data)
void
math_equation_factorize(MathEquation *equation)
{
- MPNumber x;
+ MPNumber x = mp_new();
g_return_if_fail(equation != NULL);
@@ -1393,9 +1399,10 @@ math_equation_factorize(MathEquation *equation)
if (!math_equation_get_number(equation, &x) || !mp_is_integer(&x)) {
/* Error displayed when trying to factorize a non-integer value */
math_equation_set_status(equation, _("Need an integer to factorize"));
+ mp_clear(&x);
return;
}
-
+ mp_clear(&x);
equation->priv->in_solve = true;
g_thread_new("", math_equation_factorize_real, equation);
@@ -1456,7 +1463,7 @@ math_equation_clear(MathEquation *equation)
void
math_equation_shift(MathEquation *equation, gint count)
{
- MPNumber z;
+ MPNumber z = mp_new();
g_return_if_fail(equation != NULL);
@@ -1470,13 +1477,14 @@ math_equation_shift(MathEquation *equation, gint count)
mp_shift(&z, count, &z);
math_equation_set_number(equation, &z);
+ mp_clear(&z);
}
void
math_equation_toggle_bit(MathEquation *equation, guint bit)
{
- MPNumber x;
+ MPNumber x = mp_new();
guint64 bits;
gboolean result;
@@ -1484,18 +1492,20 @@ math_equation_toggle_bit(MathEquation *equation, guint bit)
result = math_equation_get_number(equation, &x);
if (result) {
- MPNumber max;
+ MPNumber max = mp_new();
mp_set_from_unsigned_integer(G_MAXUINT64, &max);
if (mp_is_negative(&x) || mp_is_greater_than(&x, &max))
result = FALSE;
else
- bits = mp_cast_to_unsigned_int(&x);
+ bits = mp_to_unsigned_integer(&x);
+ mp_clear(&max);
}
if (!result) {
math_equation_set_status(equation,
/* Message displayed when cannot toggle bit in display*/
_("Displayed value not an integer"));
+ mp_clear(&x);
return;
}
@@ -1505,6 +1515,7 @@ math_equation_toggle_bit(MathEquation *equation, guint bit)
// FIXME: Only do this if in ans format, otherwise set text in same format as previous number
math_equation_set_number(equation, &x);
+ mp_clear(&x);
}
@@ -1807,6 +1818,7 @@ math_equation_class_init(MathEquationClass *klass)
param_types);
}
+
static void
pre_insert_text_cb(MathEquation *equation,
GtkTextIter *location,
@@ -1957,6 +1969,7 @@ math_equation_init(MathEquation *equation)
equation->priv->variables = math_variables_new();
+ equation->priv->state.ans = mp_new();
equation->priv->state.status = g_strdup("");
equation->priv->word_size = 32;
equation->priv->angle_units = MP_DEGREES;
diff --git a/src/math-variable-popup.c b/src/math-variable-popup.c
index fee5526..9c2d7ce 100644
--- a/src/math-variable-popup.c
+++ b/src/math-variable-popup.c
@@ -78,7 +78,7 @@ static void
add_variable_cb(GtkWidget *widget, MathVariablePopup *popup)
{
const gchar *name;
- MPNumber z;
+ MPNumber z = mp_new();
name = gtk_entry_get_text(GTK_ENTRY(popup->priv->variable_name_entry));
if (name[0] == '\0')
@@ -92,6 +92,7 @@ add_variable_cb(GtkWidget *widget, MathVariablePopup *popup)
g_warning("Can't add variable %s, the display is not a number", name);
gtk_widget_destroy(gtk_widget_get_toplevel(widget));
+ mp_clear(&z);
}
@@ -99,7 +100,7 @@ static void
save_variable_cb(GtkWidget *widget, MathVariablePopup *popup)
{
const gchar *name;
- MPNumber z;
+ MPNumber z = mp_new();
name = g_object_get_data(G_OBJECT(widget), "variable_name");
if (math_equation_get_number(popup->priv->equation, &z))
@@ -110,6 +111,7 @@ save_variable_cb(GtkWidget *widget, MathVariablePopup *popup)
g_warning("Can't save variable %s, the display is not a number", name);
gtk_widget_destroy(gtk_widget_get_toplevel(widget));
+ mp_clear(&z);
}
diff --git a/src/math-variables.c b/src/math-variables.c
index 225ca9c..b44b387 100644
--- a/src/math-variables.c
+++ b/src/math-variables.c
@@ -59,10 +59,14 @@ registers_load(MathVariables *variables)
value = g_strstrip(value);
t = g_malloc(sizeof(MPNumber));
+ *t = mp_new();
if (mp_set_from_string(value, 10, t) == 0)
g_hash_table_insert(variables->priv->registers, g_strdup(name), t);
else
+ {
+ mp_clear(t);
g_free(t);
+ }
}
fclose(f);
}
@@ -135,6 +139,7 @@ math_variables_set(MathVariables *variables, const char *name, const MPNumber *v
g_return_if_fail(value != NULL);
t = g_malloc(sizeof(MPNumber));
+ *t = mp_new();
mp_set_from_mp(value, t);
g_hash_table_insert(variables->priv->registers, g_strdup(name), t);
registers_save(variables);
diff --git a/src/mp-binary.c b/src/mp-binary.c
index 50483f8..800724c 100644
--- a/src/mp-binary.c
+++ b/src/mp-binary.c
@@ -12,7 +12,6 @@
#include <stdio.h>
#include "mp.h"
-#include "mp-private.h"
#include "mp-serializer.h"
// FIXME: Make dynamic
@@ -99,10 +98,15 @@ static int mp_bitwise_not(int v1, int dummy) { return v1 ^ 0xF; }
bool
mp_is_overflow (const MPNumber *x, int wordlen)
{
- MPNumber tmp1, tmp2;
+ bool is_overflow;
+ MPNumber tmp1 = mp_new();
+ MPNumber tmp2 = mp_new();
mp_set_from_integer(2, &tmp1);
mp_xpowy_integer(&tmp1, wordlen, &tmp2);
- return mp_is_greater_than (&tmp2, x);
+ is_overflow = mp_is_greater_than (&tmp2, x);
+ mp_clear(&tmp1);
+ mp_clear(&tmp2);
+ return is_overflow;
}
@@ -148,7 +152,9 @@ mp_xor(const MPNumber *x, const MPNumber *y, MPNumber *z)
void
mp_not(const MPNumber *x, int wordlen, MPNumber *z)
{
- MPNumber temp;
+ MPNumber temp = mp_new();
+
+ mp_set_from_integer(0, &temp);
if (!mp_is_positive_integer(x))
{
@@ -156,46 +162,45 @@ mp_not(const MPNumber *x, int wordlen, MPNumber *z)
mperr(_("Boolean NOT is only defined for positive integers"));
}
- mp_set_from_integer(0, &temp);
mp_bitwise(x, &temp, mp_bitwise_not, z, wordlen);
+ mp_clear(&temp);
}
void
mp_shift(const MPNumber *x, int count, MPNumber *z)
{
- int i;
- MPNumber multiplier;
-
- mp_set_from_integer(1, &multiplier);
-
if (!mp_is_integer(x)) {
/* Translators: Error displayed when bit shift attempted on non-integer values */
mperr(_("Shift is only possible on integer values"));
return;
}
+ MPNumber multiplier = mp_new();
+ mp_set_from_integer(1, &multiplier);
if (count >= 0) {
- for (i = 0; i < count; i++)
+ for (int i = 0; i < count; i++)
mp_multiply_integer(&multiplier, 2, &multiplier);
mp_multiply(x, &multiplier, z);
}
else {
- for (i = 0; i < -count; i++)
+ for (int i = 0; i < -count; i++)
mp_multiply_integer(&multiplier, 2, &multiplier);
mp_divide(x, &multiplier, z);
mp_floor(z, z);
}
+ mp_clear(&multiplier);
}
void
mp_ones_complement(const MPNumber *x, int wordlen, MPNumber *z)
{
- MPNumber t;
+ MPNumber t = mp_new();
mp_set_from_integer(0, &t);
mp_bitwise(x, &t, mp_bitwise_xor, z, wordlen);
mp_not(z, wordlen, z);
+ mp_clear(&t);
}
diff --git a/src/mp-convert.c b/src/mp-convert.c
index f41c874..1770853 100644
--- a/src/mp-convert.c
+++ b/src/mp-convert.c
@@ -17,226 +17,38 @@
#include <langinfo.h>
#include "mp.h"
-#include "mp-private.h"
void
mp_set_from_mp(const MPNumber *x, MPNumber *z)
{
- if (z != x)
- memcpy(z, x, sizeof(MPNumber));
-}
-
-
-void
-mp_set_from_float(float rx, MPNumber *z)
-{
- int i, k, ib, ie, tp;
- float rj;
-
- mp_set_from_integer(0, z);
-
- /* CHECK SIGN */
- if (rx < 0.0f) {
- z->sign = -1;
- rj = -(double)(rx);
- } else if (rx > 0.0f) {
- z->sign = 1;
- rj = rx;
- } else {
- /* IF RX = 0E0 RETURN 0 */
- mp_set_from_integer(0, z);
- return;
- }
-
- /* INCREASE IE AND DIVIDE RJ BY 16. */
- ie = 0;
- while (rj >= 1.0f) {
- ++ie;
- rj *= 0.0625f;
- }
- while (rj < 0.0625f) {
- --ie;
- rj *= 16.0f;
- }
-
- /* NOW RJ IS DY DIVIDED BY SUITABLE POWER OF 16.
- * SET EXPONENT TO 0
- */
- z->exponent = 0;
-
- /* CONVERSION LOOP (ASSUME SINGLE-PRECISION OPS. EXACT) */
- for (i = 0; i < MP_T + 4; i++) {
- rj *= (float) MP_BASE;
- z->fraction[i] = (int) rj;
- rj -= (float) z->fraction[i];
- }
-
- /* NORMALIZE RESULT */
- mp_normalize(z);
-
- /* Computing MAX */
- ib = max(MP_BASE * 7 * MP_BASE, 32767) / 16;
- tp = 1;
-
- /* NOW MULTIPLY BY 16**IE */
- if (ie < 0) {
- k = -ie;
- for (i = 1; i <= k; ++i) {
- tp <<= 4;
- if (tp <= ib && tp != MP_BASE && i < k)
- continue;
- mp_divide_integer(z, tp, z);
- tp = 1;
- }
- } else if (ie > 0) {
- for (i = 1; i <= ie; ++i) {
- tp <<= 4;
- if (tp <= ib && tp != MP_BASE && i < ie)
- continue;
- mp_multiply_integer(z, tp, z);
- tp = 1;
- }
- }
+ mpc_set(z->num, x->num, MPC_RNDNN);
}
void
mp_set_from_double(double dx, MPNumber *z)
{
- int i, k, ib, ie, tp;
- double dj;
-
- mp_set_from_integer(0, z);
-
- /* CHECK SIGN */
- if (dx < 0.0) {
- z->sign = -1;
- dj = -dx;
- } else if (dx > 0.0) {
- z->sign = 1;
- dj = dx;
- } else {
- mp_set_from_integer(0, z);
- return;
- }
-
- /* INCREASE IE AND DIVIDE DJ BY 16. */
- for (ie = 0; dj >= 1.0; ie++)
- dj *= 1.0/16.0;
-
- for ( ; dj < 1.0/16.0; ie--)
- dj *= 16.0;
-
- /* NOW DJ IS DY DIVIDED BY SUITABLE POWER OF 16
- * SET EXPONENT TO 0
- */
- z->exponent = 0;
-
- /* CONVERSION LOOP (ASSUME DOUBLE-PRECISION OPS. EXACT) */
- for (i = 0; i < MP_T + 4; i++) {
- dj *= (double) MP_BASE;
- z->fraction[i] = (int) dj;
- dj -= (double) z->fraction[i];
- }
-
- /* NORMALIZE RESULT */
- mp_normalize(z);
-
- /* Computing MAX */
- ib = max(MP_BASE * 7 * MP_BASE, 32767) / 16;
- tp = 1;
-
- /* NOW MULTIPLY BY 16**IE */
- if (ie < 0) {
- k = -ie;
- for (i = 1; i <= k; ++i) {
- tp <<= 4;
- if (tp <= ib && tp != MP_BASE && i < k)
- continue;
- mp_divide_integer(z, tp, z);
- tp = 1;
- }
- } else if (ie > 0) {
- for (i = 1; i <= ie; ++i) {
- tp <<= 4;
- if (tp <= ib && tp != MP_BASE && i < ie)
- continue;
- mp_multiply_integer(z, tp, z);
- tp = 1;
- }
- }
+ mpc_set_d(z->num, dx, MPC_RNDNN);
}
void
mp_set_from_integer(int64_t x, MPNumber *z)
{
- int i;
-
- memset(z, 0, sizeof(MPNumber));
-
- if (x == 0) {
- z->sign = 0;
- return;
- }
-
- if (x < 0) {
- x = -x;
- z->sign = -1;
- }
- else if (x > 0)
- z->sign = 1;
-
- while (x != 0) {
- z->fraction[z->exponent] = x % MP_BASE;
- z->exponent++;
- x /= MP_BASE;
- }
- for (i = 0; i < z->exponent / 2; i++) {
- int t = z->fraction[i];
- z->fraction[i] = z->fraction[z->exponent - i - 1];
- z->fraction[z->exponent - i - 1] = t;
- }
+ mpc_set_si(z->num, x, MPC_RNDNN);
}
void
mp_set_from_unsigned_integer(uint64_t x, MPNumber *z)
{
- int i;
-
- mp_set_from_integer(0, z);
-
- if (x == 0) {
- z->sign = 0;
- return;
- }
- z->sign = 1;
-
- while (x != 0) {
- z->fraction[z->exponent] = x % MP_BASE;
- x = x / MP_BASE;
- z->exponent++;
- }
- for (i = 0; i < z->exponent / 2; i++) {
- int t = z->fraction[i];
- z->fraction[i] = z->fraction[z->exponent - i - 1];
- z->fraction[z->exponent - i - 1] = t;
- }
+ mpc_set_ui(z->num, x, MPC_RNDNN);
}
void
mp_set_from_fraction(int64_t numerator, int64_t denominator, MPNumber *z)
{
- mp_gcd(&numerator, &denominator);
-
- if (denominator == 0) {
- mperr("*** J == 0 IN CALL TO MP_SET_FROM_FRACTION ***\n");
- mp_set_from_integer(0, z);
- return;
- }
-
if (denominator < 0) {
numerator = -numerator;
denominator = -denominator;
@@ -251,252 +63,54 @@ mp_set_from_fraction(int64_t numerator, int64_t denominator, MPNumber *z)
void
mp_set_from_polar(const MPNumber *r, MPAngleUnit unit, const MPNumber *theta, MPNumber *z)
{
- MPNumber x, y;
+ MPNumber x = mp_new();
+ MPNumber y = mp_new();
mp_cos(theta, unit, &x);
mp_multiply(&x, r, &x);
mp_sin(theta, unit, &y);
mp_multiply(&y, r, &y);
mp_set_from_complex(&x, &y, z);
-}
+ mp_clear(&x);
+ mp_clear(&y);
+}
void
mp_set_from_complex(const MPNumber *x, const MPNumber *y, MPNumber *z)
{
- /* NOTE: Do imaginary component first as z may be x or y */
- z->im_sign = y->sign;
- z->im_exponent = y->exponent;
- memcpy(z->im_fraction, y->fraction, sizeof(int) * MP_SIZE);
-
- z->sign = x->sign;
- z->exponent = x->exponent;
- if (z != x)
- memcpy(z->fraction, x->fraction, sizeof(int) * MP_SIZE);
+ mpc_set_fr_fr(z->num, mpc_realref(x->num), mpc_realref(y->num), MPC_RNDNN);
}
-
void
mp_set_from_random(MPNumber *z)
{
mp_set_from_double(drand48(), z);
}
-
int64_t
-mp_cast_to_int(const MPNumber *x)
+mp_to_integer(const MPNumber *x)
{
- int i;
- int64_t z = 0, v;
-
- /* |x| <= 1 */
- if (x->sign == 0 || x->exponent <= 0)
- return 0;
-
- /* Multiply digits together */
- for (i = 0; i < x->exponent; i++) {
- int64_t t;
-
- t = z;
- z = z * MP_BASE + x->fraction[i];
-
- /* Check for overflow */
- if (z <= t)
- return 0;
- }
-
- /* Validate result */
- v = z;
- for (i = x->exponent - 1; i >= 0; i--) {
- int64_t digit;
-
- /* Get last digit */
- digit = v - (v / MP_BASE) * MP_BASE;
- if (x->fraction[i] != digit)
- return 0;
-
- v /= MP_BASE;
- }
- if (v != 0)
- return 0;
-
- return x->sign * z;
+ return mpfr_get_si(mpc_realref(x->num), MPFR_RNDN);
}
uint64_t
-mp_cast_to_unsigned_int(const MPNumber *x)
-{
- int i;
- uint64_t z = 0, v;
-
- /* x <= 1 */
- if (x->sign <= 0 || x->exponent <= 0)
- return 0;
-
- /* Multiply digits together */
- for (i = 0; i < x->exponent; i++) {
- uint64_t t;
-
- t = z;
- z = z * MP_BASE + x->fraction[i];
-
- /* Check for overflow */
- if (z <= t)
- return 0;
- }
-
- /* Validate result */
- v = z;
- for (i = x->exponent - 1; i >= 0; i--) {
- uint64_t digit;
-
- /* Get last digit */
- digit = v - (v / MP_BASE) * MP_BASE;
- if (x->fraction[i] != digit)
- return 0;
-
- v /= MP_BASE;
- }
- if (v != 0)
- return 0;
-
- return z;
-}
-
-
-static double
-mppow_ri(float ap, int bp)
+mp_to_unsigned_integer(const MPNumber *x)
{
- double pow;
-
- if (bp == 0)
- return 1.0;
-
- if (bp < 0) {
- if (ap == 0)
- return 1.0;
- bp = -bp;
- ap = 1 / ap;
- }
-
- pow = 1.0;
- for (;;) {
- if (bp & 01)
- pow *= ap;
- if (bp >>= 1)
- ap *= ap;
- else
- break;
- }
-
- return pow;
+ return mpfr_get_ui(mpc_realref(x->num), MPFR_RNDN);
}
-
float
-mp_cast_to_float(const MPNumber *x)
-{
- int i;
- float rz = 0.0;
-
- if (mp_is_zero(x))
- return 0.0;
-
- for (i = 0; i < MP_T; i++) {
- rz = (float) MP_BASE * rz + (float)x->fraction[i];
-
- /* CHECK IF FULL SINGLE-PRECISION ACCURACY ATTAINED */
- if (rz + 1.0f <= rz)
- break;
- }
-
- /* NOW ALLOW FOR EXPONENT */
- rz *= mppow_ri((float) MP_BASE, x->exponent - i - 1);
-
- /* CHECK REASONABLENESS OF RESULT */
- /* LHS SHOULD BE <= 0.5, BUT ALLOW FOR SOME ERROR IN ALOG */
- if (rz <= (float)0. ||
- fabs((float) x->exponent - (log(rz) / log((float) MP_BASE) + (float).5)) > (float).6) {
- /* FOLLOWING MESSAGE INDICATES THAT X IS TOO LARGE OR SMALL -
- * TRY USING MPCMRE INSTEAD.
- */
- mperr("*** FLOATING-POINT OVER/UNDER-FLOW IN MP_CAST_TO_FLOAT ***\n");
- return 0.0;
- }
-
- if (x->sign < 0)
- rz = -(double)(rz);
-
- return rz;
-}
-
-
-static double
-mppow_di(double ap, int bp)
+mp_to_float(const MPNumber *x)
{
- double pow = 1.0;
-
- if (bp != 0) {
- if (bp < 0) {
- if (ap == 0) return(pow);
- bp = -bp;
- ap = 1/ap;
- }
- for (;;) {
- if (bp & 01) pow *= ap;
- if (bp >>= 1) ap *= ap;
- else break;
- }
- }
-
- return(pow);
+ return mpfr_get_flt(mpc_realref(x->num), MPFR_RNDN);
}
-
double
-mp_cast_to_double(const MPNumber *x)
+mp_to_double(const MPNumber *x)
{
- int i, tm = 0;
- double d__1, dz2, ret_val = 0.0;
-
- if (mp_is_zero(x))
- return 0.0;
-
- for (i = 0; i < MP_T; i++) {
- ret_val = (double) MP_BASE * ret_val + (double) x->fraction[i];
- tm = i;
-
- /* CHECK IF FULL DOUBLE-PRECISION ACCURACY ATTAINED */
- dz2 = ret_val + 1.0;
-
- /* TEST BELOW NOT ALWAYS EQUIVALENT TO - IF (DZ2.LE.DZ) GO TO 20,
- * FOR EXAMPLE ON CYBER 76.
- */
- if (dz2 - ret_val <= 0.0)
- break;
- }
-
- /* NOW ALLOW FOR EXPONENT */
- ret_val *= mppow_di((double) MP_BASE, x->exponent - tm - 1);
-
- /* CHECK REASONABLENESS OF RESULT. */
- /* LHS SHOULD BE .LE. 0.5 BUT ALLOW FOR SOME ERROR IN DLOG */
- if (ret_val <= 0. ||
- ((d__1 = (double) ((float) x->exponent) - (log(ret_val) / log((double)
- ((float) MP_BASE)) + .5), abs(d__1)) > .6)) {
- /* FOLLOWING MESSAGE INDICATES THAT X IS TOO LARGE OR SMALL -
- * TRY USING MPCMDE INSTEAD.
- */
- mperr("*** FLOATING-POINT OVER/UNDER-FLOW IN MP_CAST_TO_DOUBLE ***\n");
- return 0.0;
- }
- else
- {
- if (x->sign < 0)
- ret_val = -ret_val;
- return ret_val;
- }
+ return mpfr_get_d(mpc_realref(x->num), MPFR_RNDN);
}
static int
@@ -586,7 +200,6 @@ set_from_sexagesimal(const char *str, int length, MPNumber *z)
{
int degrees = 0, minutes = 0;
char seconds[length+1];
- MPNumber t;
int n_matched;
seconds[0] = '\0';
@@ -594,6 +207,7 @@ set_from_sexagesimal(const char *str, int length, MPNumber *z)
if (n_matched < 1)
return true;
+ MPNumber t = mp_new();
mp_set_from_integer(degrees, z);
if (n_matched > 1) {
mp_set_from_integer(minutes, &t);
@@ -605,6 +219,7 @@ set_from_sexagesimal(const char *str, int length, MPNumber *z)
mp_divide_integer(&t, 3600, &t);
mp_add(z, &t, z);
}
+ mp_clear(&t);
return false;
}
@@ -674,9 +289,10 @@ mp_set_from_string(const char *str, int default_base, MPNumber *z)
}
}
if (fractions[i] != NULL) {
- MPNumber fraction;
+ MPNumber fraction = mp_new();
mp_set_from_fraction(numerators[i], denominators[i], &fraction);
mp_add(z, &fraction, z);
+ mp_clear(&fraction);
}
if (*c == '.') {
@@ -686,7 +302,8 @@ mp_set_from_string(const char *str, int default_base, MPNumber *z)
/* Convert fractional part */
if (has_fraction) {
- MPNumber numerator, denominator;
+ MPNumber numerator = mp_new();
+ MPNumber denominator = mp_new();
mp_set_from_integer(0, &numerator);
mp_set_from_integer(1, &denominator);
@@ -697,6 +314,8 @@ mp_set_from_string(const char *str, int default_base, MPNumber *z)
}
mp_divide(&numerator, &denominator, &numerator);
mp_add(z, &numerator, z);
+ mp_clear(&numerator);
+ mp_clear(&denominator);
}
if (c != end) {
@@ -704,10 +323,11 @@ mp_set_from_string(const char *str, int default_base, MPNumber *z)
}
if (multiplier != 0) {
- MPNumber t;
+ MPNumber t = mp_new();
mp_set_from_integer(10, &t);
mp_xpowy_integer(&t, multiplier, &t);
mp_multiply(z, &t, z);
+ mp_clear(&t);
}
if (negate == 1)
diff --git a/src/mp-private.h b/src/mp-private.h
deleted file mode 100644
index 3211ea8..0000000
--- a/src/mp-private.h
+++ /dev/null
@@ -1,39 +0,0 @@
-/*
- * 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 the Free Software
- * Foundation, either version 2 of the License, or (at your option) any later
- * version. See http://www.gnu.org/copyleft/gpl.html the full text of the
- * license.
- */
-
-#ifndef MP_INTERNAL_H
-#define MP_INTERNAL_H
-
-#include <glib/gi18n.h>
-
-/* If we're not using GNU C, elide __attribute__ */
-#ifndef __GNUC__
-# define __attribute__(x) /*NOTHING*/
-#endif
-
-#define min(a, b) ((a) <= (b) ? (a) : (b))
-#define max(a, b) ((a) >= (b) ? (a) : (b))
-
-//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 100
-
-void mperr(const char *format, ...) __attribute__((format(printf, 1, 2)));
-void mp_gcd(int64_t *, int64_t *);
-void mp_normalize(MPNumber *);
-void convert_to_radians(const MPNumber *x, MPAngleUnit unit, MPNumber *z);
-void convert_from_radians(const MPNumber *x, MPAngleUnit unit, MPNumber *z);
-
-#endif /* MP_INTERNAL_H */
diff --git a/src/mp-serializer.c b/src/mp-serializer.c
index 4067016..58f9030 100644
--- a/src/mp-serializer.c
+++ b/src/mp-serializer.c
@@ -56,10 +56,13 @@ mp_serializer_new(MpDisplayFormat format, int base, int trailing_digits)
static void
-mp_cast_to_string_real(MpSerializer *serializer, const MPNumber *x, int base, gboolean force_sign, int *n_digits, GString *string)
+mp_to_string_real(MpSerializer *serializer, const MPNumber *x, int base, gboolean force_sign, int *n_digits, GString *string)
{
static gchar digits[] = "0123456789ABCDEF";
- MPNumber number, integer_component, fractional_component, temp;
+ MPNumber number = mp_new();
+ MPNumber integer_component = mp_new();
+ MPNumber fractional_component = mp_new();
+ MPNumber temp = mp_new();
int i, last_non_zero;
if (mp_is_negative(x))
@@ -86,7 +89,9 @@ mp_cast_to_string_real(MpSerializer *serializer, const MPNumber *x, int base, gb
mp_set_from_mp(&integer_component, &temp);
i = 0;
do {
- MPNumber t, t2, t3;
+ MPNumber t = mp_new();
+ MPNumber t2 = mp_new();
+ MPNumber t3 = mp_new();
int64_t d;
if (serializer->priv->base == 10 && serializer->priv->show_tsep && i == serializer->priv->tsep_count) {
@@ -101,11 +106,21 @@ mp_cast_to_string_real(MpSerializer *serializer, const MPNumber *x, int base, gb
mp_subtract(&temp, &t2, &t3);
- d = mp_cast_to_int(&t3);
- g_string_prepend_c(string, d < 16 ? digits[d] : '?');
+ d = mp_to_integer(&t3);
+ if (d < 16 && d >= 0)
+ g_string_prepend_c(string, digits[d]);
+ else
+ { /* FIXME: Should show an overflow error message in fixed mode */
+ g_string_assign(string, "0");
+ *n_digits = serializer->priv->leading_digits+1;
+ break;
+ }
(*n_digits)++;
mp_set_from_mp(&t, &temp);
+ mp_clear(&t);
+ mp_clear(&t2);
+ mp_clear(&t3);
} while (!mp_is_zero(&temp));
last_non_zero = string->len;
@@ -116,17 +131,18 @@ mp_cast_to_string_real(MpSerializer *serializer, const MPNumber *x, int base, gb
mp_set_from_mp(&fractional_component, &temp);
for (i = serializer->priv->trailing_digits; i > 0 && !mp_is_zero(&temp); i--) {
int d;
- MPNumber digit;
+ MPNumber digit = mp_new();
mp_multiply_integer(&temp, base, &temp);
mp_floor(&temp, &digit);
- d = mp_cast_to_int(&digit);
+ d = mp_to_integer(&digit);
g_string_append_c(string, digits[d]);
if(d != 0)
last_non_zero = string->len;
mp_subtract(&temp, &digit, &temp);
+ mp_clear(&digit);
}
/* Strip trailing zeroes */
@@ -157,24 +173,26 @@ mp_cast_to_string_real(MpSerializer *serializer, const MPNumber *x, int base, gb
b -= d * multiplier;
}
}
+ mp_clear(&number); mp_clear(&integer_component); mp_clear(&fractional_component);
+ mp_clear(&temp);
}
static gchar *
-mp_cast_to_string(MpSerializer *serializer, const MPNumber *x, int *n_digits)
+mp_to_string(MpSerializer *serializer, const MPNumber *x, int *n_digits)
{
GString *string;
- MPNumber x_real;
+ MPNumber x_real = mp_new();
gchar *result;
string = g_string_sized_new(1024);
mp_real_component(x, &x_real);
- mp_cast_to_string_real(serializer, &x_real, serializer->priv->base, FALSE, n_digits, string);
+ mp_to_string_real(serializer, &x_real, serializer->priv->base, FALSE, n_digits, string);
if (mp_is_complex(x)) {
GString *s;
gboolean force_sign = TRUE;
- MPNumber x_im;
+ MPNumber x_im = mp_new();
int n_complex_digits = 0;
mp_imaginary_component(x, &x_im);
@@ -185,7 +203,7 @@ mp_cast_to_string(MpSerializer *serializer, const MPNumber *x, int *n_digits)
}
s = g_string_sized_new(1024);
- mp_cast_to_string_real(serializer, &x_im, 10, force_sign, &n_complex_digits, s);
+ mp_to_string_real(serializer, &x_im, 10, force_sign, &n_complex_digits, s);
if (n_complex_digits > *n_digits)
*n_digits = n_complex_digits;
if (strcmp(s->str, "0") == 0 || strcmp(s->str, "+0") == 0 || strcmp(s->str, "−0") == 0) {
@@ -209,20 +227,28 @@ mp_cast_to_string(MpSerializer *serializer, const MPNumber *x, int *n_digits)
g_string_append(string, "i");
}
g_string_free(s, TRUE);
+ mp_clear(&x_im);
}
result = g_strndup(string->str, string->len + 1);
g_string_free(string, TRUE);
+ mp_clear(&x_real);
return result;
}
static int
-mp_cast_to_exponential_string_real(MpSerializer *serializer, const MPNumber *x, GString *string, gboolean eng_format, int *n_digits)
+mp_to_exponential_string_real(MpSerializer *serializer, const MPNumber *x, GString *string, gboolean eng_format, int *n_digits)
{
gchar *fixed;
- MPNumber t, z, base, base3, base10, base10inv, mantissa;
+ MPNumber t = mp_new();
+ MPNumber z = mp_new();
+ MPNumber base = mp_new();
+ MPNumber base3 = mp_new();
+ MPNumber base10 = mp_new();
+ MPNumber base10inv = mp_new();
+ MPNumber mantissa = mp_new();
int exponent = 0;
mp_abs(x, &z);
@@ -260,9 +286,11 @@ mp_cast_to_exponential_string_real(MpSerializer *serializer, const MPNumber *x,
}
}
- fixed = mp_cast_to_string(serializer, &mantissa, n_digits);
+ fixed = mp_to_string(serializer, &mantissa, n_digits);
g_string_append(string, fixed);
g_free(fixed);
+ mp_clear(&t); mp_clear(&z); mp_clear(&base); mp_clear(&base3);
+ mp_clear(&base10); mp_clear(&base10inv); mp_clear(&mantissa);
return exponent;
}
@@ -291,22 +319,22 @@ append_exponent(GString *string, int exponent)
static gchar *
-mp_cast_to_exponential_string(MpSerializer *serializer, const MPNumber *x, gboolean eng_format, int *n_digits)
+mp_to_exponential_string(MpSerializer *serializer, const MPNumber *x, gboolean eng_format, int *n_digits)
{
GString *string;
- MPNumber x_real;
+ MPNumber x_real = mp_new();
gchar *result;
int exponent;
string = g_string_sized_new(1024);
mp_real_component(x, &x_real);
- exponent = mp_cast_to_exponential_string_real(serializer, &x_real, string, eng_format, n_digits);
+ exponent = mp_to_exponential_string_real(serializer, &x_real, string, eng_format, n_digits);
append_exponent(string, exponent);
if (mp_is_complex(x)) {
GString *s;
- MPNumber x_im;
+ MPNumber x_im = mp_new();
int n_complex_digits = 0;
mp_imaginary_component(x, &x_im);
@@ -317,7 +345,7 @@ mp_cast_to_exponential_string(MpSerializer *serializer, const MPNumber *x, gbool
g_string_append(string, "+");
s = g_string_sized_new(1024);
- exponent = mp_cast_to_exponential_string_real(serializer, &x_im, s, eng_format, &n_complex_digits);
+ exponent = mp_to_exponential_string_real(serializer, &x_im, s, eng_format, &n_complex_digits);
if (n_complex_digits > *n_digits)
*n_digits = n_complex_digits;
if (strcmp(s->str, "0") == 0 || strcmp(s->str, "+0") == 0 || strcmp(s->str, "−0") == 0) {
@@ -341,11 +369,13 @@ mp_cast_to_exponential_string(MpSerializer *serializer, const MPNumber *x, gbool
g_string_append(string, "i");
}
g_string_free(s, TRUE);
+ mp_clear(&x_im);
append_exponent(string, exponent);
}
result = g_strndup(string->str, string->len + 1);
g_string_free(string, TRUE);
+ mp_clear(&x_real);
return result;
}
@@ -354,31 +384,25 @@ mp_cast_to_exponential_string(MpSerializer *serializer, const MPNumber *x, gbool
gchar *
mp_serializer_to_string(MpSerializer *serializer, const MPNumber *x)
{
- MPNumber cmp, xcmp;
gchar *s0;
int n_digits = 0;
- mp_set_from_integer(10, &cmp);
- mp_xpowy_integer(&cmp, -(serializer->priv->trailing_digits), &cmp);
- mp_real_component(x, &xcmp);
- mp_abs(&xcmp, &xcmp);
switch(serializer->priv->format) {
default:
case MP_DISPLAY_FORMAT_AUTOMATIC:
- s0 = mp_cast_to_string(serializer, x, &n_digits);
- if (n_digits <= serializer->priv->leading_digits &&
- mp_is_greater_equal(&xcmp, &cmp))
+ s0 = mp_to_string(serializer, x, &n_digits);
+ if (n_digits <= serializer->priv->leading_digits)
return s0;
else {
g_free (s0);
- return mp_cast_to_exponential_string(serializer, x, FALSE, &n_digits);
+ return mp_to_exponential_string(serializer, x, FALSE, &n_digits);
}
break;
case MP_DISPLAY_FORMAT_FIXED:
- return mp_cast_to_string(serializer, x, &n_digits);
+ return mp_to_string(serializer, x, &n_digits);
case MP_DISPLAY_FORMAT_SCIENTIFIC:
- return mp_cast_to_exponential_string(serializer, x, FALSE, &n_digits);
+ return mp_to_exponential_string(serializer, x, FALSE, &n_digits);
case MP_DISPLAY_FORMAT_ENGINEERING:
- return mp_cast_to_exponential_string(serializer, x, TRUE, &n_digits);
+ return mp_to_exponential_string(serializer, x, TRUE, &n_digits);
}
}
diff --git a/src/mp-trigonometric.c b/src/mp-trigonometric.c
index 72aba0c..ca19ca4 100644
--- a/src/mp-trigonometric.c
+++ b/src/mp-trigonometric.c
@@ -15,614 +15,258 @@
#include <libintl.h>
#include "mp.h"
-#include "mp-private.h"
-
-static MPNumber pi;
-static gboolean have_pi = FALSE;
-
-static int
-mp_compare_mp_to_int(const MPNumber *x, int i)
-{
- MPNumber t;
- mp_set_from_integer(i, &t);
- return mp_compare_mp_to_mp(x, &t);
-}
-
/* Convert x to radians */
void
convert_to_radians(const MPNumber *x, MPAngleUnit unit, MPNumber *z)
{
- MPNumber t1, t2;
+ int i;
switch(unit) {
default:
case MP_RADIANS:
mp_set_from_mp(x, z);
- break;
+ return;
case MP_DEGREES:
- mp_get_pi(&t1);
- mp_multiply(x, &t1, &t2);
- mp_divide_integer(&t2, 180, z);
+ i = 180;
break;
case MP_GRADIANS:
- mp_get_pi(&t1);
- mp_multiply(x, &t1, &t2);
- mp_divide_integer(&t2, 200, z);
+ i = 200;
break;
}
+ mpfr_t scale;
+ mpfr_init2(scale, PRECISION);
+ mpfr_const_pi(scale, MPFR_RNDN);
+ mpfr_div_si(scale, scale, i, MPFR_RNDN);
+ mpc_mul_fr(z->num, x->num, scale, MPC_RNDNN);
+ mpfr_clear(scale);
}
-
-void
-mp_get_pi(MPNumber *z)
-{
- if (mp_is_zero(&pi)) {
- mp_set_from_string("3.1415926535897932384626433832795028841971693993751058209749445923078164062862089986280348253421170679", 10, &pi);
- have_pi = TRUE;
- }
- mp_set_from_mp(&pi, z);
-}
-
-
void
convert_from_radians(const MPNumber *x, MPAngleUnit unit, MPNumber *z)
{
- MPNumber t1, t2;
+ int i;
- switch (unit) {
+ switch(unit) {
default:
case MP_RADIANS:
mp_set_from_mp(x, z);
- break;
+ return;
case MP_DEGREES:
- mp_multiply_integer(x, 180, &t2);
- mp_get_pi(&t1);
- mp_divide(&t2, &t1, z);
+ i = 180;
break;
case MP_GRADIANS:
- mp_multiply_integer(x, 200, &t2);
- mp_get_pi(&t1);
- mp_divide(&t2, &t1, z);
+ i = 200;
break;
}
+ mpfr_t scale;
+ mpfr_init2(scale, PRECISION);
+ mpfr_const_pi(scale, MPFR_RNDN);
+ mpfr_si_div(scale, i, scale, MPFR_RNDN);
+ mpc_mul_fr(z->num, x->num, scale, MPC_RNDNN);
+ mpfr_clear(scale);
}
-/* z = sin(x) -1 >= x >= 1, do_sin = 1
- * z = cos(x) -1 >= x >= 1, do_sin = 0
- */
-static void
-mpsin1(const MPNumber *x, MPNumber *z, int do_sin)
-{
- int i, b2;
- MPNumber t1, t2;
-
- /* sin(0) = 0, cos(0) = 1 */
- if (mp_is_zero(x)) {
- 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 ***");
- }
-
- 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;
- }
-
- /* Taylor series */
- /* POWER SERIES LOOP. REDUCE T IF POSSIBLE */
- b2 = 2 * max(MP_BASE, 64);
- do {
- if (MP_T + t1.exponent <= 0)
- break;
-
- /* IF I*(I+1) IS NOT REPRESENTABLE AS AN INTEGER, THE FOLLOWING
- * DIVISION BY I*(I+1) HAS TO BE SPLIT UP.
- */
- mp_multiply(&t2, &t1, &t1);
- if (i > b2) {
- mp_divide_integer(&t1, -i, &t1);
- mp_divide_integer(&t1, i + 1, &t1);
- } else {
- mp_divide_integer(&t1, -i * (i + 1), &t1);
- }
- mp_add(&t1, z, z);
-
- i += 2;
- } while (t1.sign != 0);
-
- if (do_sin == 0)
- mp_add_integer(z, 1, z);
-}
-
-
-static void
-mp_sin_real(const MPNumber *x, MPAngleUnit unit, MPNumber *z)
-{
- int xs;
- MPNumber x_radians;
-
- /* sin(0) = 0 */
- if (mp_is_zero(x)) {
- mp_set_from_integer(0, z);
- return;
- }
-
- convert_to_radians(x, unit, &x_radians);
-
- xs = x_radians.sign;
- mp_abs(&x_radians, &x_radians);
-
- /* USE MPSIN1 IF ABS(X) <= 1 */
- if (mp_compare_mp_to_int(&x_radians, 1) <= 0) {
- mpsin1(&x_radians, z, 1);
- }
- /* FIND ABS(X) MODULO 2PI */
- else {
- 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(&x_radians, -1, 2, &x_radians);
- xs = -xs * x_radians.sign;
- if (xs == 0) {
- mp_set_from_integer(0, z);
- return;
- }
-
- x_radians.sign = 1;
- mp_multiply_integer(&x_radians, 4, &x_radians);
-
- /* IF NOT LESS THAN 1, SUBTRACT FROM 2 */
- if (x_radians.exponent > 0)
- mp_add_integer(&x_radians, -2, &x_radians);
-
- if (mp_is_zero(&x_radians)) {
- mp_set_from_integer(0, z);
- return;
- }
-
- 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 (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(&x_radians, z, &x_radians);
- mpsin1(&x_radians, z, 1);
- }
- }
-
- z->sign = xs;
-}
-
-
-static void
-mp_cos_real(const MPNumber *x, MPAngleUnit unit, MPNumber *z)
+void
+mp_get_pi (MPNumber *z)
{
- /* cos(0) = 1 */
- if (mp_is_zero(x)) {
- mp_set_from_integer(1, z);
- return;
- }
-
- convert_to_radians(x, unit, z);
-
- /* Use power series if |x| <= 1 */
- mp_abs(z, z);
- if (mp_compare_mp_to_int(z, 1) <= 0) {
- mpsin1(z, z, 0);
- } else {
- MPNumber t;
-
- /* cos(x) = sin(π/2 - |x|) */
- mp_get_pi(&t);
- mp_divide_integer(&t, 2, &t);
- mp_subtract(&t, z, z);
- mp_sin(z, MP_RADIANS, z);
- }
+ mpfr_const_pi(mpc_realref(z->num), MPFR_RNDN);
+ mpfr_set_zero(mpc_imagref(z->num), 0);
}
void
mp_sin(const MPNumber *x, MPAngleUnit unit, MPNumber *z)
{
- if (mp_is_complex(x)) {
- MPNumber x_real, x_im, z_real, z_im, t;
-
- mp_real_component(x, &x_real);
- mp_imaginary_component(x, &x_im);
-
- mp_sin_real(&x_real, unit, &z_real);
- mp_cosh(&x_im, &t);
- mp_multiply(&z_real, &t, &z_real);
-
- mp_cos_real(&x_real, unit, &z_im);
- mp_sinh(&x_im, &t);
- mp_multiply(&z_im, &t, &z_im);
-
- mp_set_from_complex(&z_real, &z_im, z);
- }
+ if (mp_is_complex(x))
+ mp_set_from_mp(x, z);
else
- mp_sin_real(x, unit, z);
+ convert_to_radians(x, unit, z);
+ mpc_sin(z->num, z->num, MPC_RNDNN);
}
void
mp_cos(const MPNumber *x, MPAngleUnit unit, MPNumber *z)
{
- if (mp_is_complex(x)) {
- MPNumber x_real, x_im, z_real, z_im, t;
-
- mp_real_component(x, &x_real);
- mp_imaginary_component(x, &x_im);
-
- mp_cos_real(&x_real, unit, &z_real);
- mp_cosh(&x_im, &t);
- mp_multiply(&z_real, &t, &z_real);
-
- mp_sin_real(&x_real, unit, &z_im);
- mp_sinh(&x_im, &t);
- mp_multiply(&z_im, &t, &z_im);
- mp_invert_sign(&z_im, &z_im);
-
- mp_set_from_complex(&z_real, &z_im, z);
- }
+ if (mp_is_complex(x))
+ mp_set_from_mp(x, z);
else
- mp_cos_real(x, unit, z);
+ convert_to_radians(x, unit, z);
+ mpc_cos(z->num, z->num, MPC_RNDNN);
}
void
mp_tan(const MPNumber *x, MPAngleUnit unit, MPNumber *z)
{
- MPNumber cos_x, sin_x;
+ MPNumber x_radians = mp_new();
+ MPNumber pi = mp_new();
+ MPNumber t1 = mp_new();
+
+ convert_to_radians(x, unit, &x_radians);
+ mp_get_pi(&pi);
+ mp_divide_integer(&pi, 2, &t1);
+ mp_subtract(&x_radians, &t1, &t1);
+ mp_divide(&t1, &pi, &t1);
- /* Check for undefined values */
- mp_cos(x, unit, &cos_x);
- if (mp_is_zero(&cos_x)) {
+ if (mp_is_integer(&t1)) {
/* Translators: Error displayed when tangent value is undefined */
mperr(_("Tangent is undefined for angles that are multiples of π (180°) from π∕2 (90°)"));
mp_set_from_integer(0, z);
return;
}
- /* tan(x) = sin(x) / cos(x) */
- mp_sin(x, unit, &sin_x);
- mp_divide(&sin_x, &cos_x, z);
+ if (mp_is_complex(x))
+ mp_set_from_mp(x, z);
+ else
+ mp_set_from_mp(&x_radians, z);
+ mpc_tan(z->num, z->num, MPC_RNDNN);
+ mp_clear(&x_radians);
+ mp_clear(&pi);
+ mp_clear(&t1);
}
void
mp_asin(const MPNumber *x, MPAngleUnit unit, MPNumber *z)
{
- MPNumber t1, t2;
+ MPNumber x_max = mp_new();
+ MPNumber x_min = mp_new();
+ mp_set_from_integer(1, &x_max);
+ mp_set_from_integer(-1, &x_min);
- /* asin⁻¹(0) = 0 */
- if (mp_is_zero(x)) {
+ if (mp_compare(x, &x_max) > 0 || mp_compare(x, &x_min) < 0)
+ {
+ /* Translators: Error displayed when inverse sine value is undefined */
+ mperr(_("Inverse sine is undefined for values outside [-1, 1]"));
mp_set_from_integer(0, z);
return;
}
-
- /* sin⁻¹(x) = tan⁻¹(x / √(1 - x²)), |x| < 1 */
- if (x->exponent <= 0) {
- mp_set_from_integer(1, &t1);
- mp_set_from_mp(&t1, &t2);
- mp_subtract(&t1, x, &t1);
- mp_add(&t2, x, &t2);
- mp_multiply(&t1, &t2, &t2);
- mp_root(&t2, -2, &t2);
- mp_multiply(x, &t2, z);
- mp_atan(z, unit, z);
- return;
- }
-
- /* sin⁻¹(1) = π/2, sin⁻¹(-1) = -π/2 */
- mp_set_from_integer(x->sign, &t2);
- if (mp_is_equal(x, &t2)) {
- mp_get_pi(z);
- mp_divide_integer(z, 2 * t2.sign, z);
+ mpc_asin(z->num, x->num, MPC_RNDNN);
+ if (!mp_is_complex(z))
convert_from_radians(z, unit, z);
- return;
- }
-
- /* Translators: Error displayed when inverse sine value is undefined */
- mperr(_("Inverse sine is undefined for values outside [-1, 1]"));
- mp_set_from_integer(0, z);
+ mp_clear(&x_max);
+ mp_clear(&x_min);
}
void
mp_acos(const MPNumber *x, MPAngleUnit unit, MPNumber *z)
{
- MPNumber t1, t2;
- MPNumber MPn1, pi, MPy;
-
- mp_get_pi(&pi);
- mp_set_from_integer(1, &t1);
- mp_set_from_integer(-1, &MPn1);
+ MPNumber x_max = mp_new();
+ MPNumber x_min = mp_new();
+ mp_set_from_integer(1, &x_max);
+ mp_set_from_integer(-1, &x_min);
- if (mp_is_greater_than(x, &t1) || mp_is_less_than(x, &MPn1)) {
- /* Translators: Error displayed when inverse cosine value is undefined */
+ if (mp_compare(x, &x_max) > 0 || mp_compare(x, &x_min) < 0)
+ {
+ /* Translators: Error displayed when inverse sine value is undefined */
mperr(_("Inverse cosine is undefined for values outside [-1, 1]"));
mp_set_from_integer(0, z);
- } else if (mp_is_zero(x)) {
- 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(&pi, z);
- } else {
- /* cos⁻¹(x) = tan⁻¹(√(1 - x²) / x) */
- 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, &pi, z);
- }
+ return;
}
-
- convert_from_radians(z, unit, z);
+ mpc_acos(z->num, x->num, MPC_RNDNN);
+ if (!mp_is_complex(z))
+ convert_from_radians(z, unit, z);
+ mp_clear(&x_max);
+ mp_clear(&x_min);
}
void
mp_atan(const MPNumber *x, MPAngleUnit unit, MPNumber *z)
{
- int i, q;
- float rx = 0.0;
- MPNumber t1, t2;
+ MPNumber i = mp_new();
+ MPNumber minus_i = mp_new();
+ mpc_set_si_si(i.num, 0, 1, MPC_RNDNN);
+ mpc_set_si_si(minus_i.num, 0, -1, MPC_RNDNN);
- if (mp_is_zero(x)) {
+ /* Check x != i and x != -i */
+ if (mp_is_equal(x, &i) || mp_is_equal(x, &minus_i))
+ {
+ /* Translators: Error displayed when inverse sine value is undefined */
+ mperr(_("Arctangent function is undefined for values i and -i"));
mp_set_from_integer(0, z);
return;
}
-
- mp_set_from_mp(x, &t2);
- if (abs(x->exponent) <= 2)
- rx = mp_cast_to_float(x);
-
- /* REDUCE ARGUMENT IF NECESSARY BEFORE USING SERIES */
- q = 1;
- while (t2.exponent >= 0)
- {
- if (t2.exponent == 0 && 2 * (t2.fraction[0] + 1) <= MP_BASE)
- break;
-
- q *= 2;
-
- /* t = t / (√(t² + 1) + 1) */
- mp_multiply(&t2, &t2, z);
- mp_add_integer(z, 1, z);
- mp_sqrt(z, z);
- mp_add_integer(z, 1, z);
- mp_divide(&t2, z, &t2);
- }
-
- /* USE POWER SERIES NOW ARGUMENT IN (-0.5, 0.5) */
- mp_set_from_mp(&t2, z);
- mp_multiply(&t2, &t2, &t1);
-
- /* SERIES LOOP. REDUCE T IF POSSIBLE. */
- for (i = 1; ; i += 2) {
- if (MP_T + 2 + t2.exponent <= 1)
- break;
-
- mp_multiply(&t2, &t1, &t2);
- mp_multiply_fraction(&t2, -i, i + 2, &t2);
-
- mp_add(z, &t2, z);
- if (mp_is_zero(&t2))
- break;
- }
-
- /* CORRECT FOR ARGUMENT REDUCTION */
- mp_multiply_integer(z, q, z);
-
- /* CHECK THAT RELATIVE ERROR LESS THAN 0.01 UNLESS EXPONENT
- * OF X IS LARGE (WHEN ATAN MIGHT NOT WORK)
- */
- if (abs(x->exponent) <= 2) {
- float ry = mp_cast_to_float(z);
- /* THE FOLLOWING MESSAGE MAY INDICATE THAT B**(T-1) IS TOO SMALL. */
- if (fabs(ry - atan(rx)) >= fabs(ry) * 0.01)
- mperr("*** ERROR OCCURRED IN MP_ATAN, RESULT INCORRECT ***");
- }
-
- convert_from_radians(z, unit, z);
+ mpc_atan(z->num, x->num, MPC_RNDNN);
+ if (!mp_is_complex(z))
+ convert_from_radians(z, unit, z);
+ mp_clear(&i);
+ mp_clear(&minus_i);
}
void
mp_sinh(const MPNumber *x, MPNumber *z)
{
- MPNumber abs_x;
-
- /* sinh(0) = 0 */
- if (mp_is_zero(x)) {
- mp_set_from_integer(0, z);
- return;
- }
-
- /* WORK WITH ABS(X) */
- mp_abs(x, &abs_x);
-
- /* 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;
-
- /* 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);
+ mpc_sinh(z->num, x->num, MPC_RNDNN);
}
void
mp_cosh(const MPNumber *x, MPNumber *z)
{
- MPNumber t;
-
- /* cosh(0) = 1 */
- if (mp_is_zero(x)) {
- mp_set_from_integer(1, z);
- return;
- }
-
- /* cosh(x) = (e^x + e^-x) / 2 */
- mp_abs(x, &t);
- mp_epowy(&t, &t);
- mp_reciprocal(&t, z);
- mp_add(&t, z, z);
- mp_divide_integer(z, 2, z);
+ mpc_cosh(z->num, x->num, MPC_RNDNN);
}
void
mp_tanh(const MPNumber *x, MPNumber *z)
{
- float r__1;
- MPNumber t;
-
- /* tanh(0) = 0 */
- if (mp_is_zero(x)) {
- mp_set_from_integer(0, z);
- return;
- }
-
- mp_abs(x, &t);
-
- /* SEE IF ABS(X) SO LARGE THAT RESULT IS +-1 */
- r__1 = (float) MP_T * 0.5 * log((float) MP_BASE);
- mp_set_from_float(r__1, z);
- if (mp_compare_mp_to_mp(&t, z) > 0) {
- mp_set_from_integer(x->sign, z);
- return;
- }
-
- /* 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) {
- mp_epowy(&t, &t);
- mp_add_integer(&t, -1, z);
- mp_add_integer(&t, 1, &t);
- mp_divide(z, &t, z);
- } else {
- mp_epowy(&t, &t);
- mp_add_integer(&t, 1, z);
- mp_add_integer(&t, -1, &t);
- mp_divide(&t, z, z);
- }
-
- /* Restore sign */
- z->sign = x->sign * z->sign;
+ mpc_tanh(z->num, x->num, MPC_RNDNN);
}
void
mp_asinh(const MPNumber *x, MPNumber *z)
{
- MPNumber t;
-
- /* sinh⁻¹(x) = ln(x + √(x² + 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);
+ mpc_asinh(z->num, x->num, MPC_RNDNN);
}
void
mp_acosh(const MPNumber *x, MPNumber *z)
{
- MPNumber t;
+ MPNumber t = mp_new();
/* Check x >= 1 */
mp_set_from_integer(1, &t);
- if (mp_is_less_than(x, &t)) {
+ if (mp_is_less_than(x, &t))
+ {
/* Translators: Error displayed when inverse hyperbolic cosine value is undefined */
mperr(_("Inverse hyperbolic cosine is undefined for values less than one"));
mp_set_from_integer(0, z);
return;
}
- /* cosh⁻¹(x) = ln(x + √(x² - 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);
+ mpc_acosh(z->num, x->num, MPC_RNDNN);
+ mp_clear(&t);
}
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)) {
- /* Translators: Error displayed when inverse hyperbolic tangent value is undefined */
- mperr(_("Inverse hyperbolic tangent is undefined for values outside [-1, 1]"));
+ MPNumber x_max = mp_new();
+ MPNumber x_min = mp_new();
+ mp_set_from_integer(1, &x_max);
+ mp_set_from_integer(-1, &x_min);
+
+ if (mp_compare(x, &x_max) >= 0 || mp_compare(x, &x_min) <= 0)
+ {
+ /* Translators: Error displayed when inverse hyperbolic tangent value is undefined */
+ mperr(_("Inverse hyperbolic tangent is undefined for values outside (-1, 1)"));
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);
+ mpc_atanh(z->num, x->num, MPC_RNDNN);
+ mp_clear(&x_max);
+ mp_clear(&x_min);
}
diff --git a/src/mp.c b/src/mp.c
index dcf4d92..a25280b 100644
--- a/src/mp.c
+++ b/src/mp.c
@@ -11,16 +11,11 @@
#include <stdlib.h>
#include <stdio.h>
+#include <stdint.h>
#include <math.h>
#include <errno.h>
#include "mp.h"
-#include "mp-private.h"
-
-static MPNumber eulers_number;
-static gboolean have_eulers_number = FALSE;
-
-// FIXME: Re-add overflow and underflow detection
char *mp_error = NULL;
@@ -58,818 +53,219 @@ void mp_clear_error()
}
-/* ROUTINE CALLED BY MP_DIVIDE AND MP_SQRT TO ENSURE THAT
- * RESULTS ARE REPRESENTED EXACTLY IN T-2 DIGITS IF THEY
- * CAN BE. X IS AN MP NUMBER, I AND J ARE INTEGERS.
- */
-static void
-mp_ext(int i, int j, MPNumber *x)
+MPNumber
+mp_new(void)
{
- int q, s;
+ MPNumber z;
+ mpc_init2(z.num, PRECISION);
+ return z;
+}
- if (mp_is_zero(x) || MP_T <= 2 || i == 0)
- return;
- /* COMPUTE MAXIMUM POSSIBLE ERROR IN THE LAST PLACE */
- q = (j + 1) / i + 1;
- s = MP_BASE * x->fraction[MP_T - 2] + x->fraction[MP_T - 1];
+MPNumber *
+mp_new_ptr(void)
+{
+ MPNumber *z = malloc(sizeof(MPNumber));
+ mpc_init2(z->num, PRECISION);
+ return z;
+}
- /* SET LAST TWO DIGITS TO ZERO */
- if (s <= q) {
- x->fraction[MP_T - 2] = 0;
- x->fraction[MP_T - 1] = 0;
- return;
- }
- if (s + q < MP_BASE * MP_BASE)
- return;
+void
+mp_clear(MPNumber *z)
+{
+ if (z != NULL)
+ mpc_clear(z->num);
+}
- /* ROUND UP HERE */
- 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);
+void
+mp_free(MPNumber *z)
+{
+ if (z != NULL)
+ {
+ mpc_clear(z->num);
+ free(z);
+ }
}
void
mp_get_eulers(MPNumber *z)
{
- if (!have_eulers_number) {
- MPNumber t;
- mp_set_from_integer(1, &t);
- mp_epowy(&t, &eulers_number);
- have_eulers_number = TRUE;
- }
- mp_set_from_mp(&eulers_number, z);
+ /* e^1, since mpfr doesn't have a function to return e */
+ mpfr_set_ui(mpc_realref(z->num), 1, MPFR_RNDN);
+ mpfr_exp(mpc_realref(z->num), mpc_realref(z->num), MPFR_RNDN);
+ mpfr_set_zero(mpc_imagref(z->num), 0);
}
void
mp_get_i(MPNumber *z)
{
- mp_set_from_integer(0, z);
- z->im_sign = 1;
- z->im_exponent = 1;
- z->im_fraction[0] = 1;
+ mpc_set_si_si(z->num, 0, 1, MPC_RNDNN);
}
void
mp_abs(const MPNumber *x, MPNumber *z)
{
- if (mp_is_complex(x)){
- MPNumber x_real, x_im;
-
- mp_real_component(x, &x_real);
- mp_imaginary_component(x, &x_im);
-
- mp_multiply(&x_real, &x_real, &x_real);
- mp_multiply(&x_im, &x_im, &x_im);
- mp_add(&x_real, &x_im, z);
- mp_root(z, 2, z);
- }
- else {
- mp_set_from_mp(x, z);
- if (z->sign < 0)
- z->sign = -z->sign;
- }
+ mpfr_set_zero(mpc_imagref(z->num), MPFR_RNDN);
+ mpc_abs(mpc_realref(z->num), x->num, MPC_RNDNN);
}
void
mp_arg(const MPNumber *x, MPAngleUnit unit, MPNumber *z)
{
- MPNumber x_real, x_im, pi;
-
- if (mp_is_zero(x)) {
+ if (mp_is_zero(x))
+ {
/* Translators: Error display when attempting to take argument of zero */
mperr(_("Argument not defined for zero"));
mp_set_from_integer(0, z);
return;
}
- mp_real_component(x, &x_real);
- mp_imaginary_component(x, &x_im);
- mp_get_pi(&pi);
-
- if (mp_is_zero(&x_im)) {
- if (mp_is_negative(&x_real))
- convert_from_radians(&pi, MP_RADIANS, z);
- else
- mp_set_from_integer(0, z);
- }
- else if (mp_is_zero(&x_real)) {
- mp_set_from_mp(&pi, z);
- if (mp_is_negative(&x_im))
- mp_divide_integer(z, -2, z);
- else
- mp_divide_integer(z, 2, z);
- }
- else if (mp_is_negative(&x_real)) {
- mp_divide(&x_im, &x_real, z);
- mp_atan(z, MP_RADIANS, z);
- if (mp_is_negative(&x_im))
- mp_subtract(z, &pi, z);
- else
- mp_add(z, &pi, z);
- }
- else {
- mp_divide(&x_im, &x_real, z);
- mp_atan(z, MP_RADIANS, z);
- }
-
+ mpfr_set_zero(mpc_imagref(z->num), MPFR_RNDN);
+ mpc_arg(mpc_realref(z->num), x->num, MPC_RNDNN);
convert_from_radians(z, unit, z);
+ // MPC returns -π for the argument of negative real numbers if
+ // their imaginary part is -0, we want +π for all real negative
+ // numbers
+ if (!mp_is_complex (x) && mp_is_negative (x))
+ mpfr_abs(mpc_realref(z->num), mpc_realref(z->num), MPFR_RNDN);
}
void
mp_conjugate(const MPNumber *x, MPNumber *z)
{
- mp_set_from_mp(x, z);
- z->im_sign = -z->im_sign;
+ mpc_conj(z->num, x->num, MPC_RNDNN);
}
void
mp_real_component(const MPNumber *x, MPNumber *z)
{
- mp_set_from_mp(x, z);
-
- /* Clear imaginary component */
- z->im_sign = 0;
- z->im_exponent = 0;
- memset(z->im_fraction, 0, sizeof(int) * MP_SIZE);
+ mpc_set_fr(z->num, mpc_realref(x->num), MPC_RNDNN);
}
void
mp_imaginary_component(const MPNumber *x, MPNumber *z)
{
- /* Copy imaginary component to real component */
- z->sign = x->im_sign;
- z->exponent = x->im_exponent;
- memcpy(z->fraction, x->im_fraction, sizeof(int) * MP_SIZE);
-
- /* Clear (old) imaginary component */
- z->im_sign = 0;
- z->im_exponent = 0;
- memset(z->im_fraction, 0, sizeof(int) * MP_SIZE);
+ mpc_set_fr(z->num, mpc_imagref(x->num), MPC_RNDNN);
}
-
-static void
-mp_add_real(const MPNumber *x, int y_sign, const MPNumber *y, MPNumber *z)
-{
- int sign_prod, i, c;
- int exp_diff, med;
- bool x_largest = false;
- const int *big_fraction, *small_fraction;
- MPNumber x_copy, y_copy;
-
- /* 0 + y = y */
- if (mp_is_zero(x)) {
- mp_set_from_mp(y, z);
- z->sign = y_sign;
- return;
- }
- /* x + 0 = x */
- else if (mp_is_zero(y)) {
- mp_set_from_mp(x, z);
- return;
- }
-
- sign_prod = y_sign * x->sign;
- exp_diff = x->exponent - y->exponent;
- med = abs(exp_diff);
- if (exp_diff < 0) {
- x_largest = false;
- } else if (exp_diff > 0) {
- x_largest = true;
- } else {
- /* EXPONENTS EQUAL SO COMPARE SIGNS, THEN FRACTIONS IF NEC. */
- if (sign_prod < 0) {
- /* Signs are not equal. find out which mantissa is larger. */
- int j;
- for (j = 0; j < MP_T; j++) {
- int i = x->fraction[j] - y->fraction[j];
- if (i == 0)
- continue;
- if (i < 0)
- x_largest = false;
- else if (i > 0)
- x_largest = true;
- break;
- }
-
- /* Both mantissas equal, so result is zero. */
- if (j >= MP_T) {
- mp_set_from_integer(0, z);
- return;
- }
- }
- }
-
- mp_set_from_mp(x, &x_copy);
- mp_set_from_mp(y, &y_copy);
- mp_set_from_integer(0, z);
-
- if (x_largest) {
- z->sign = x_copy.sign;
- z->exponent = x_copy.exponent;
- big_fraction = x_copy.fraction;
- small_fraction = y_copy.fraction;
- } else {
- z->sign = y_sign;
- z->exponent = y_copy.exponent;
- big_fraction = y_copy.fraction;
- small_fraction = x_copy.fraction;
- }
-
- /* CLEAR GUARD DIGITS TO RIGHT OF X DIGITS */
- for(i = 3; i >= med; i--)
- z->fraction[MP_T + i] = 0;
-
- if (sign_prod >= 0) {
- /* This is probably insufficient overflow detection, but it makes us
- * not crash at least.
- */
- if (MP_T + 3 < med) {
- mperr(_("Overflow: the result couldn't be calculated"));
- mp_set_from_integer(0, z);
- return;
- }
-
- /* HERE DO ADDITION, EXPONENT(Y) >= EXPONENT(X) */
- for (i = MP_T + 3; i >= MP_T; i--)
- z->fraction[i] = small_fraction[i - med];
-
- c = 0;
- for (; i >= med; i--) {
- c = big_fraction[i] + small_fraction[i - med] + c;
-
- if (c < MP_BASE) {
- /* NO CARRY GENERATED HERE */
- z->fraction[i] = c;
- c = 0;
- } else {
- /* CARRY GENERATED HERE */
- z->fraction[i] = c - MP_BASE;
- c = 1;
- }
- }
-
- for (; i >= 0; i--)
- {
- c = big_fraction[i] + c;
- if (c < MP_BASE) {
- z->fraction[i] = c;
- i--;
-
- /* NO CARRY POSSIBLE HERE */
- for (; i >= 0; i--)
- z->fraction[i] = big_fraction[i];
-
- c = 0;
- break;
- }
-
- z->fraction[i] = 0;
- c = 1;
- }
-
- /* MUST SHIFT RIGHT HERE AS CARRY OFF END */
- if (c != 0) {
- for (i = MP_T + 3; i > 0; i--)
- z->fraction[i] = z->fraction[i - 1];
- z->fraction[0] = 1;
- z->exponent++;
- }
- }
- else {
- c = 0;
- for (i = MP_T + med - 1; i >= MP_T; i--) {
- /* HERE DO SUBTRACTION, ABS(Y) > ABS(X) */
- z->fraction[i] = c - small_fraction[i - med];
- c = 0;
-
- /* BORROW GENERATED HERE */
- if (z->fraction[i] < 0) {
- c = -1;
- z->fraction[i] += MP_BASE;
- }
- }
-
- for(; i >= med; i--) {
- c = big_fraction[i] + c - small_fraction[i - med];
- if (c >= 0) {
- /* NO BORROW GENERATED HERE */
- z->fraction[i] = c;
- c = 0;
- } else {
- /* BORROW GENERATED HERE */
- z->fraction[i] = c + MP_BASE;
- c = -1;
- }
- }
-
- for (; i >= 0; i--) {
- c = big_fraction[i] + c;
-
- if (c >= 0) {
- z->fraction[i] = c;
- i--;
-
- /* NO CARRY POSSIBLE HERE */
- for (; i >= 0; i--)
- z->fraction[i] = big_fraction[i];
-
- break;
- }
-
- z->fraction[i] = c + MP_BASE;
- c = -1;
- }
- }
-
- mp_normalize(z);
-}
-
-
-static void
-mp_add_with_sign(const MPNumber *x, int y_sign, const MPNumber *y, MPNumber *z)
-{
- if (mp_is_complex(x) || mp_is_complex(y)) {
- MPNumber real_x, real_y, im_x, im_y, real_z, im_z;
-
- mp_real_component(x, &real_x);
- mp_imaginary_component(x, &im_x);
- mp_real_component(y, &real_y);
- mp_imaginary_component(y, &im_y);
-
- mp_add_real(&real_x, y_sign * y->sign, &real_y, &real_z);
- mp_add_real(&im_x, y_sign * y->im_sign, &im_y, &im_z);
-
- mp_set_from_complex(&real_z, &im_z, z);
- }
- else
- mp_add_real(x, y_sign * y->sign, y, z);
-}
-
-
void
mp_add(const MPNumber *x, const MPNumber *y, MPNumber *z)
{
- mp_add_with_sign(x, 1, y, z);
+ mpc_add(z->num, x->num, y->num, MPC_RNDNN);
}
void
mp_add_integer(const MPNumber *x, int64_t y, MPNumber *z)
{
- MPNumber t;
- mp_set_from_integer(y, &t);
- mp_add(x, &t, z);
+ mpc_add_si(z->num, x->num, y, MPC_RNDNN);
}
-
-void
-mp_add_fraction(const MPNumber *x, int64_t i, int64_t j, MPNumber *y)
-{
- MPNumber t;
- mp_set_from_fraction(i, j, &t);
- mp_add(x, &t, y);
-}
-
-
void
mp_subtract(const MPNumber *x, const MPNumber *y, MPNumber *z)
{
- mp_add_with_sign(x, -1, y, z);
+ mpc_sub(z->num, x->num, y->num, MPC_RNDNN);
}
-
void
mp_sgn(const MPNumber *x, MPNumber *z)
{
- if (mp_is_zero(x))
- mp_set_from_integer(0, z);
- else if (mp_is_negative(x))
- mp_set_from_integer(-1, z);
- else
- mp_set_from_integer(1, z);
+ mpc_set_si(z->num, mpfr_sgn(mpc_realref(x->num)), MPC_RNDNN);
}
void
mp_integer_component(const MPNumber *x, MPNumber *z)
{
- int i;
-
- /* Clear fraction */
- mp_set_from_mp(x, z);
- for (i = z->exponent; i < MP_SIZE; i++)
- z->fraction[i] = 0;
- z->im_sign = 0;
- z->im_exponent = 0;
- memset(z->im_fraction, 0, sizeof(int) * MP_SIZE);
+ mpfr_set_zero(mpc_imagref(z->num), MPFR_RNDN);
+ mpfr_trunc(mpc_realref(z->num), mpc_realref(x->num));
}
void
mp_fractional_component(const MPNumber *x, MPNumber *z)
{
- int i, shift;
-
- /* Fractional component of zero is 0 */
- if (mp_is_zero(x)) {
- mp_set_from_integer(0, z);
- return;
- }
-
- /* All fractional */
- if (x->exponent <= 0) {
- mp_set_from_mp(x, z);
- return;
- }
-
- /* 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;
-
- z->im_sign = 0;
- z->im_exponent = 0;
- memset(z->im_fraction, 0, sizeof(int) * MP_SIZE);
+ mpfr_set_zero(mpc_imagref(z->num), MPFR_RNDN);
+ mpfr_frac(mpc_realref(z->num), mpc_realref(x->num), MPFR_RNDN);
}
void
mp_fractional_part(const MPNumber *x, MPNumber *z)
{
- MPNumber f;
+ MPNumber f = mp_new();
mp_floor(x, &f);
mp_subtract(x, &f, z);
+ mp_clear(&f);
}
void
mp_floor(const MPNumber *x, MPNumber *z)
{
- int i;
- bool have_fraction = false, is_negative;
-
- /* Integer component of zero = 0 */
- if (mp_is_zero(x)) {
- mp_set_from_mp(x, z);
- return;
- }
-
- /* If all fractional then no integer component */
- if (x->exponent <= 0) {
- mp_set_from_integer(0, z);
- return;
- }
-
- is_negative = mp_is_negative(x);
-
- /* Clear fraction */
- mp_set_from_mp(x, z);
- for (i = z->exponent; i < MP_SIZE; i++) {
- if (z->fraction[i])
- have_fraction = true;
- z->fraction[i] = 0;
- }
- z->im_sign = 0;
- z->im_exponent = 0;
- memset(z->im_fraction, 0, sizeof(int) * MP_SIZE);
-
- if (have_fraction && is_negative)
- mp_add_integer(z, -1, z);
+ mpfr_set_zero(mpc_imagref(z->num), MPFR_RNDN);
+ mpfr_floor(mpc_realref(z->num), mpc_realref(x->num));
}
void
mp_ceiling(const MPNumber *x, MPNumber *z)
{
- MPNumber f;
-
- mp_floor(x, z);
- mp_fractional_component(x, &f);
- if (mp_is_zero(&f))
- return;
- mp_add_integer(z, 1, z);
+ mpfr_set_zero(mpc_imagref(z->num), MPFR_RNDN);
+ mpfr_ceil(mpc_realref(z->num), mpc_realref(x->num));
}
void
mp_round(const MPNumber *x, MPNumber *z)
{
- MPNumber f, one;
- bool do_floor;
-
- do_floor = !mp_is_negative(x);
-
- mp_fractional_component(x, &f);
- mp_multiply_integer(&f, 2, &f);
- mp_abs(&f, &f);
- mp_set_from_integer(1, &one);
- if (mp_is_greater_equal(&f, &one))
- do_floor = !do_floor;
-
- if (do_floor)
- mp_floor(x, z);
- else
- mp_ceiling(x, z);
+ mpfr_set_zero(mpc_imagref(z->num), MPFR_RNDN);
+ mpfr_round(mpc_realref(z->num), mpc_realref(x->num));
}
int
-mp_compare_mp_to_mp(const MPNumber *x, const MPNumber *y)
+mp_compare(const MPNumber *x, const MPNumber *y)
{
- int i;
-
- if (x->sign != y->sign) {
- if (x->sign > y->sign)
- return 1;
- else
- return -1;
- }
-
- /* x = y = 0 */
- if (mp_is_zero(x))
- return 0;
-
- /* See if numbers are of different magnitude */
- if (x->exponent != y->exponent) {
- if (x->exponent > y->exponent)
- return x->sign;
- else
- return -x->sign;
- }
-
- /* Compare fractions */
- for (i = 0; i < MP_SIZE; i++) {
- if (x->fraction[i] == y->fraction[i])
- continue;
-
- if (x->fraction[i] > y->fraction[i])
- return x->sign;
- else
- return -x->sign;
- }
-
- /* x = y */
- return 0;
+ return mpfr_cmp(mpc_realref(x->num), mpc_realref(y->num));
}
-
void
mp_divide(const MPNumber *x, const MPNumber *y, MPNumber *z)
{
- int i, ie;
- MPNumber t;
-
- /* x/0 */
- if (mp_is_zero(y)) {
- /* Translators: Error displayed attempted to divide by zero */
- mperr(_("Division by zero is undefined"));
- mp_set_from_integer(0, z);
- return;
- }
-
- /* 0/y = 0 */
- if (mp_is_zero(x)) {
- mp_set_from_integer(0, z);
- return;
- }
-
- /* z = x × y⁻¹ */
- /* FIXME: Set exponent to zero to avoid overflow in mp_multiply??? */
- mp_reciprocal(y, &t);
- ie = t.exponent;
- t.exponent = 0;
- i = t.fraction[0];
- mp_multiply(x, &t, z);
- mp_ext(i, z->fraction[0], z);
- z->exponent += ie;
-}
-
-
-static void
-mp_divide_integer_real(const MPNumber *x, int64_t y, MPNumber *z)
-{
- int c, i, k, b2, c2, j1, j2;
- MPNumber x_copy;
-
- /* x/0 */
- if (y == 0) {
+ if (mp_is_zero(y))
+ {
/* Translators: Error displayed attempted to divide by zero */
mperr(_("Division by zero is undefined"));
mp_set_from_integer(0, z);
return;
}
-
- /* 0/y = 0 */
- if (mp_is_zero(x)) {
- mp_set_from_integer(0, z);
- return;
- }
-
- /* Division by -1 or 1 just changes sign */
- if (y == 1 || y == -1) {
- if (y < 0)
- mp_invert_sign(x, z);
- else
- mp_set_from_mp(x, z);
- return;
- }
-
- /* Copy x as z may also refer to x */
- mp_set_from_mp(x, &x_copy);
- mp_set_from_integer(0, z);
-
- if (y < 0) {
- y = -y;
- z->sign = -x_copy.sign;
- }
- else
- z->sign = x_copy.sign;
- z->exponent = x_copy.exponent;
-
- c = 0;
- i = 0;
-
- /* IF y*B NOT REPRESENTABLE AS AN INTEGER HAVE TO SIMULATE
- * LONG DIVISION. ASSUME AT LEAST 16-BIT WORD.
- */
-
- /* Computing MAX */
- b2 = max(MP_BASE << 3, 32767 / MP_BASE);
- if (y < b2) {
- int kh, r1;
-
- /* LOOK FOR FIRST NONZERO DIGIT IN QUOTIENT */
- do {
- c = MP_BASE * c;
- if (i < MP_T)
- c += x_copy.fraction[i];
- i++;
- r1 = c / y;
- if (r1 < 0)
- goto L210;
- } while (r1 == 0);
-
- /* ADJUST EXPONENT AND GET T+4 DIGITS IN QUOTIENT */
- z->exponent += 1 - i;
- z->fraction[0] = r1;
- c = MP_BASE * (c - y * r1);
- kh = 1;
- if (i < MP_T) {
- kh = MP_T + 1 - i;
- for (k = 1; k < kh; k++) {
- c += x_copy.fraction[i];
- z->fraction[k] = c / y;
- c = MP_BASE * (c - y * z->fraction[k]);
- i++;
- }
- if (c < 0)
- goto L210;
- }
-
- for (k = kh; k < MP_T + 4; k++) {
- z->fraction[k] = c / y;
- c = MP_BASE * (c - y * z->fraction[k]);
- }
- if (c < 0)
- goto L210;
-
- mp_normalize(z);
- return;
- }
-
- /* HERE NEED SIMULATED DOUBLE-PRECISION DIVISION */
- j1 = y / MP_BASE;
- j2 = y - j1 * MP_BASE;
-
- /* LOOK FOR FIRST NONZERO DIGIT */
- c2 = 0;
- do {
- c = MP_BASE * c + c2;
- c2 = i < MP_T ? x_copy.fraction[i] : 0;
- i++;
- } while (c < j1 || (c == j1 && c2 < j2));
-
- /* COMPUTE T+4 QUOTIENT DIGITS */
- z->exponent += 1 - i;
- i--;
-
- /* MAIN LOOP FOR LARGE ABS(y) CASE */
- for (k = 1; k <= MP_T + 4; k++) {
- int ir, iq, iqj;
-
- /* GET APPROXIMATE QUOTIENT FIRST */
- ir = c / (j1 + 1);
-
- /* NOW REDUCE SO OVERFLOW DOES NOT OCCUR */
- iq = c - ir * j1;
- if (iq >= b2) {
- /* HERE IQ*B WOULD POSSIBLY OVERFLOW SO INCREASE IR */
- ++ir;
- iq -= j1;
- }
-
- iq = iq * MP_BASE - ir * j2;
- if (iq < 0) {
- /* HERE IQ NEGATIVE SO IR WAS TOO LARGE */
- ir--;
- iq += y;
- }
-
- if (i < MP_T)
- iq += x_copy.fraction[i];
- i++;
- iqj = iq / y;
-
- /* R(K) = QUOTIENT, C = REMAINDER */
- z->fraction[k - 1] = iqj + ir;
- c = iq - y * iqj;
-
- if (c < 0)
- goto L210;
- }
-
- mp_normalize(z);
-
-L210:
- /* CARRY NEGATIVE SO OVERFLOW MUST HAVE OCCURRED */
- mperr("*** INTEGER OVERFLOW IN MP_DIVIDE_INTEGER, B TOO LARGE ***");
- mp_set_from_integer(0, z);
+ mpc_div(z->num, x->num, y->num, MPC_RNDNN);
}
-
void
mp_divide_integer(const MPNumber *x, int64_t y, MPNumber *z)
{
- if (mp_is_complex(x)) {
- MPNumber re_z, im_z;
-
- mp_real_component(x, &re_z);
- mp_imaginary_component(x, &im_z);
- mp_divide_integer_real(&re_z, y, &re_z);
- mp_divide_integer_real(&im_z, y, &im_z);
- mp_set_from_complex(&re_z, &im_z, z);
- }
- else
- mp_divide_integer_real(x, y, z);
-}
+ MPNumber t1 = mp_new();
+ mp_set_from_integer(y, &t1);
+ mp_divide(x, &t1, z);
+ mp_clear(&t1);
+}
bool
mp_is_integer(const MPNumber *x)
{
- MPNumber t1, t2, t3;
-
if (mp_is_complex(x))
return false;
- /* 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_floor() routine in mp.c, when the exponent is less than 1.
- */
- mp_set_from_integer(10000, &t3);
- mp_multiply(x, &t3, &t1);
- mp_divide(&t1, &t3, &t1);
- mp_floor(&t1, &t2);
- return mp_is_equal(&t1, &t2);
-
- /* Correct way to check for integer */
- /*int i;
-
- // Zero is an integer
- if (mp_is_zero(x))
- return true;
-
- // Fractional
- if (x->exponent <= 0)
- return false;
-
- // Look for fractional components
- for (i = x->exponent; i < MP_SIZE; i++) {
- if (x->fraction[i] != 0)
- return false;
- }
-
- return true;*/
+ return mpfr_integer_p(mpc_realref(x->num)) != 0;
}
@@ -879,7 +275,7 @@ mp_is_positive_integer(const MPNumber *x)
if (mp_is_complex(x))
return false;
else
- return x->sign >= 0 && mp_is_integer(x);
+ return mpfr_sgn(mpc_realref(x->num)) >= 0 && mp_is_integer(x);
}
@@ -889,389 +285,77 @@ mp_is_natural(const MPNumber *x)
if (mp_is_complex(x))
return false;
else
- return x->sign > 0 && mp_is_integer(x);
+ return mpfr_sgn(mpc_realref(x->num)) > 0 && mp_is_integer(x);
}
bool
mp_is_complex(const MPNumber *x)
{
- return x->im_sign != 0;
+ return !mpfr_zero_p(mpc_imagref(x->num));
}
bool
mp_is_equal(const MPNumber *x, const MPNumber *y)
{
- return mp_compare_mp_to_mp(x, y) == 0;
-}
-
-
-/* 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.
- */
-static void
-mp_exp(const MPNumber *x, MPNumber *z)
-{
- int i, q;
- float rlb;
- MPNumber t1, t2;
-
- /* e^0 = 1 */
- if (mp_is_zero(x)) {
- mp_set_from_integer(1, z);
- return;
- }
-
- /* Only defined for |x| < 1 */
- if (x->exponent > 0) {
- mperr("*** ABS(X) NOT LESS THAN 1 IN CALL TO MP_EXP ***");
- 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) {
- int ib, ic;
-
- ib = MP_BASE << 2;
- ic = 1;
- for (i = 1; i <= q; ++i) {
- ic *= 2;
- if (ic < ib && ic != MP_BASE && i < q)
- continue;
- mp_divide_integer(&t1, ic, &t1);
- ic = 1;
- }
- }
-
- if (mp_is_zero(&t1)) {
- mp_set_from_integer(0, z);
- return;
- }
-
- /* Sum series, reducing t where possible */
- mp_set_from_mp(&t1, z);
- mp_set_from_mp(&t1, &t2);
- for (i = 2; MP_T + t2.exponent - z->exponent > 0; i++) {
- mp_multiply(&t1, &t2, &t2);
- mp_divide_integer(&t2, i, &t2);
- mp_add(&t2, z, z);
- if (mp_is_zero(&t2))
- 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);
- }
-
- mp_add_integer(z, 1, z);
-}
-
-
-static void
-mp_epowy_real(const MPNumber *x, MPNumber *z)
-{
- float r__1;
- int i, ix, xs, tss;
- float rx, rz;
- MPNumber t1, t2;
-
- /* e^0 = 1 */
- if (mp_is_zero(x)) {
- mp_set_from_integer(1, z);
- return;
- }
-
- /* If |x| < 1 use mp_exp */
- if (x->exponent <= 0) {
- mp_exp(x, z);
- return;
- }
-
- /* NOW SAFE TO CONVERT X TO REAL */
- rx = mp_cast_to_float(x);
-
- /* SAVE SIGN AND WORK WITH ABS(X) */
- xs = x->sign;
- mp_abs(x, &t2);
-
- /* GET FRACTIONAL AND INTEGER PARTS OF ABS(X) */
- ix = mp_cast_to_int(&t2);
- mp_fractional_component(&t2, &t2);
-
- /* ATTACH SIGN TO FRACTIONAL PART AND COMPUTE EXP OF IT */
- t2.sign *= xs;
- mp_exp(&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;
- else
- tss = MP_T + 2;
-
- /* LOOP FOR E COMPUTATION. DECREASE T IF POSSIBLE. */
- /* Computing MIN */
- mp_set_from_integer(xs, &t1);
-
- t2.sign = 0;
- for (i = 2 ; ; i++) {
- if (min(tss, tss + 2 + t1.exponent) <= 2)
- break;
-
- mp_divide_integer(&t1, i * xs, &t1);
- mp_add(&t2, &t1, &t2);
- if (mp_is_zero(&t1))
- break;
- }
-
- /* RAISE E OR 1/E TO POWER IX */
- if (xs > 0)
- mp_add_integer(&t2, 2, &t2);
- mp_xpowy_integer(&t2, ix, &t2);
-
- /* MULTIPLY EXPS OF INTEGER AND FRACTIONAL PARTS */
- mp_multiply(z, &t2, z);
-
- /* CHECK THAT RELATIVE ERROR LESS THAN 0.01 UNLESS ABS(X) LARGE
- * (WHEN EXP MIGHT OVERFLOW OR UNDERFLOW)
- */
- if (fabs(rx) > 10.0f)
- return;
-
- rz = mp_cast_to_float(z);
- r__1 = rz - exp(rx);
- if (fabs(r__1) < rz * 0.01f)
- return;
-
- /* THE FOLLOWING MESSAGE MAY INDICATE THAT
- * B**(T-1) IS TOO SMALL, OR THAT M IS TOO SMALL SO THE
- * RESULT UNDERFLOWED.
- */
- mperr("*** ERROR OCCURRED IN MP_EPOWY, RESULT INCORRECT ***");
+ int res = mpc_cmp(x->num, y->num);
+ return MPC_INEX_RE(res) == 0 && MPC_INEX_IM(res) == 0;
}
void
mp_epowy(const MPNumber *x, MPNumber *z)
{
- /* e^0 = 1 */
- if (mp_is_zero(x)) {
- mp_set_from_integer(1, z);
- return;
- }
-
- if (mp_is_complex(x)) {
- MPNumber x_real, r, theta;
-
- mp_real_component(x, &x_real);
- mp_imaginary_component(x, &theta);
-
- mp_epowy_real(&x_real, &r);
- mp_set_from_polar(&r, MP_RADIANS, &theta, z);
- }
- else
- mp_epowy_real(x, z);
-}
-
-
-/* RETURNS K = K/GCD AND L = L/GCD, WHERE GCD IS THE
- * GREATEST COMMON DIVISOR OF K AND L.
- * SAVE INPUT PARAMETERS IN LOCAL VARIABLES
- */
-void
-mp_gcd(int64_t *k, int64_t *l)
-{
- int64_t i, j;
-
- i = abs(*k);
- j = abs(*l);
- if (j == 0) {
- /* IF J = 0 RETURN (1, 0) UNLESS I = 0, THEN (0, 0) */
- *k = 1;
- *l = 0;
- if (i == 0)
- *k = 0;
- return;
- }
-
- /* EUCLIDEAN ALGORITHM LOOP */
- do {
- i %= j;
- if (i == 0) {
- *k = *k / j;
- *l = *l / j;
- return;
- }
- j %= i;
- } while (j != 0);
-
- /* HERE J IS THE GCD OF K AND L */
- *k = *k / i;
- *l = *l / i;
+ mpc_exp(z->num, x->num, MPC_RNDNN);
}
-
bool
-mp_is_zero(const MPNumber *x)
+mp_is_zero (const MPNumber *x)
{
- return x->sign == 0 && x->im_sign == 0;
+ int res = mpc_cmp_si_si(x->num, 0, 0);
+ return MPC_INEX_RE(res) == 0 && MPC_INEX_IM(res) == 0;
}
-
bool
mp_is_negative(const MPNumber *x)
{
- return x->sign < 0;
+ return mpfr_sgn(mpc_realref(x->num)) < 0;
}
bool
mp_is_greater_equal(const MPNumber *x, const MPNumber *y)
{
- return mp_compare_mp_to_mp(x, y) >= 0;
+ return mp_compare(x, y) >= 0;
}
bool
mp_is_greater_than(const MPNumber *x, const MPNumber *y)
{
- return mp_compare_mp_to_mp(x, y) > 0;
+ return mp_compare(x, y) > 0;
}
bool
mp_is_less_equal(const MPNumber *x, const MPNumber *y)
{
- return mp_compare_mp_to_mp(x, y) <= 0;
-}
-
-
-/* RETURNS MP Y = LN(1+X) IF X IS AN MP NUMBER SATISFYING THE
- * CONDITION ABS(X) < 1/B, ERROR OTHERWISE.
- * USES NEWTONS METHOD TO SOLVE THE EQUATION
- * EXP1(-Y) = X, THEN REVERSES SIGN OF Y.
- */
-static void
-mp_lns(const MPNumber *x, MPNumber *z)
-{
- int t, it0;
- MPNumber t1, t2, t3;
-
- /* ln(1+0) = 0 */
- if (mp_is_zero(x)) {
- mp_set_from_integer(0, z);
- return;
- }
-
- /* Get starting approximation -ln(1+x) ~= -x + x^2/2 - x^3/3 + x^4/4 */
- mp_set_from_mp(x, &t2);
- mp_divide_integer(x, 4, &t1);
- mp_add_fraction(&t1, -1, 3, &t1);
- mp_multiply(x, &t1, &t1);
- mp_add_fraction(&t1, 1, 2, &t1);
- mp_multiply(x, &t1, &t1);
- mp_add_integer(&t1, -1, &t1);
- mp_multiply(x, &t1, z);
-
- /* Solve using Newtons method */
- it0 = t = 5;
- while(1)
- {
- int ts2, ts3;
-
- /* t3 = (e^t3 - 1) */
- /* z = z - (t2 + t3 + (t2 * 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(z, &t3, z);
- if (t >= MP_T)
- break;
-
- /* FOLLOWING LOOP COMPUTES NEXT VALUE OF T TO USE.
- * BECAUSE NEWTONS METHOD HAS 2ND ORDER CONVERGENCE,
- * WE CAN ALMOST DOUBLE T EACH TIME.
- */
- ts3 = t;
- t = MP_T;
- do {
- ts2 = t;
- t = (t + it0) / 2;
- } while (t > ts3);
- t = ts2;
- }
-
- /* CHECK THAT NEWTON ITERATION WAS CONVERGING AS EXPECTED */
- if (t3.sign != 0 && t3.exponent << 1 > it0 - MP_T) {
- mperr("*** ERROR OCCURRED IN MP_LNS, NEWTON ITERATION NOT CONVERGING PROPERLY ***");
- }
-
- z->sign = -z->sign;
+ return mp_compare(x, y) <= 0;
}
-
-static void
-mp_ln_real(const MPNumber *x, MPNumber *z)
+bool
+mp_is_less_than(const MPNumber *x, const MPNumber *y)
{
- int e, k;
- double rx, rlx;
- MPNumber t1, t2;
-
- /* LOOP TO GET APPROXIMATE LN(X) USING SINGLE-PRECISION */
- mp_set_from_mp(x, &t1);
- mp_set_from_integer(0, z);
- for(k = 0; k < 10; k++)
- {
- /* COMPUTE FINAL CORRECTION ACCURATELY USING MP_LNS */
- mp_add_integer(&t1, -1, &t2);
- if (mp_is_zero(&t2) || t2.exponent + 1 <= 0) {
- mp_lns(&t2, &t2);
- mp_add(z, &t2, z);
- return;
- }
-
- /* REMOVE EXPONENT TO AVOID FLOATING-POINT OVERFLOW */
- e = t1.exponent;
- t1.exponent = 0;
- rx = mp_cast_to_double(&t1);
- t1.exponent = e;
- rlx = log(rx) + e * log(MP_BASE);
- mp_set_from_double(-(double)rlx, &t2);
-
- /* UPDATE Z AND COMPUTE ACCURATE EXP OF APPROXIMATE LOG */
- mp_subtract(z, &t2, z);
- mp_epowy(&t2, &t2);
-
- /* COMPUTE RESIDUAL WHOSE LOG IS STILL TO BE FOUND */
- mp_multiply(&t1, &t2, &t1);
- }
-
- mperr("*** ERROR IN MP_LN, ITERATION NOT CONVERGING ***");
+ return mp_compare(x, y) < 0;
}
-
void
mp_ln(const MPNumber *x, MPNumber *z)
{
/* ln(0) undefined */
- if (mp_is_zero(x)) {
+ if (mp_is_zero(x))
+ {
/* Translators: Error displayed when attempting to take logarithm of zero */
mperr(_("Logarithm of zero is undefined"));
mp_set_from_integer(0, z);
@@ -1287,28 +371,20 @@ mp_ln(const MPNumber *x, MPNumber *z)
return;
}*/
- if (mp_is_complex(x) || mp_is_negative(x)) {
- MPNumber r, theta, z_real;
-
- /* ln(re^iθ) = e^(ln(r)+iθ) */
- mp_abs(x, &r);
- mp_arg(x, MP_RADIANS, &theta);
-
- mp_ln_real(&r, &z_real);
- mp_set_from_complex(&z_real, &theta, z);
- }
- else
- mp_ln_real(x, z);
+ mpc_log(z->num, x->num, MPC_RNDNN);
+ // MPC returns -π for the imaginary part of the log of
+ // negative real numbers if their imaginary part is -0, we want +π
+ if (!mp_is_complex (x) && mp_is_negative (x))
+ mpfr_abs(mpc_imagref(z->num), mpc_imagref(z->num), MPFR_RNDN);
}
void
mp_logarithm(int64_t n, const MPNumber *x, MPNumber *z)
{
- MPNumber t1, t2;
-
/* log(0) undefined */
- if (mp_is_zero(x)) {
+ if (mp_is_zero(x))
+ {
/* Translators: Error displayed when attempting to take logarithm of zero */
mperr(_("Logarithm of zero is undefined"));
mp_set_from_integer(0, z);
@@ -1316,781 +392,339 @@ mp_logarithm(int64_t n, const MPNumber *x, MPNumber *z)
}
/* logn(x) = ln(x) / ln(n) */
+ MPNumber t1 = mp_new();
mp_set_from_integer(n, &t1);
mp_ln(&t1, &t1);
- mp_ln(x, &t2);
- mp_divide(&t2, &t1, z);
-}
-
-
-bool
-mp_is_less_than(const MPNumber *x, const MPNumber *y)
-{
- return mp_compare_mp_to_mp(x, y) < 0;
+ mp_ln(x, z);
+ mp_divide(z, &t1, z);
+ mp_clear(&t1);
}
-
-static void
-mp_multiply_real(const MPNumber *x, const MPNumber *y, MPNumber *z)
-{
- int c, i, xi;
- MPNumber r;
-
- /* x*0 = 0*y = 0 */
- if (x->sign == 0 || y->sign == 0) {
- mp_set_from_integer(0, z);
- return;
- }
-
- z->sign = x->sign * y->sign;
- z->exponent = x->exponent + y->exponent;
- memset(&r, 0, sizeof(MPNumber));
-
- /* PERFORM MULTIPLICATION */
- c = 8;
- for (i = 0; i < MP_T; i++) {
- int j;
-
- xi = x->fraction[i];
-
- /* FOR SPEED, PUT THE NUMBER WITH MANY ZEROS FIRST */
- if (xi == 0)
- continue;
-
- /* Computing MIN */
- for (j = 0; j < min(MP_T, MP_T + 3 - i); 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_BASE) {
- mperr("*** ILLEGAL BASE B DIGIT IN CALL TO MP_MULTIPLY, POSSIBLE OVERWRITING PROBLEM ***");
- mp_set_from_integer(0, z);
- return;
- }
-
- /* PROPAGATE CARRIES AT END AND EVERY EIGHTH TIME,
- * FASTER THAN DOING IT EVERY TIME.
- */
- for (j = MP_T + 3; j >= 0; j--) {
- int ri = r.fraction[j] + c;
- if (ri < 0) {
- mperr("*** INTEGER OVERFLOW IN MP_MULTIPLY, B TOO LARGE ***");
- mp_set_from_integer(0, z);
- return;
- }
- 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 ***");
- mp_set_from_integer(0, z);
- return;
- }
- c = 8;
- }
-
- if (c != 8) {
- if (xi < 0 || xi >= MP_BASE) {
- mperr("*** ILLEGAL BASE B DIGIT IN CALL TO MP_MULTIPLY, POSSIBLE OVERWRITING PROBLEM ***");
- mp_set_from_integer(0, z);
- return;
- }
-
- c = 0;
- for (i = MP_T + 3; i >= 0; i--) {
- int ri = r.fraction[i] + c;
- if (ri < 0) {
- mperr("*** INTEGER OVERFLOW IN MP_MULTIPLY, B TOO LARGE ***");
- mp_set_from_integer(0, z);
- return;
- }
- c = ri / MP_BASE;
- r.fraction[i] = ri - MP_BASE * c;
- }
-
- if (c != 0) {
- mperr("*** ILLEGAL BASE B DIGIT IN CALL TO MP_MULTIPLY, POSSIBLE OVERWRITING PROBLEM ***");
- mp_set_from_integer(0, z);
- return;
- }
- }
-
- /* Clear complex part */
- z->im_sign = 0;
- z->im_exponent = 0;
- memset(z->im_fraction, 0, sizeof(int) * MP_SIZE);
-
- /* NORMALIZE AND ROUND RESULT */
- // FIXME: Use stack variable because of mp_normalize brokeness
- for (i = 0; i < MP_SIZE; i++)
- z->fraction[i] = r.fraction[i];
- mp_normalize(z);
-}
-
-
void
mp_multiply(const MPNumber *x, const MPNumber *y, MPNumber *z)
{
- /* x*0 = 0*y = 0 */
- if (mp_is_zero(x) || mp_is_zero(y)) {
- mp_set_from_integer(0, z);
- return;
- }
-
- /* (a+bi)(c+di) = (ac-bd)+(ad+bc)i */
- if (mp_is_complex(x) || mp_is_complex(y)) {
- MPNumber real_x, real_y, im_x, im_y, t1, t2, real_z, im_z;
-
- mp_real_component(x, &real_x);
- mp_imaginary_component(x, &im_x);
- mp_real_component(y, &real_y);
- mp_imaginary_component(y, &im_y);
-
- mp_multiply_real(&real_x, &real_y, &t1);
- mp_multiply_real(&im_x, &im_y, &t2);
- mp_subtract(&t1, &t2, &real_z);
-
- mp_multiply_real(&real_x, &im_y, &t1);
- mp_multiply_real(&im_x, &real_y, &t2);
- mp_add(&t1, &t2, &im_z);
-
- mp_set_from_complex(&real_z, &im_z, z);
- }
- else {
- mp_multiply_real(x, y, z);
- }
-}
-
-
-static void
-mp_multiply_integer_real(const MPNumber *x, int64_t y, MPNumber *z)
-{
- int c, i;
- MPNumber x_copy;
-
- /* x*0 = 0*y = 0 */
- if (mp_is_zero(x) || y == 0) {
- mp_set_from_integer(0, z);
- return;
- }
-
- /* x*1 = x, x*-1 = -x */
- // FIXME: Why is this not working? mp_ext is using this function to do a normalization
- /*if (y == 1 || y == -1) {
- if (y < 0)
- mp_invert_sign(x, z);
- else
- mp_set_from_mp(x, z);
- return;
- }*/
-
- /* Copy x as z may also refer to x */
- mp_set_from_mp(x, &x_copy);
- mp_set_from_integer(0, z);
-
- if (y < 0) {
- y = -y;
- z->sign = -x_copy.sign;
- }
- else
- z->sign = x_copy.sign;
- z->exponent = x_copy.exponent + 4;
-
- /* FORM PRODUCT IN ACCUMULATOR */
- c = 0;
-
- /* IF y*B NOT REPRESENTABLE AS AN INTEGER WE HAVE TO SIMULATE
- * DOUBLE-PRECISION MULTIPLICATION.
- */
-
- /* Computing MAX */
- if (y >= max(MP_BASE << 3, 32767 / MP_BASE)) {
- int64_t j1, j2;
-
- /* HERE J IS TOO LARGE FOR SINGLE-PRECISION MULTIPLICATION */
- j1 = y / MP_BASE;
- j2 = y - j1 * MP_BASE;
-
- /* FORM PRODUCT */
- for (i = MP_T + 3; i >= 0; i--) {
- int64_t c1, c2, is, ix, t;
-
- c1 = c / MP_BASE;
- c2 = c - MP_BASE * c1;
- ix = 0;
- if (i > 3)
- ix = x_copy.fraction[i - 4];
-
- t = j2 * ix + c2;
- is = t / MP_BASE;
- c = j1 * ix + c1 + is;
- z->fraction[i] = t - MP_BASE * is;
- }
- }
- else
- {
- int64_t ri = 0;
-
- for (i = MP_T + 3; i >= 4; i--) {
- ri = y * x_copy.fraction[i - 4] + c;
- c = ri / MP_BASE;
- z->fraction[i] = ri - MP_BASE * c;
- }
-
- /* CHECK FOR INTEGER OVERFLOW */
- if (ri < 0) {
- mperr("*** INTEGER OVERFLOW IN mp_multiply_integer, B TOO LARGE ***");
- mp_set_from_integer(0, z);
- return;
- }
-
- /* HAVE TO TREAT FIRST FOUR WORDS OF R SEPARATELY */
- for (i = 3; i >= 0; i--) {
- int t;
-
- t = c;
- c = t / MP_BASE;
- z->fraction[i] = t - MP_BASE * c;
- }
- }
-
- /* HAVE TO SHIFT RIGHT HERE AS CARRY OFF END */
- while (c != 0) {
- int64_t t;
-
- for (i = MP_T + 3; i >= 1; i--)
- z->fraction[i] = z->fraction[i - 1];
- t = c;
- c = t / MP_BASE;
- z->fraction[0] = t - MP_BASE * c;
- z->exponent++;
- }
-
- if (c < 0) {
- mperr("*** INTEGER OVERFLOW IN mp_multiply_integer, B TOO LARGE ***");
- mp_set_from_integer(0, z);
- return;
- }
-
- z->im_sign = 0;
- z->im_exponent = 0;
- memset(z->im_fraction, 0, sizeof(int) * MP_SIZE);
- mp_normalize(z);
+ mpc_mul(z->num, x->num, y->num, MPC_RNDNN);
}
void
mp_multiply_integer(const MPNumber *x, int64_t y, MPNumber *z)
{
- if (mp_is_complex(x)) {
- MPNumber re_z, im_z;
- mp_real_component(x, &re_z);
- mp_imaginary_component(x, &im_z);
- mp_multiply_integer_real(&re_z, y, &re_z);
- mp_multiply_integer_real(&im_z, y, &im_z);
- mp_set_from_complex(&re_z, &im_z, z);
- }
- else
- mp_multiply_integer_real(x, y, z);
-}
-
-
-void
-mp_multiply_fraction(const MPNumber *x, int64_t numerator, int64_t denominator, MPNumber *z)
-{
- if (denominator == 0) {
- mperr(_("Division by zero is undefined"));
- mp_set_from_integer(0, z);
- return;
- }
-
- if (numerator == 0) {
- mp_set_from_integer(0, z);
- return;
- }
-
- /* Reduce to lowest terms */
- mp_gcd(&numerator, &denominator);
- mp_divide_integer(x, denominator, z);
- mp_multiply_integer(z, numerator, z);
+ mpc_mul_si(z->num, x->num, (long) y, MPC_RNDNN);
}
void
mp_invert_sign(const MPNumber *x, MPNumber *z)
{
- mp_set_from_mp(x, z);
- z->sign = -z->sign;
- z->im_sign = -z->im_sign;
-}
-
-
-// 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)
-void
-mp_normalize(MPNumber *x)
-{
- int start_index;
-
- /* Find first non-zero digit */
- for (start_index = 0; start_index < MP_SIZE && x->fraction[start_index] == 0; start_index++);
-
- /* Mark as zero */
- if (start_index >= MP_SIZE) {
- x->sign = 0;
- x->exponent = 0;
- return;
- }
-
- /* Shift left so first digit is non-zero */
- if (start_index > 0) {
- int i;
-
- x->exponent -= start_index;
- for (i = 0; (i + start_index) < MP_SIZE; i++)
- x->fraction[i] = x->fraction[i + start_index];
- for (; i < MP_SIZE; i++)
- x->fraction[i] = 0;
- }
-}
-
-
-static void
-mp_pwr(const MPNumber *x, const MPNumber *y, MPNumber *z)
-{
- MPNumber t;
-
- /* (-x)^y imaginary */
- /* FIXME: Make complex numbers optional */
- /*if (x->sign < 0) {
- mperr(_("The power of negative numbers is only defined for integer exponents"));
- mp_set_from_integer(0, z);
- return;
- }*/
-
- /* 0^y = 0, 0^-y undefined */
- if (mp_is_zero(x)) {
- mp_set_from_integer(0, z);
- if (y->sign < 0)
- mperr(_("The power of zero is undefined for a negative exponent"));
- return;
- }
-
- /* x^0 = 1 */
- if (mp_is_zero(y)) {
- mp_set_from_integer(1, z);
- return;
- }
-
- mp_ln(x, &t);
- mp_multiply(y, &t, z);
- mp_epowy(z, z);
-}
-
-
-static void
-mp_reciprocal_real(const MPNumber *x, MPNumber *z)
-{
- MPNumber t1, t2;
- int it0, t;
-
- /* 1/0 invalid */
- if (mp_is_zero(x)) {
- mperr(_("Reciprocal of zero is undefined"));
- mp_set_from_integer(0, z);
- return;
- }
-
- /* Start by approximating value using floating point */
- mp_set_from_mp(x, &t1);
- t1.exponent = 0;
- mp_set_from_float(1.0f / mp_cast_to_float(&t1), &t1);
- t1.exponent -= x->exponent;
-
- it0 = t = 3;
- while(1) {
- int ts2, ts3;
-
- /* t1 = t1 - (t1 * ((x * t1) - 1)) (2*t1 - t1^2*x) */
- mp_multiply(x, &t1, &t2);
- mp_add_integer(&t2, -1, &t2);
- mp_multiply(&t1, &t2, &t2);
- mp_subtract(&t1, &t2, &t1);
- if (t >= MP_T)
- break;
-
- /* FOLLOWING LOOP ALMOST DOUBLES T (POSSIBLE
- * BECAUSE NEWTONS METHOD HAS 2ND ORDER CONVERGENCE).
- */
- ts3 = t;
- t = MP_T;
- do {
- ts2 = t;
- t = (t + it0) / 2;
- } while (t > ts3);
- t = min(ts2, MP_T);
- }
-
- /* RETURN IF NEWTON ITERATION WAS CONVERGING */
- 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.
- */
- mperr("*** ERROR OCCURRED IN MP_RECIPROCAL, NEWTON ITERATION NOT CONVERGING PROPERLY ***");
- }
-
- mp_set_from_mp(&t1, z);
+ mpc_neg(z->num, x->num, MPC_RNDNN);
}
void
mp_reciprocal(const MPNumber *x, MPNumber *z)
{
- if (mp_is_complex(x)) {
- MPNumber t1, t2;
- MPNumber real_x, im_x;
-
- mp_real_component(x, &real_x);
- mp_imaginary_component(x, &im_x);
-
- /* 1/(a+bi) = (a-bi)/(a+bi)(a-bi) = (a-bi)/(a²+b²) */
- mp_multiply(&real_x, &real_x, &t1);
- mp_multiply(&im_x, &im_x, &t2);
- mp_add(&t1, &t2, &t1);
- mp_reciprocal_real(&t1, z);
- mp_conjugate(x, &t1);
- mp_multiply(&t1, z, z);
- }
- else
- mp_reciprocal_real(x, z);
+ mpc_set_si(z->num, 1, MPC_RNDNN);
+ mpc_fr_div(z->num, mpc_realref(z->num), x->num, MPC_RNDNN);
}
-
-static void
-mp_root_real(const MPNumber *x, int64_t n, MPNumber *z)
+void
+mp_root(const MPNumber *x, int64_t n, MPNumber *z)
{
- float approximation;
- int ex, np, it0, t;
- MPNumber t1, t2;
-
- /* x^(1/1) = x */
- if (n == 1) {
- mp_set_from_mp(x, z);
- return;
- }
+ uint64_t p;
- /* x^(1/0) invalid */
- if (n == 0) {
- mperr(_("Root must be non-zero"));
- 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_BASE, 64)) {
- mperr("*** ABS(N) TOO LARGE IN CALL TO MP_ROOT ***");
- mp_set_from_integer(0, z);
- return;
+ if (n < 0)
+ {
+ if (n == INT64_MIN)
+ p = (uint64_t) INT64_MAX + 1;
+ else
+ p = -n;
}
-
- /* 0^(1/n) = 0 for positive n */
- if (mp_is_zero(x)) {
- mp_set_from_integer(0, z);
- if (n <= 0)
- mperr(_("Negative root of zero is undefined"));
- return;
+ else if (n > 0)
+ {
+ mp_set_from_mp(x, z);
+ p = n;
}
-
- // FIXME: Imaginary root
- if (x->sign < 0 && np % 2 == 0) {
- mperr(_("nth root of negative number is undefined for even n"));
+ else
+ { /* Translators: Error displayed when attempting to take zeroth root */
+ mperr(_("The zeroth root of a number is undefined"));
mp_set_from_integer(0, z);
return;
}
-
- /* DIVIDE EXPONENT BY NP */
- ex = x->exponent / np;
-
- /* Get initial approximation */
- mp_set_from_mp(x, &t1);
- t1.exponent = 0;
- approximation = exp(((float)(np * ex - x->exponent) * log((float)MP_BASE) -
- log((fabs(mp_cast_to_float(&t1))))) / (float)np);
- mp_set_from_float(approximation, &t1);
- t1.sign = x->sign;
- t1.exponent -= ex;
-
- /* MAIN ITERATION LOOP */
- it0 = t = 3;
- while(1) {
- int ts2, ts3;
-
- /* t1 = t1 - ((t1 * ((x * t1^np) - 1)) / np) */
- mp_xpowy_integer(&t1, np, &t2);
- mp_multiply(x, &t2, &t2);
- mp_add_integer(&t2, -1, &t2);
- mp_multiply(&t1, &t2, &t2);
- mp_divide_integer(&t2, np, &t2);
- mp_subtract(&t1, &t2, &t1);
-
- /* FOLLOWING LOOP ALMOST DOUBLES T (POSSIBLE BECAUSE
- * NEWTONS METHOD HAS 2ND ORDER CONVERGENCE).
- */
- if (t >= MP_T)
- break;
-
- ts3 = t;
- t = MP_T;
- do {
- ts2 = t;
- t = (t + it0) / 2;
- } while (t > ts3);
- 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) {
- /* 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.
- */
- mperr("*** ERROR OCCURRED IN MP_ROOT, NEWTON ITERATION NOT CONVERGING PROPERLY ***");
- }
-
- if (n >= 0) {
- mp_xpowy_integer(&t1, n - 1, &t1);
- mp_multiply(x, &t1, z);
- return;
- }
-
- mp_set_from_mp(&t1, z);
-}
-
-
-void
-mp_root(const MPNumber *x, int64_t n, MPNumber *z)
-{
- if (!mp_is_complex(x) && mp_is_negative(x) && n % 2 == 1) {
- mp_abs(x, z);
- mp_root_real(z, n, z);
- mp_invert_sign(z, z);
- }
- else if (mp_is_complex(x) || mp_is_negative(x)) {
- MPNumber r, theta;
-
- mp_abs(x, &r);
- mp_arg(x, MP_RADIANS, &theta);
-
- mp_root_real(&r, n, &r);
- mp_divide_integer(&theta, n, &theta);
- mp_set_from_polar(&r, MP_RADIANS, &theta, z);
+ if (!mp_is_complex(x) && (!mp_is_negative(x) || (p & 1) == 1))
+ {
+ mpfr_rootn_ui(mpc_realref(z->num), mpc_realref(z->num), (uint64_t) p, MPFR_RNDN);
+ mpfr_set_zero(mpc_imagref(z->num), MPFR_RNDN);
}
else
- mp_root_real(x, n, z);
+ {
+ mpfr_t tmp;
+ mpfr_init2(tmp, PRECISION);
+ mpfr_set_ui(tmp, (uint64_t) p, MPFR_RNDN);
+ mpfr_ui_div(tmp, 1, tmp, MPFR_RNDN);
+ mpc_pow_fr(z->num, z->num, tmp, MPC_RNDNN);
+ mpfr_clear(tmp);
+ }
}
void
mp_sqrt(const MPNumber *x, MPNumber *z)
{
- if (mp_is_zero(x))
- mp_set_from_integer(0, z);
- /* FIXME: Make complex numbers optional */
- /*else if (x->sign < 0) {
- mperr(_("Square root is undefined for negative values"));
- mp_set_from_integer(0, z);
- }*/
- else {
- MPNumber t;
-
- mp_root(x, -2, &t);
- mp_multiply(x, &t, z);
- mp_ext(t.fraction[0], z->fraction[0], z);
- }
+ mp_root(x, 2, z);
}
void
mp_factorial(const MPNumber *x, MPNumber *z)
{
- int i, value;
-
/* 0! == 1 */
- if (mp_is_zero(x)) {
- mp_set_from_integer(1, z);
+ if (mp_is_zero(x))
+ {
+ mpc_set_si(z->num, 1, MPC_RNDNN);
return;
}
- if (!mp_is_natural(x)) {
- /* Translators: Error displayed when attempted take the factorial of a fractional number */
- mperr(_("Factorial is only defined for natural numbers"));
- mp_set_from_integer(0, z);
- return;
+ if (!mp_is_natural(x))
+ {
+ /* Factorial Not defined for Complex or for negative numbers */
+ if (mp_is_negative(x) || mp_is_complex(x))
+ { /* Translators: Error displayed when attempted take the factorial of a negative or complex number */
+ mperr(_("Factorial is only defined for non-negative real numbers"));
+ mp_set_from_integer(0, z);
+ return;
+ }
+ MPNumber tmp = mp_new();
+ mpfr_t tmp2;
+ mpfr_init2(tmp2, PRECISION);
+ mp_set_from_integer(1, &tmp);
+ mp_add(&tmp, x, &tmp);
+
+ /* Factorial(x) = Gamma(x+1) - This is the formula used to calculate Factorial of positive real numbers.*/
+ mpfr_gamma(tmp2, mpc_realref(tmp.num), MPFR_RNDN);
+ mpc_set_fr(z->num, tmp2, MPC_RNDNN);
+ mp_clear(&tmp);
+ mpfr_clear(tmp2);
+ }
+ else
+ {
+ /* Convert to integer - if couldn't be converted then the factorial would be too big anyway */
+ int64_t value = mp_to_integer(x);
+ mpfr_fac_ui(mpc_realref(z->num), value, MPFR_RNDN);
+ mpfr_set_zero(mpc_imagref(z->num), MPFR_RNDN);
}
-
- /* Convert to integer - if couldn't be converted then the factorial would be too big anyway */
- value = mp_cast_to_int(x);
- mp_set_from_mp(x, z);
- for (i = 2; i < value; i++)
- mp_multiply_integer(z, i, z);
}
void
mp_modulus_divide(const MPNumber *x, const MPNumber *y, MPNumber *z)
{
- MPNumber t1, t2;
-
- if (!mp_is_integer(x) || !mp_is_integer(y)) {
- /* Translators: Error displayed when attemping to do a modulus division on non-integer numbers */
+ if (!mp_is_integer(x) || !mp_is_integer(y))
+ { /* Translators: Error displayed when attemping to do a modulus division on non-integer numbers */
mperr(_("Modulus division is only defined for integers"));
mp_set_from_integer(0, z);
+ return;
}
- mp_divide(x, y, &t1);
- mp_floor(&t1, &t1);
- mp_multiply(&t1, y, &t2);
- mp_subtract(x, &t2, z);
+ mp_set_from_mp(x, z);
+ mp_divide(x, y, z);
+ mp_floor(z, z);
+ mp_multiply(z, y, z);
+ mp_subtract(x, z, z);
+ MPNumber t1 = mp_new();
mp_set_from_integer(0, &t1);
- if ((mp_is_less_than(y, &t1) && mp_is_greater_than(z, &t1)) || mp_is_less_than(z, &t1))
+ if ((mp_compare(y, &t1) < 0 && mp_compare(z, &t1) > 0) ||
+ (mp_compare(y, &t1) > 0 && mp_compare(z, &t1) < 0))
mp_add(z, y, z);
+
+ mp_clear(&t1);
}
void
mp_xpowy(const MPNumber *x, const MPNumber *y, MPNumber *z)
{
- if (mp_is_integer(y)) {
- mp_xpowy_integer(x, mp_cast_to_int(y), z);
- } else {
- MPNumber reciprocal;
+ /* 0^-n invalid */
+ if (mp_is_zero(x) && mp_is_negative(y))
+ { /* Translators: Error displayed when attempted to raise 0 to a negative exponent */
+ mperr(_("The power of zero is undefined for a negative exponent"));
+ mp_set_from_integer(0, z);
+ return;
+ }
+
+ if (!mp_is_complex(x) && !mp_is_complex(y) && !mp_is_integer(y))
+ {
+ MPNumber reciprocal = mp_new();
mp_reciprocal(y, &reciprocal);
if (mp_is_integer(&reciprocal))
- mp_root(x, mp_cast_to_int(&reciprocal), z);
- else
- mp_pwr(x, y, z);
+ {
+ mp_root(x, mp_to_integer(&reciprocal), z);
+ mp_clear(&reciprocal);
+ return;
+ }
+ mp_clear(&reciprocal);
}
+
+ mpc_pow(z->num, x->num, y->num, MPC_RNDNN);
}
void
mp_xpowy_integer(const MPNumber *x, int64_t n, MPNumber *z)
{
- int i;
- MPNumber t;
-
/* 0^-n invalid */
- if (mp_is_zero(x) && n < 0) {
- /* Translators: Error displayed when attempted to raise 0 to a negative exponent */
+ if (mp_is_zero(x) && n < 0)
+ { /* Translators: Error displayed when attempted to raise 0 to a negative re_exponent */
mperr(_("The power of zero is undefined for a negative exponent"));
mp_set_from_integer(0, z);
return;
}
- /* x^0 = 1 */
- if (n == 0) {
- mp_set_from_integer(1, z);
- return;
- }
-
- /* 0^n = 0 */
- if (mp_is_zero(x)) {
- mp_set_from_integer(0, z);
- return;
- }
-
- /* x^1 = x */
- if (n == 1) {
- mp_set_from_mp(x, z);
- return;
- }
-
- if (n < 0) {
- mp_reciprocal(x, &t);
- n = -n;
- }
- else
- mp_set_from_mp(x, &t);
-
- /* Multply x n times */
- // FIXME: Can do mp_multiply(z, z, z) until close to answer (each call doubles number of multiples) */
- mp_set_from_integer(1, z);
- for (i = 0; i < n; i++)
- mp_multiply(z, &t, z);
+ mpc_pow_si(z->num, x->num, (long) n, MPC_RNDNN);
}
-
GList*
mp_factorize(const MPNumber *x)
{
GList *list = NULL;
MPNumber *factor = g_slice_alloc0(sizeof(MPNumber));
- MPNumber value, tmp, divisor, root;
+ MPNumber value = mp_new();
+ MPNumber tmp = mp_new();
+ MPNumber divisor = mp_new();
+ MPNumber root = mp_new();
+ MPNumber int_max = mp_new();
+ mpc_init2(factor->num, PRECISION);
mp_abs(x, &value);
- if (mp_is_zero(&value)) {
+ if (mp_is_zero(&value))
+ {
mp_set_from_mp(&value, factor);
list = g_list_append(list, factor);
return list;
}
mp_set_from_integer(1, &tmp);
- if (mp_is_equal(&value, &tmp)) {
+ if (mp_is_equal(&value, &tmp))
+ {
mp_set_from_mp(x, factor);
list = g_list_append(list, factor);
return list;
}
+ /* If value < 2^64-1, call for factorize_uint64 function which deals in integers */
+ uint64_t num = 1;
+ num = num << 63;
+ num += (num - 1);
+ mp_set_from_unsigned_integer(num, &int_max);
+ if (mp_is_less_equal(x, &int_max))
+ {
+ list = mp_factorize_unit64(mp_to_unsigned_integer(&value));
+ if (mp_is_negative(x))
+ mp_invert_sign(list->data, list->data);
+ mp_clear(&int_max);
+ return list;
+ }
+
mp_set_from_integer(2, &divisor);
- while (TRUE) {
+ while (TRUE)
+ {
mp_divide(&value, &divisor, &tmp);
- if (mp_is_integer(&tmp)) {
- value = tmp;
+ if (mp_is_integer(&tmp))
+ {
+ mp_set_from_mp(&tmp, &value);
mp_set_from_mp(&divisor, factor);
list = g_list_append(list, factor);
factor = g_slice_alloc0(sizeof(MPNumber));
- } else {
- break;
+ mpc_init2(factor->num, PRECISION);
}
+ else
+ break;
}
mp_set_from_integer(3, &divisor);
mp_sqrt(&value, &root);
- while (mp_is_less_equal(&divisor, &root)) {
+ while (mp_is_less_equal(&divisor, &root))
+ {
mp_divide(&value, &divisor, &tmp);
- if (mp_is_integer(&tmp)) {
- value = tmp;
+ if (mp_is_integer(&tmp))
+ {
+ mp_set_from_mp(&tmp, &value);
mp_sqrt(&value, &root);
mp_set_from_mp(&divisor, factor);
list = g_list_append(list, factor);
factor = g_slice_alloc0(sizeof(MPNumber));
- } else {
+ mpc_init2(factor->num, PRECISION);
+ }
+ else
+ {
mp_add_integer(&divisor, 2, &tmp);
- divisor = tmp;
+ mp_set_from_mp(&tmp, &divisor);
}
}
mp_set_from_integer(1, &tmp);
- if (mp_is_greater_than(&value, &tmp)) {
+ if (mp_is_greater_than(&value, &tmp))
+ {
mp_set_from_mp(&value, factor);
list = g_list_append(list, factor);
- } else {
+ }
+ else
+ {
+ mpc_clear(factor->num);
g_slice_free(MPNumber, factor);
}
- if (mp_is_negative(x)) {
+ if (mp_is_negative(x))
mp_invert_sign(list->data, list->data);
+
+ /* MPNumbers in GList will be freed in math_equation_factorize_real */
+ mp_clear(&value); mp_clear(&tmp); mp_clear(&divisor); mp_clear(&root);
+
+ return list;
+}
+
+GList*
+mp_factorize_unit64(uint64_t n)
+{
+ GList *list = NULL;
+ MPNumber *factor = g_slice_alloc0(sizeof(MPNumber));
+ MPNumber tmp = mp_new();
+ mpc_init2(factor->num, PRECISION);
+
+ mp_set_from_unsigned_integer(2, &tmp);
+ while (n % 2 == 0)
+ {
+ n /= 2;
+ mp_set_from_mp(&tmp, factor);
+ list = g_list_append(list, factor);
+ factor = g_slice_alloc0(sizeof(MPNumber));
+ mpc_init2(factor->num, PRECISION);
+ }
+
+ for (uint64_t divisor = 3; divisor <= n / divisor; divisor +=2)
+ {
+ while (n % divisor == 0)
+ {
+ n /= divisor;
+ mp_set_from_unsigned_integer(divisor, factor);
+ list = g_list_append(list, factor);
+ factor = g_slice_alloc0(sizeof(MPNumber));
+ mpc_init2(factor->num, PRECISION);
+ }
+ }
+
+ if (n > 1)
+ {
+ mp_set_from_unsigned_integer(n, factor);
+ list = g_list_append(list, factor);
+ }
+ else
+ {
+ mpc_clear(factor->num);
+ g_slice_free(MPNumber, factor);
}
+ mp_clear(&tmp);
return list;
}
diff --git a/src/mp.h b/src/mp.h
index 53943b8..1748c03 100644
--- a/src/mp.h
+++ b/src/mp.h
@@ -33,27 +33,21 @@
#include <stdbool.h>
#include <stdint.h>
#include <glib.h>
+#include <glib/gi18n.h>
+#include <mpfr.h>
+#include <mpc.h>
-/* Size of the multiple precision values */
-#define MP_SIZE 1000
+/* If we're not using GNU C, elide __attribute__ */
+#ifndef __GNUC__
+# define __attribute__(x) /*NOTHING*/
+#endif
-/* Base for numbers */
-#define MP_BASE 10000
+/* Precision of mpfr_t and mpc_t type objects */
+#define PRECISION 1000
-/* 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, im_sign;
-
- /* Exponent (to base MP_BASE) */
- int exponent, im_exponent;
-
- /* Normalized fraction */
- int fraction[MP_SIZE], im_fraction[MP_SIZE];
+ mpc_t num;
} MPNumber;
typedef enum
@@ -65,17 +59,28 @@ typedef enum
/* Returns error string or NULL if no error */
// FIXME: Global variable
-const char *mp_get_error(void);
+const char *mp_get_error(void);
/* Clear any current error */
-void mp_clear_error(void);
+void mp_clear_error(void);
+
+void mperr(const char *format, ...) __attribute__((format(printf, 1, 2)));
+
+/* Returns initialized MPNumber object */
+MPNumber mp_new(void);
+
+MPNumber* mp_new_ptr(void);
+
+void mp_clear(MPNumber *z);
+
+void mp_free(MPNumber *z);
/* Returns:
* 0 if x == y
* <0 if x < y
* >0 if x > y
*/
-int mp_compare_mp_to_mp(const MPNumber *x, const MPNumber *y);
+int mp_compare(const MPNumber *x, const MPNumber *y);
/* Return true if the value is x == 0 */
bool mp_is_zero(const MPNumber *x);
@@ -134,9 +139,6 @@ void mp_add(const MPNumber *x, const MPNumber *y, MPNumber *z);
/* Sets z = x + y */
void mp_add_integer(const MPNumber *x, int64_t y, MPNumber *z);
-/* Sets z = x + numerator ÷ denominator */
-void mp_add_fraction(const MPNumber *x, int64_t numerator, int64_t denominator, MPNumber *z);
-
/* Sets z = x − y */
void mp_subtract(const MPNumber *x, const MPNumber *y, MPNumber *z);
@@ -146,9 +148,6 @@ void mp_multiply(const MPNumber *x, const MPNumber *y, MPNumber *z);
/* Sets z = x × y */
void mp_multiply_integer(const MPNumber *x, int64_t y, MPNumber *z);
-/* Sets z = x × numerator ÷ denominator */
-void mp_multiply_fraction(const MPNumber *x, int64_t numerator, int64_t denominator, MPNumber *z);
-
/* Sets z = x ÷ y */
void mp_divide(const MPNumber *x, const MPNumber *y, MPNumber *z);
@@ -217,6 +216,8 @@ void mp_epowy(const MPNumber *x, MPNumber *z);
/* Returns a list of all prime factors in x as MPNumbers */
GList* mp_factorize(const MPNumber *x);
+GList* mp_factorize_unit64 (uint64_t n);
+
/* Sets z = x */
void mp_set_from_mp(const MPNumber *x, MPNumber *z);
@@ -250,16 +251,16 @@ void mp_set_from_random(MPNumber *z);
bool mp_set_from_string(const char *text, int default_base, MPNumber *z);
/* Returns x as a native single-precision floating point number */
-float mp_cast_to_float(const MPNumber *x);
+float mp_to_float(const MPNumber *x);
/* Returns x as a native double-precision floating point number */
-double mp_cast_to_double(const MPNumber *x);
+double mp_to_double(const MPNumber *x);
/* Returns x as a native integer */
-int64_t mp_cast_to_int(const MPNumber *x);
+int64_t mp_to_integer(const MPNumber *x);
/* Returns x as a native unsigned integer */
-uint64_t mp_cast_to_unsigned_int(const MPNumber *x);
+uint64_t mp_to_unsigned_integer(const MPNumber *x);
/* Sets z = sin x */
void mp_sin(const MPNumber *x, MPAngleUnit unit, MPNumber *z);
@@ -324,4 +325,8 @@ void mp_ones_complement(const MPNumber *x, int wordlen, MPNumber *z);
/* Sets z to be the twos complement of x for word of length 'wordlen' */
void mp_twos_complement(const MPNumber *x, int wordlen, MPNumber *z);
+void convert_to_radians(const MPNumber *x, MPAngleUnit unit, MPNumber *z);
+
+void convert_from_radians(const MPNumber *x, MPAngleUnit unit, MPNumber *z);
+
#endif /* MP_H */
diff --git a/src/parser.c b/src/parser.c
index b4f90fc..633bc2f 100644
--- a/src/parser.c
+++ b/src/parser.c
@@ -276,8 +276,9 @@ p_parse(ParserState* state)
ans = (MPNumber *) (*(state->root->evaluate))(state->root);
if(ans)
{
+ state->ret = mp_new();
mp_set_from_mp(ans, &state->ret);
- free(ans);
+ mp_free(ans);
return PARSER_ERR_NONE;
}
return PARSER_ERR_INVALID;
@@ -330,10 +331,11 @@ p_check_variable(ParserState* state, gchar* name)
{
return FALSE;
}
-
+ temp = mp_new();
/* If defined, then get the variable */
if((*(state->get_variable))(state, name, &temp))
{
+ mp_clear(&temp);
return TRUE;
}
@@ -354,6 +356,7 @@ p_check_variable(ParserState* state, gchar* name)
}
free(buffer);
}
+ mp_clear(&temp);
if(!result)
{
return FALSE;
diff --git a/src/parserfunc.c b/src/parserfunc.c
index f35e6c5..8f153aa 100644
--- a/src/parserfunc.c
+++ b/src/parserfunc.c
@@ -31,7 +31,7 @@ pf_set_var(ParseNode* self)
if(!val || !(self->state->set_variable))
{
if(val)
- free(val);
+ mp_free(val);
return NULL;
}
(*(self->state->set_variable))(self->state, self->left->token->string, val);
@@ -46,9 +46,8 @@ pf_convert_number(ParseNode* self)
gchar* to;
gint free_from = 0;
gint free_to = 0;
- MPNumber tmp;
- MPNumber* ans;
- ans = (MPNumber *) malloc(sizeof(MPNumber));
+ MPNumber tmp = mp_new();
+ MPNumber* ans = mp_new_ptr();
if(self->left->value)
{
from = (gchar*) self->left->value;
@@ -66,20 +65,20 @@ pf_convert_number(ParseNode* self)
if(mp_set_from_string(self->left->left->token->string, self->state->options->base, &tmp) != 0)
{
- free(ans);
+ mp_free(ans);
ans = NULL;
goto END_PF_CONVERT_NUMBER;
}
if(!(self->state->convert))
{
- free(ans);
+ mp_free(ans);
ans = NULL;
goto END_PF_CONVERT_NUMBER;
}
if(!(*(self->state->convert))(self->state, &tmp, from, to, ans))
{
set_error(self->state, PARSER_ERR_UNKNOWN_CONVERSION, NULL);
- free(ans);
+ mp_free(ans);
ans = NULL;
}
END_PF_CONVERT_NUMBER:
@@ -93,6 +92,7 @@ END_PF_CONVERT_NUMBER:
g_free(self->right->value);
self->right->value = NULL;
}
+ mp_clear(&tmp);
return ans;
}
@@ -104,9 +104,8 @@ pf_convert_1(ParseNode* self )
gchar* to;
gint free_from = 0;
gint free_to = 0;
- MPNumber tmp;
- MPNumber* ans;
- ans = (MPNumber *) malloc(sizeof(MPNumber));
+ MPNumber tmp = mp_new();
+ MPNumber* ans = mp_new_ptr();
if(self->left->value)
{
from = (gchar*) self->left->value;
@@ -124,13 +123,13 @@ pf_convert_1(ParseNode* self )
mp_set_from_integer(1, &tmp);
if(!(self->state->convert))
{
- free(ans);
+ mp_free(ans);
return NULL;
}
if(!(*(self->state->convert))(self->state, &tmp, from, to, ans))
{
set_error(self->state, PARSER_ERR_UNKNOWN_CONVERSION, NULL);
- free(ans);
+ mp_free(ans);
ans = NULL;
}
if(free_from)
@@ -143,6 +142,7 @@ pf_convert_1(ParseNode* self )
g_free(self->right->value);
self->right->value = NULL;
}
+ mp_clear(&tmp);
return ans;
}
@@ -170,11 +170,10 @@ pf_get_variable(ParseNode* self)
const gchar *c, *next;
gchar *buffer;
- MPNumber value;
+ MPNumber value = mp_new();
- MPNumber t;
- MPNumber* ans;
- ans = (MPNumber*) malloc(sizeof(MPNumber));
+ MPNumber t = mp_new();
+ MPNumber* ans = mp_new_ptr();
if(!(self->state->get_variable))
{
@@ -227,11 +226,10 @@ pf_get_variable_with_power(ParseNode* self)
const gchar *c, *next;
gchar *buffer;
- MPNumber value;
+ MPNumber value = mp_new();
- MPNumber t;
- MPNumber* ans;
- ans = (MPNumber*) malloc(sizeof(MPNumber));
+ MPNumber t = mp_new();
+ MPNumber* ans = mp_new_ptr();
pow = super_atoi(((LexerToken*) self->value)->string);
/* No need to free the memory. It is allocated and freed somewhere else. */
@@ -289,8 +287,7 @@ void*
pf_apply_func(ParseNode* self)
{
MPNumber* val;
- MPNumber* ans;
- ans = (MPNumber*) malloc(sizeof(MPNumber));
+ MPNumber* ans = mp_new_ptr();
val = (MPNumber*) (*(self->right->evaluate))(self->right);
if(!(self->state->get_function))
{
@@ -319,40 +316,38 @@ void*
pf_apply_func_with_power(ParseNode* self)
{
MPNumber* val;
- MPNumber* tmp;
- MPNumber* ans;
+ MPNumber* tmp = mp_new_ptr();
+ MPNumber* ans = mp_new_ptr();
gint pow;
- tmp = (MPNumber*) malloc(sizeof(MPNumber));
- ans = (MPNumber*) malloc(sizeof(MPNumber));
val = (MPNumber*) (*(self->right->evaluate))(self->right);
if(!(self->state->get_function))
{
- free(tmp);
- free(ans);
- free(val);
+ mp_free(tmp);
+ mp_free(ans);
+ mp_free(val);
self->value = NULL;
return NULL;
}
if(!val)
{
- free(tmp);
- free(ans);
+ mp_free(tmp);
+ mp_free(ans);
self->value = NULL;
return NULL;
}
if(!(*(self->state->get_function))(self->state, self->token->string, val, tmp))
{
- free(tmp);
- free(ans);
- free(val);
+ mp_free(tmp);
+ mp_free(ans);
+ mp_free(val);
self->value = NULL;
set_error(self->state, PARSER_ERR_UNKNOWN_FUNCTION, self->token->string);
return NULL;
}
pow = super_atoi(((LexerToken*) self->value)->string);
mp_xpowy_integer(tmp, pow, ans);
- free(val);
- free(tmp);
+ mp_free(val);
+ mp_free(tmp);
self->value = NULL;
return ans;
}
@@ -362,37 +357,35 @@ void*
pf_apply_func_with_npower(ParseNode* self)
{
MPNumber* val;
- MPNumber* tmp;
- MPNumber* ans;
+ MPNumber* tmp = mp_new_ptr();
+ MPNumber* ans = mp_new_ptr();
gint pow;
gchar* inv_name;
inv_name = (gchar*) malloc(sizeof(gchar) * strlen(self->token->string) + strlen("⁻¹") + 1);
strcpy(inv_name, self->token->string);
strcat(inv_name, "⁻¹");
- tmp = (MPNumber*) malloc(sizeof(MPNumber));
- ans = (MPNumber*) malloc(sizeof(MPNumber));
val = (MPNumber*) (*(self->right->evaluate))(self->right);
if(!val)
{
- free(tmp);
+ mp_free(tmp);
free(inv_name);
- free(ans);
+ mp_free(ans);
self->value = NULL;
return NULL;
}
if(!(self->state->get_function))
{
- free(tmp);
- free(ans);
+ mp_free(tmp);
+ mp_free(ans);
free(inv_name);
self->value = NULL;
return NULL;
}
if(!(*(self->state->get_function))(self->state, inv_name, val, tmp))
{
- free(tmp);
- free(ans);
- free(val);
+ mp_free(tmp);
+ mp_free(ans);
+ mp_free(val);
free(inv_name);
self->value = NULL;
set_error(self->state, PARSER_ERR_UNKNOWN_FUNCTION, self->token->string);
@@ -400,8 +393,8 @@ pf_apply_func_with_npower(ParseNode* self)
}
pow = super_atoi(((LexerToken*) self->value)->string);
mp_xpowy_integer(tmp, -pow, ans);
- free(val);
- free(tmp);
+ mp_free(val);
+ mp_free(tmp);
free(inv_name);
self->value = NULL;
return ans;
@@ -413,18 +406,17 @@ pf_do_nth_root(ParseNode* self)
{
MPNumber* val;
gint pow;
- MPNumber* ans;
+ MPNumber* ans = mp_new_ptr();
pow = sub_atoi(((LexerToken*) self->value)->string);
self->value = NULL;
- ans = (MPNumber*) malloc(sizeof(MPNumber));
val = (MPNumber*) (*(self->right->evaluate))(self->right);
if(!val)
{
- free(ans);
+ mp_free(ans);
return NULL;
}
mp_root(val, pow, ans);
- free(val);
+ mp_free(val);
return ans;
}
@@ -433,8 +425,7 @@ void*
pf_do_sqrt(ParseNode* self)
{
MPNumber* val;
- MPNumber* ans;
- ans = (MPNumber*) malloc(sizeof(MPNumber));
+ MPNumber* ans = mp_new_ptr();
val = (MPNumber*) (*(self->right->evaluate))(self->right);
if(!val)
{
@@ -451,16 +442,15 @@ void*
pf_do_root_3(ParseNode* self)
{
MPNumber* val;
- MPNumber* ans;
- ans = (MPNumber*) malloc(sizeof(MPNumber));
+ MPNumber* ans = mp_new_ptr();
val = (MPNumber*) (*(self->right->evaluate))(self->right);
if(!val)
{
- free(ans);
+ mp_free(ans);
return NULL;
}
mp_root(val, 3, ans);
- free(val);
+ mp_free(val);
return ans;
}
@@ -469,16 +459,15 @@ void*
pf_do_root_4(ParseNode* self)
{
MPNumber* val;
- MPNumber* ans;
- ans = (MPNumber*) malloc(sizeof(MPNumber));
+ MPNumber* ans = mp_new_ptr();
val = (MPNumber*) (*(self->right->evaluate))(self->right);
if(!val)
{
- free(ans);
+ mp_free(ans);
return NULL;
}
mp_root(val, 4, ans);
- free(val);
+ mp_free(val);
return ans;
}
@@ -487,16 +476,15 @@ void*
pf_do_floor(ParseNode* self)
{
MPNumber* val;
- MPNumber* ans;
- ans = (MPNumber*) malloc(sizeof(MPNumber));
+ MPNumber* ans = mp_new_ptr();
val = (MPNumber*) (*(self->right->evaluate))(self->right);
if(!val)
{
- free(ans);
+ mp_free(ans);
return NULL;
}
mp_floor(val, ans);
- free(val);
+ mp_free(val);
return ans;
}
@@ -504,16 +492,15 @@ pf_do_floor(ParseNode* self)
void* pf_do_ceiling (ParseNode* self)
{
MPNumber* val;
- MPNumber* ans;
- ans = (MPNumber*) malloc(sizeof(MPNumber));
+ MPNumber* ans = mp_new_ptr();
val = (MPNumber*) (*(self->right->evaluate))(self->right);
if(!val)
{
- free(ans);
+ mp_free(ans);
return NULL;
}
mp_ceiling(val, ans);
- free(val);
+ mp_free(val);
return ans;
}
@@ -522,16 +509,15 @@ void*
pf_do_round(ParseNode* self)
{
MPNumber* val;
- MPNumber* ans;
- ans = (MPNumber*) malloc(sizeof(MPNumber));
+ MPNumber* ans = mp_new_ptr();
val = (MPNumber*) (*(self->right->evaluate))(self->right);
if(!val)
{
- free(ans);
+ mp_free(ans);
return NULL;
}
mp_round(val, ans);
- free(val);
+ mp_free(val);
return ans;
}
@@ -540,16 +526,15 @@ void*
pf_do_fraction(ParseNode* self)
{
MPNumber* val;
- MPNumber* ans;
- ans = (MPNumber*) malloc(sizeof(MPNumber));
+ MPNumber* ans = mp_new_ptr();
val = (MPNumber*) (*(self->right->evaluate))(self->right);
if(!val)
{
- free(ans);
+ mp_free(ans);
return NULL;
}
mp_fractional_part(val, ans);
- free(val);
+ mp_free(val);
return ans;
}
@@ -558,16 +543,15 @@ void*
pf_do_abs(ParseNode* self)
{
MPNumber* val;
- MPNumber* ans;
- ans = (MPNumber*) malloc(sizeof(MPNumber));
+ MPNumber* ans = mp_new_ptr();
val = (MPNumber*) (*(self->right->evaluate))(self->right);
if(!val)
{
- free(ans);
+ mp_free(ans);
return NULL;
}
mp_abs(val, ans);
- free(val);
+ mp_free(val);
return ans;
}
@@ -577,22 +561,21 @@ pf_do_x_pow_y(ParseNode* self)
{
MPNumber* val;
MPNumber* pow;
- MPNumber* ans;
- ans = (MPNumber*) malloc(sizeof(MPNumber));
+ MPNumber* ans = mp_new_ptr();
val = (MPNumber*) (*(self->left->evaluate))(self->left);
pow = (MPNumber*) (*(self->right->evaluate))(self->right);
if(!val || !pow)
{
if(val)
- free(val);
+ mp_free(val);
if(pow)
- free(pow);
- free(ans);
+ mp_free(pow);
+ mp_free(ans);
return NULL;
}
mp_xpowy(val, pow, ans);
- free(val);
- free(pow);
+ mp_free(val);
+ mp_free(pow);
return ans;
}
@@ -602,17 +585,16 @@ pf_do_x_pow_y_int(ParseNode* self)
{
MPNumber* val;
gint pow;
- MPNumber* ans;
- ans = (MPNumber*) malloc(sizeof(MPNumber));
+ MPNumber* ans = mp_new_ptr();
val = (MPNumber*) (*(self->left->evaluate))(self->left);
pow = super_atoi(self->right->token->string);
if(!val)
{
- free(ans);
+ mp_free(ans);
return NULL;
}
mp_xpowy_integer(val, pow, ans);
- free(val);
+ mp_free(val);
return ans;
}
@@ -621,16 +603,15 @@ void*
pf_do_factorial(ParseNode* self)
{
MPNumber* val;
- MPNumber* ans;
- ans = (MPNumber*) malloc(sizeof(MPNumber));
+ MPNumber* ans = mp_new_ptr();
val = (MPNumber*) (*(self->right->evaluate))(self->right);
if(!val)
{
- free(ans);
+ mp_free(ans);
return NULL;
}
mp_factorial(val, ans);
- free(val);
+ mp_free(val);
return ans;
}
@@ -639,16 +620,15 @@ void*
pf_unary_minus(ParseNode* self)
{
MPNumber* val;
- MPNumber* ans;
- ans = (MPNumber*) malloc(sizeof(MPNumber));
+ MPNumber* ans = mp_new_ptr();
val = (MPNumber*) (*(self->right->evaluate))(self->right);
if(!val)
{
- free(ans);
+ mp_free(ans);
return NULL;
}
mp_invert_sign(val, ans);
- free(val);
+ mp_free(val);
return ans;
}
@@ -658,22 +638,21 @@ pf_do_divide(ParseNode* self)
{
MPNumber* left;
MPNumber* right;
- MPNumber* ans;
- ans = (MPNumber*) malloc(sizeof(MPNumber));
+ MPNumber* ans = mp_new_ptr();
left = (MPNumber*) (*(self->left->evaluate))(self->left);
right = (MPNumber*) (*(self->right->evaluate))(self->right);
if(!left || !right)
{
if(left)
- free(left);
+ mp_free(left);
if(right)
- free(right);
- free(ans);
+ mp_free(right);
+ mp_free(ans);
return NULL;
}
mp_divide(left, right, ans);
- free(left);
- free(right);
+ mp_free(left);
+ mp_free(right);
return ans;
}
@@ -683,22 +662,21 @@ pf_do_mod(ParseNode* self)
{
MPNumber* left;
MPNumber* right;
- MPNumber* ans;
- ans = (MPNumber*) malloc(sizeof(MPNumber));
+ MPNumber* ans = mp_new_ptr();
left = (MPNumber*) (*(self->left->evaluate))(self->left);
right = (MPNumber*) (*(self->right->evaluate))(self->right);
if(!left || !right)
{
if(left)
- free(left);
+ mp_free(left);
if(right)
- free(right);
+ mp_free(right);
free(ans);
return NULL;
}
mp_modulus_divide(left, right, ans);
- free(left);
- free(right);
+ mp_free(left);
+ mp_free(right);
return ans;
}
@@ -708,22 +686,21 @@ pf_do_multiply(ParseNode* self)
{
MPNumber* left;
MPNumber* right;
- MPNumber* ans;
- ans = (MPNumber*) malloc(sizeof(MPNumber));
+ MPNumber* ans = mp_new_ptr();
left = (MPNumber*) (*(self->left->evaluate))(self->left);
right = (MPNumber*) (*(self->right->evaluate))(self->right);
if(!left || !right)
{
if(left)
- free(left);
+ mp_free(left);
if(right)
- free(right);
- free(ans);
+ mp_free(right);
+ mp_free(ans);
return NULL;
}
mp_multiply(left, right, ans);
- free(left);
- free(right);
+ mp_free(left);
+ mp_free(right);
return ans;
}
@@ -733,22 +710,21 @@ pf_do_subtract(ParseNode* self)
{
MPNumber* left;
MPNumber* right;
- MPNumber* ans;
- ans = (MPNumber*) malloc(sizeof(MPNumber));
+ MPNumber* ans = mp_new_ptr();
left = (MPNumber*) (*(self->left->evaluate))(self->left);
right = (MPNumber*) (*(self->right->evaluate))(self->right);
if(!left || !right)
{
if(left)
- free(left);
+ mp_free(left);
if(right)
- free (right);
+ mp_free(right);
free(ans);
return NULL;
}
mp_subtract(left, right, ans);
- free(left);
- free(right);
+ mp_free(left);
+ mp_free(right);
return ans;
}
@@ -758,22 +734,21 @@ pf_do_add(ParseNode* self)
{
MPNumber* left;
MPNumber* right;
- MPNumber* ans;
- ans = (MPNumber*) malloc(sizeof(MPNumber));
+ MPNumber* ans = mp_new_ptr();
left = (MPNumber*) (*(self->left->evaluate))(self->left);
right = (MPNumber*) (*(self->right->evaluate))(self->right);
if(!left || !right)
{
if(left)
- free(left);
+ mp_free(left);
if(right)
- free(right);
+ mp_free(right);
free(ans);
return NULL;
}
mp_add(left, right, ans);
- free(left);
- free(right);
+ mp_free(left);
+ mp_free(right);
return ans;
}
@@ -781,26 +756,25 @@ pf_do_add(ParseNode* self)
void*
pf_do_add_percent(ParseNode* self)
{
- MPNumber* ans;
+ MPNumber* ans = mp_new_ptr();
MPNumber* val;
MPNumber* per;
- ans = (MPNumber*) malloc(sizeof(MPNumber));
val = (MPNumber*) (*(self->left->evaluate))(self->left);
per = (MPNumber*) (*(self->right->evaluate))(self->right);
if(!val || !per)
{
if(val)
- free(val);
+ mp_free(val);
if(per)
- free(per);
- free(ans);
+ mp_free(per);
+ mp_free(ans);
return NULL;
}
mp_add_integer(per, 100, per);
mp_divide_integer(per, 100, per);
mp_multiply(val, per, ans);
- free(val);
- free(per);
+ mp_free(val);
+ mp_free(per);
return ans;
}
@@ -808,26 +782,25 @@ pf_do_add_percent(ParseNode* self)
void*
pf_do_subtract_percent(ParseNode* self)
{
- MPNumber* ans;
+ MPNumber* ans = mp_new_ptr();
MPNumber* val;
MPNumber* per;
- ans = (MPNumber*) malloc(sizeof(MPNumber));
val = (MPNumber*) (*(self->left->evaluate))(self->left);
per = (MPNumber*) (*(self->right->evaluate))(self->right);
if(!val || !per)
{
if(val)
- free(val);
+ mp_free(val);
if(per)
- free(per);
+ mp_free(per);
free(ans);
return NULL;
}
mp_add_integer(per, -100, per);
mp_divide_integer(per, -100, per);
mp_multiply(val, per, ans);
- free(val);
- free(per);
+ mp_free(val);
+ mp_free(per);
return ans;
}
@@ -836,16 +809,15 @@ void*
pf_do_percent(ParseNode* self)
{
MPNumber* val;
- MPNumber* ans;
- ans = (MPNumber*) malloc(sizeof(MPNumber));
+ MPNumber* ans = mp_new_ptr();
val = (MPNumber*) (*(self->right->evaluate))(self->right);
if(!val)
{
- free(ans);
+ mp_free(ans);
return NULL;
}
mp_divide_integer(val, 100, ans);
- free(val);
+ mp_free(val);
return ans;
}
@@ -854,23 +826,22 @@ void*
pf_do_not(ParseNode* self)
{
MPNumber* val;
- MPNumber* ans;
- ans = (MPNumber*) malloc(sizeof(MPNumber));
+ MPNumber* ans = mp_new_ptr();
val = (MPNumber*) (*(self->right->evaluate))(self->right);
if(!val)
{
- free(ans);
+ mp_free(ans);
return NULL;
}
if(!mp_is_overflow(val, self->state->options->wordlen))
{
set_error(self->state, PARSER_ERR_OVERFLOW, NULL);
- free(ans);
- free(val);
+ mp_free(ans);
+ mp_free(val);
return NULL;
}
mp_not(val, self->state->options->wordlen, ans);
- free(val);
+ mp_free(val);
return ans;
}
@@ -880,22 +851,21 @@ pf_do_and(ParseNode* self)
{
MPNumber* left;
MPNumber* right;
- MPNumber* ans;
- ans = (MPNumber*) malloc(sizeof(MPNumber));
+ MPNumber* ans = mp_new_ptr();
left = (MPNumber*) (*(self->left->evaluate))(self->left);
right = (MPNumber*) (*(self->right->evaluate))(self->right);
if(!left || !right)
{
if(left)
- free(left);
+ mp_free(left);
if(right)
- free(right);
- free(ans);
+ mp_free(right);
+ mp_free(ans);
return NULL;
}
mp_and(left, right, ans);
- free(left);
- free(right);
+ mp_free(left);
+ mp_free(right);
return ans;
}
@@ -905,22 +875,21 @@ pf_do_or(ParseNode* self)
{
MPNumber* left;
MPNumber* right;
- MPNumber* ans;
- ans = (MPNumber*) malloc(sizeof(MPNumber));
+ MPNumber* ans = mp_new_ptr();
left = (MPNumber*) (*(self->left->evaluate))(self->left);
right = (MPNumber*) (*(self->right->evaluate))(self->right);
if(!left || !right)
{
if(left)
- free(left);
+ mp_free(left);
if(right)
- free(right);
- free(ans);
+ mp_free(right);
+ mp_free(ans);
return NULL;
}
mp_or(left, right, ans);
- free(left);
- free(right);
+ mp_free(left);
+ mp_free(right);
return ans;
}
@@ -930,22 +899,21 @@ pf_do_xor(ParseNode* self)
{
MPNumber* left;
MPNumber* right;
- MPNumber* ans;
- ans = (MPNumber*) malloc(sizeof(MPNumber));
+ MPNumber* ans = mp_new_ptr();
left = (MPNumber*) (*(self->left->evaluate))(self->left);
right = (MPNumber*) (*(self->right->evaluate))(self->right);
if(!left || !right)
{
if(left)
- free(left);
+ mp_free(left);
if(right)
- free(right);
+ mp_free(right);
free(ans);
return NULL;
}
mp_xor(left, right, ans);
- free(left);
- free(right);
+ mp_free(left);
+ mp_free(right);
return ans;
}
@@ -953,13 +921,12 @@ pf_do_xor(ParseNode* self)
void*
pf_constant(ParseNode* self)
{
- MPNumber* ans;
- ans = (MPNumber*) malloc(sizeof(MPNumber));
+ MPNumber* ans = mp_new_ptr();
if(mp_set_from_string(self->token->string, self->state->options->base, ans) != 0)
{
/* This should never happen, as l_check_if_number() has already passed the string once. */
/* If the code reaches this point, something is really wrong. X( */
- free(ans);
+ mp_free(ans);
set_error(self->state, PARSER_ERR_INVALID, self->token->string);
return NULL;
}
diff --git a/src/preferences.ui b/src/preferences.ui
index 56635b4..ca1e2da 100644
--- a/src/preferences.ui
+++ b/src/preferences.ui
@@ -11,7 +11,7 @@
</columns>
</object>
<object class="GtkAdjustment" id="decimal_places_adjustment">
- <property name="upper">15</property>
+ <property name="upper">100</property>
<property name="step_increment">1</property>
<property name="page_increment">1</property>
</object>
diff --git a/src/test-mp-equation.c b/src/test-mp-equation.c
index 0ff3a53..ac9b000 100644
--- a/src/test-mp-equation.c
+++ b/src/test-mp-equation.c
@@ -77,7 +77,7 @@ static void
Test(char *expression, char *expected, int expected_error, int trailing_digits)
{
MPErrorCode error;
- MPNumber result;
+ MPNumber result = mp_new();
error = mp_equation_parse(expression, &options, &result, NULL);
@@ -107,6 +107,7 @@ Test(char *expression, char *expected, int expected_error, int trailing_digits)
fail("'%s' -> error %s, expected error %s", expression,
error_code_to_string(error), error_code_to_string(expected_error));
}
+ mp_clear(&result);
}
@@ -425,7 +426,7 @@ test_equations()
test("1!", "1", 0);
test("5!", "120", 0);
test("69!", "171122452428141311372468338881272839092270544893520369393648040923257279754140647424000000000000000", 0);
- test("0.1!", "", PARSER_ERR_MP);
+ test("0.1!", "0.95135077", 0);
test("−1!", "−1", 0);
test("(−1)!", "", PARSER_ERR_MP);
test("−(1!)", "−1", 0);
diff --git a/src/test-mp.c b/src/test-mp.c
index fa7408f..78b1a5b 100644
--- a/src/test-mp.c
+++ b/src/test-mp.c
@@ -15,7 +15,6 @@
#include <locale.h>
#include "mp.h"
-#include "mp-private.h"
static int fails = 0;
static int passes = 0;
@@ -58,48 +57,39 @@ static void fail(const char *format, ...)
static void
print_number(MPNumber *x)
{
- int i, j;
-
- printf("sign=%d exponent=%d fraction=%d", x->sign, x->exponent, x->fraction[0]);
- for (i = 1; i < MP_SIZE; i++) {
- for (j = i; j < MP_SIZE && x->fraction[j] == 0; j++);
- if (j == MP_SIZE) {
- printf(",...");
- break;
- }
- printf(",%d", x->fraction[i]);
- }
+ mpc_out_str(stdout, 10, 5, x->num, MPC_RNDNN);
}
static void
test_string(const char *number)
{
- MPNumber t;
+ MPNumber t = mp_new();
mp_set_from_string(number, 10, &t);
printf("MPNumber(%s) -> {", number);
print_number(&t);
printf("}\n");
+ mp_clear(&t);
}
static void
test_integer(int number)
{
- MPNumber t;
+ MPNumber t = mp_new();
mp_set_from_integer(number, &t);
printf("MPNumber(%d) -> {", number);
print_number(&t);
printf("}\n");
+ mp_clear(&t);
}
static void
test_numbers()
{
- printf("base=%d\n", MP_BASE);
test_integer(0);
test_integer(1);
test_integer(-1);
@@ -143,8 +133,9 @@ try(const char *string, bool result, bool expected)
static void
test_mp()
{
- MPNumber zero, one, minus_one;
-
+ MPNumber zero = mp_new();
+ MPNumber one = mp_new();
+ MPNumber minus_one = mp_new();
mp_set_from_integer(0, &zero);
mp_set_from_integer(1, &one);
mp_set_from_integer(-1, &minus_one);
@@ -222,6 +213,8 @@ test_mp()
try("mp_is_less_equal(-1, -1)", mp_is_less_equal (&minus_one, &minus_one), true);
try("mp_is_less_equal(-1, 0)", mp_is_less_equal (&minus_one, &zero), true);
try("mp_is_less_equal(-1, 1)", mp_is_less_equal (&minus_one, &one), true);
+
+ mp_clear(&zero); mp_clear(&one); mp_clear(&minus_one);
}
diff --git a/src/unit-category.c b/src/unit-category.c
index b8fa303..72f50d5 100644
--- a/src/unit-category.c
+++ b/src/unit-category.c
@@ -106,18 +106,23 @@ unit_category_get_units(UnitCategory *category)
gboolean
unit_category_convert(UnitCategory *category, const MPNumber *x, Unit *x_units, Unit *z_units, MPNumber *z)
{
- MPNumber t;
-
g_return_val_if_fail (category != NULL, FALSE);
g_return_val_if_fail (x_units != NULL, FALSE);
g_return_val_if_fail (z_units != NULL, FALSE);
g_return_val_if_fail (z != NULL, FALSE);
+ MPNumber t = mp_new();
if (!unit_convert_from(x_units, x, &t))
+ {
+ mp_clear(&t);
return FALSE;
+ }
if (!unit_convert_to(z_units, &t, z))
+ {
+ mp_clear(&t);
return FALSE;
-
+ }
+ mp_clear(&t);
return TRUE;
}
diff --git a/src/unittest.c b/src/unittest.c
index d02ceea..4bbb436 100644
--- a/src/unittest.c
+++ b/src/unittest.c
@@ -77,7 +77,7 @@ static const char* error_code_to_string(MPErrorCode error)
static void test(char* expression, char* expected, int expected_error)
{
MPErrorCode error;
- MPNumber result;
+ MPNumber result = mp_new();
char result_str[1024] = "";
error = mp_equation_parse(expression, &options, &result, NULL);
@@ -110,6 +110,7 @@ static void test(char* expression, char* expected, int expected_error)
fail("'%s' -> error %s, expected error %s", expression, error_code_to_string(error), error_code_to_string(expected_error));
}
}
+ mp_clear(result);
}
@@ -594,54 +595,35 @@ static void test_equations()
static void print_number(MPNumber* x)
{
- int i, j;
-
- printf("sign=%d exponent=%d fraction=%d", x->sign, x->exponent, x->fraction[0]);
-
- for (i = 1; i < MP_SIZE; i++)
- {
- for (j = i; j < MP_SIZE && x->fraction[j] == 0; j++)
- {
- /* nothing*/
- }
-
- if (j == MP_SIZE)
- {
- printf(",...");
- break;
- }
-
- printf(",%d", x->fraction[i]);
- }
+ mpc_out_str(stdout, 10, 5, x->num, MPC_RNDNN);
}
static void test_string(const char* number)
{
- MPNumber t;
+ MPNumber t = mp_new();
mp_set_from_string(number, 10, &t);
printf("MPNumber(%s) -> {", number);
print_number(&t);
printf("}\n");
+ mp_clear(&t);
}
static void test_integer(int number)
{
- MPNumber t;
+ MPNumber t = mp_new();
mp_set_from_integer(number, &t);
printf("MPNumber(%d) -> {", number);
print_number(&t);
printf("}\n");
+ mp_clear(&t);
}
-#include "mp-private.h"
static void test_numbers()
{
- printf("base=%d\n", MP_BASE);
-
test_integer(0);
test_integer(1);
test_integer(-1);
@@ -687,11 +669,12 @@ static void try(const char* string, bool result, bool expected)
static void test_mp()
{
- MPNumber zero, one, minus_one;
-
- mp_set_from_integer(0, &zero);
- mp_set_from_integer(1, &one);
- mp_set_from_integer(-1, &minus_one);
+ MPNumber zero = mp_new();
+ MPNumber one = mp_new();
+ MPNumber minus_one = mp_new();
+ mp_set_from_integer(0, &zero);
+ mp_set_from_integer(1, &one);
+ mp_set_from_integer(-1, &minus_one);
try("mp_is_zero(-1)", mp_is_zero(&minus_one), false);
try("mp_is_zero(0)", mp_is_zero(&zero), true);
@@ -766,6 +749,8 @@ static void test_mp()
try("mp_is_less_equal(-1, -1)", mp_is_less_equal (&minus_one, &minus_one), true);
try("mp_is_less_equal(-1, 0)", mp_is_less_equal (&minus_one, &zero), true);
try("mp_is_less_equal(-1, 1)", mp_is_less_equal (&minus_one, &one), true);
+
+ mp_clear(&zero); mp_clear(&one); mp_clear(&minus_one);
}
generated by cgit v1.2.1 (git 2.18.0) at 2024-11-05 08:35:00 +0000