Port to GNU MPFR/MPC Library

For further information please visit:
https://www.mpfr.org/
http://www.multiprecision.org/mpc

26 files changed, 877 insertions, 2929 deletions

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