diff options
| -rw-r--r-- | configure.ac | 3 | ||||
| -rw-r--r-- | data/org.mate.calc.gschema.xml | 2 | ||||
| -rw-r--r-- | src/Makefile.am | 16 | ||||
| -rw-r--r-- | src/currency-manager.c | 13 | ||||
| -rw-r--r-- | src/currency.c | 1 | ||||
| -rw-r--r-- | src/financial.c | 68 | ||||
| -rw-r--r-- | src/lexer.c | 5 | ||||
| -rw-r--r-- | src/mate-calc-cmd.c | 3 | ||||
| -rw-r--r-- | src/mate-calc.c | 5 | ||||
| -rw-r--r-- | src/math-buttons.c | 20 | ||||
| -rw-r--r-- | src/math-converter.c | 10 | ||||
| -rw-r--r-- | src/math-equation.c | 33 | ||||
| -rw-r--r-- | src/math-variable-popup.c | 6 | ||||
| -rw-r--r-- | src/math-variables.c | 5 | ||||
| -rw-r--r-- | src/mp-binary.c | 31 | ||||
| -rw-r--r-- | src/mp-convert.c | 436 | ||||
| -rw-r--r-- | src/mp-private.h | 39 | ||||
| -rw-r--r-- | src/mp-serializer.c | 88 | ||||
| -rw-r--r-- | src/mp-trigonometric.c | 578 | ||||
| -rw-r--r-- | src/mp.c | 1962 | ||||
| -rw-r--r-- | src/mp.h | 63 | ||||
| -rw-r--r-- | src/parser.c | 7 | ||||
| -rw-r--r-- | src/parserfunc.c | 327 | ||||
| -rw-r--r-- | src/preferences.ui | 2 | ||||
| -rw-r--r-- | src/test-mp-equation.c | 5 | ||||
| -rw-r--r-- | src/test-mp.c | 27 | ||||
| -rw-r--r-- | src/unit-category.c | 11 | ||||
| -rw-r--r-- | src/unittest.c | 45 |
28 files changed, 881 insertions, 2930 deletions
diff --git a/configure.ac b/configure.ac index f64c607..8f57249 100644 --- a/configure.ac +++ b/configure.ac @@ -41,6 +41,9 @@ GLIB_MKENUMS=`$PKG_CONFIG --variable=glib_mkenums glib-2.0` AC_SUBST(GLIB_MKENUMS) AC_CHECK_LIB(m, log) +AC_CHECK_LIB(gmp, log) +AC_CHECK_LIB(mpfr, log) +AC_CHECK_LIB(mpc, log) dnl ########################################################################### dnl Internationalization diff --git a/data/org.mate.calc.gschema.xml b/data/org.mate.calc.gschema.xml index 507af68..34b17c3 100644 --- a/data/org.mate.calc.gschema.xml +++ b/data/org.mate.calc.gschema.xml @@ -21,7 +21,7 @@ <schema path="/org/mate/calc/" id="org.mate.calc" gettext-domain="mate-calc"> <key type="i" name="accuracy"> <default>9</default> - <range min="0" max="15"/> + <range min="0" max="100"/> <summary>Accuracy value</summary> <description>The number of digits displayed after the numeric point</description> </key> 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); } @@ -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; } @@ -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); } |