Add modular exponentiation ability and add acccording tests

4 files changed, 87 insertions, 14 deletions

diff --git a/src/mp.c b/src/mp.c
index c66285e..d7fcc75 100644
--- a/src/mp.c
+++ b/src/mp.c
@@ -383,7 +383,7 @@ void
 mp_logarithm(int64_t n, const MPNumber *x, MPNumber *z)
 {
     /* 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"));
@@ -525,19 +525,65 @@ mp_modulus_divide(const MPNumber *x, const MPNumber *y, MPNumber *z)
         return;
     }
 
-    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();
+    MPNumber t2 = mp_new();
+
+    mp_divide(x, y, &t1);
+    mp_floor(&t1, &t1);
+    mp_multiply(&t1, y, &t2);
+    mp_subtract(x, &t2, z);
+
     mp_set_from_integer(0, &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);
+    mp_clear(&t2);
+}
+
+
+void
+mp_modular_exponentiation(const MPNumber *x, const MPNumber *y, const MPNumber *p, MPNumber *z)
+{
+    MPNumber base_value = mp_new();
+    MPNumber exp_value = mp_new();
+    MPNumber ans = mp_new();
+    MPNumber two = mp_new();
+    MPNumber tmp = mp_new();
+
+    mp_set_from_integer(1, &ans);
+    mp_set_from_integer(2, &two);
+    mp_abs(y, &exp_value);
+
+    if (mp_is_negative(y))
+        mp_reciprocal(x, &base_value);
+    else
+        mp_set_from_mp(x, &base_value);
+
+    while (!mp_is_zero(&exp_value))
+    {
+        mp_modulus_divide(&exp_value, &two, &tmp);
+
+        bool is_even = mp_is_zero(&tmp);
+        if (!is_even)
+        {
+            mp_multiply(&ans, &base_value, &ans);
+            mp_modulus_divide(&ans, p, &ans);
+        }
+        mp_multiply(&base_value, &base_value, &base_value);
+        mp_modulus_divide(&base_value, p, &base_value);
+        mp_divide_integer(&exp_value, 2, &exp_value);
+        mp_floor(&exp_value, &exp_value);
+    }
+
+    mp_modulus_divide(&ans, p, z);
+
+    mp_clear(&base_value);
+    mp_clear(&exp_value);
+    mp_clear(&ans);
+    mp_clear(&two);
+    mp_clear(&tmp);
 }
 
 
@@ -694,7 +740,7 @@ mp_factorize(const MPNumber *x)
             list = g_list_append(list, factor);
             factor = g_slice_alloc0(sizeof(MPNumber));
             mpc_init2(factor->num, PRECISION);
-        } 
+        }
         else
         {
             mp_add_integer(&divisor, 2, &tmp);
@@ -707,8 +753,8 @@ mp_factorize(const MPNumber *x)
     {
         mp_set_from_mp(&value, factor);
         list = g_list_append(list, factor);
-    } 
-    else 
+    }
+    else
     {
         mpc_clear(factor->num);
         g_slice_free(MPNumber, factor);
@@ -758,7 +804,7 @@ mp_factorize_unit64(uint64_t n)
         mp_set_from_unsigned_integer(n, factor);
         list = g_list_append(list, factor);
     }
-    else 
+    else
     {
         mpc_clear(factor->num);
         g_slice_free(MPNumber, factor);
diff --git a/src/mp.h b/src/mp.h
index f57ff16..3e03089 100644
--- a/src/mp.h
+++ b/src/mp.h
@@ -204,6 +204,9 @@ void   mp_factorial(const MPNumber *x, MPNumber *z);
 /* Sets z = x mod y */
 void   mp_modulus_divide(const MPNumber *x, const MPNumber *y, MPNumber *z);
 
+/* Sets z = x ^ y mod p */
+void mp_modular_exponentiation(const MPNumber *x, const MPNumber *y, const MPNumber *p, MPNumber *z);
+
 /* Sets z = x^y */
 void   mp_xpowy(const MPNumber *x, const MPNumber *y, MPNumber *z);
 
diff --git a/src/parserfunc.c b/src/parserfunc.c
index ff5c6f4..c63ff96 100644
--- a/src/parserfunc.c
+++ b/src/parserfunc.c
@@ -688,10 +688,29 @@ pf_do_mod(ParseNode* self)
             mp_free(left);
         if(right)
             mp_free(right);
-        free(ans);
+        mp_free(ans);
         return NULL;
     }
-    mp_modulus_divide(left, right, ans);
+    if (self->left->evaluate == pf_do_x_pow_y)
+    {
+        MPNumber* base_value = (MPNumber*) (*(self->left->left->evaluate))(self->left->left);
+        MPNumber* exponent = (MPNumber*) (*(self->left->right->evaluate))(self->left->right);
+        if(!base_value || !exponent)
+        {
+            if(base_value)
+                mp_free(base_value);
+            if(exponent)
+                mp_free(exponent);
+            mp_free(ans);
+            return NULL;
+        }
+        mp_modular_exponentiation(base_value, exponent, right, ans);
+        mp_free(base_value);
+        mp_free(exponent);
+    }
+    else
+        mp_modulus_divide(left, right, ans);
+
     mp_free(left);
     mp_free(right);
     return ans;
@@ -769,7 +788,7 @@ pf_do_add(ParseNode* self)
     return ans;
 }
 
-/*Add (x) Percentage to value. */
+/* Add (x) Percentage to value. */
 void*
 pf_do_add_percent(ParseNode* self)
 {
diff --git a/src/test-mp-equation.c b/src/test-mp-equation.c
index c1c136a..c6ab98c 100644
--- a/src/test-mp-equation.c
+++ b/src/test-mp-equation.c
@@ -482,6 +482,11 @@ test_equations(void)
     test("8 mod 7", "1", 0);
     test("−1 mod 7", "6", 0);
 
+    test("123^456 mod 78", "27", 0);
+    test("2^2 mod 2", "0", 0);
+    test("1^0 mod 2", "1", 0);
+    test("2^1 mod 1", "0", 0);
+
     test("sgn 0", "0", 0);
     test("sgn 3", "1", 0);
     test("sgn −3", "−1", 0);