Improve factorization speed

- implements Miller-Rabin primality test
- implements Pollard's rho algorithm to find prime factors
- adds mp_new_from_integer function for convenience

2 files changed, 231 insertions, 39 deletions

diff --git a/src/mp.c b/src/mp.c
index 7b224c6..6563b90 100644
--- a/src/mp.c
+++ b/src/mp.c
@@ -61,6 +61,14 @@ mp_new(void)
     return z;
 }
 
+MPNumber
+mp_new_from_integer(uint64_t x)
+{
+    MPNumber z;
+    mpc_init2(z.num, PRECISION);
+    mpc_set_si(z.num, x, MPC_RNDNN);
+    return z;
+}
 
 MPNumber *
 mp_new_ptr(void)
@@ -70,7 +78,6 @@ mp_new_ptr(void)
     return z;
 }
 
-
 void
 mp_clear(MPNumber *z)
 {
@@ -78,7 +85,6 @@ mp_clear(MPNumber *z)
         mpc_clear(z->num);
 }
 
-
 void
 mp_free(MPNumber *z)
 {
@@ -89,7 +95,6 @@ mp_free(MPNumber *z)
     }
 }
 
-
 void
 mp_get_eulers(MPNumber *z)
 {
@@ -99,14 +104,12 @@ mp_get_eulers(MPNumber *z)
     mpfr_set_zero(mpc_imagref(z->num), 0);
 }
 
-
 void
 mp_get_i(MPNumber *z)
 {
     mpc_set_si_si(z->num, 0, 1, MPC_RNDNN);
 }
 
-
 void
 mp_abs(const MPNumber *x, MPNumber *z)
 {
@@ -114,7 +117,6 @@ mp_abs(const MPNumber *x, MPNumber *z)
     mpc_abs(mpc_realref(z->num), x->num, MPC_RNDNN);
 }
 
-
 void
 mp_arg(const MPNumber *x, MPAngleUnit unit, MPNumber *z)
 {
@@ -136,21 +138,18 @@ mp_arg(const MPNumber *x, MPAngleUnit unit, MPNumber *z)
         mpfr_abs(mpc_realref(z->num), mpc_realref(z->num), MPFR_RNDN);
 }
 
-
 void
 mp_conjugate(const MPNumber *x, MPNumber *z)
 {
     mpc_conj(z->num, x->num, MPC_RNDNN);
 }
 
-
 void
 mp_real_component(const MPNumber *x, MPNumber *z)
 {
     mpc_set_fr(z->num, mpc_realref(x->num), MPC_RNDNN);
 }
 
-
 void
 mp_imaginary_component(const MPNumber *x, MPNumber *z)
 {
@@ -163,7 +162,6 @@ mp_add(const MPNumber *x, const MPNumber *y, MPNumber *z)
     mpc_add(z->num, x->num, y->num, MPC_RNDNN);
 }
 
-
 void
 mp_add_integer(const MPNumber *x, int64_t y, MPNumber *z)
 {
@@ -196,7 +194,6 @@ mp_fractional_component(const MPNumber *x, MPNumber *z)
     mpfr_frac(mpc_realref(z->num), mpc_realref(x->num), MPFR_RNDN);
 }
 
-
 void
 mp_fractional_part(const MPNumber *x, MPNumber *z)
 {
@@ -206,7 +203,6 @@ mp_fractional_part(const MPNumber *x, MPNumber *z)
     mp_clear(&f);
 }
 
-
 void
 mp_floor(const MPNumber *x, MPNumber *z)
 {
@@ -214,7 +210,6 @@ mp_floor(const MPNumber *x, MPNumber *z)
     mpfr_floor(mpc_realref(z->num), mpc_realref(x->num));
 }
 
-
 void
 mp_ceiling(const MPNumber *x, MPNumber *z)
 {
@@ -222,7 +217,6 @@ mp_ceiling(const MPNumber *x, MPNumber *z)
     mpfr_ceil(mpc_realref(z->num), mpc_realref(x->num));
 }
 
-
 void
 mp_round(const MPNumber *x, MPNumber *z)
 {
@@ -268,7 +262,6 @@ mp_is_integer(const MPNumber *x)
     return mpfr_integer_p(mpc_realref(x->num)) != 0;
 }
 
-
 bool
 mp_is_positive_integer(const MPNumber *x)
 {
@@ -278,7 +271,6 @@ mp_is_positive_integer(const MPNumber *x)
         return mpfr_sgn(mpc_realref(x->num)) >= 0 && mp_is_integer(x);
 }
 
-
 bool
 mp_is_natural(const MPNumber *x)
 {
@@ -288,14 +280,12 @@ mp_is_natural(const MPNumber *x)
         return mpfr_sgn(mpc_realref(x->num)) > 0 && mp_is_integer(x);
 }
 
-
 bool
 mp_is_complex(const MPNumber *x)
 {
     return !mpfr_zero_p(mpc_imagref(x->num));
 }
 
-
 bool
 mp_is_equal(const MPNumber *x, const MPNumber *y)
 {
@@ -303,7 +293,6 @@ mp_is_equal(const MPNumber *x, const MPNumber *y)
     return MPC_INEX_RE(res) == 0 && MPC_INEX_IM(res) == 0;
 }
 
-
 void
 mp_epowy(const MPNumber *x, MPNumber *z)
 {
@@ -323,21 +312,18 @@ mp_is_negative(const MPNumber *x)
     return mpfr_sgn(mpc_realref(x->num)) < 0;
 }
 
-
 bool
 mp_is_greater_equal(const MPNumber *x, const MPNumber *y)
 {
     return mp_compare(x, y) >= 0;
 }
 
-
 bool
 mp_is_greater_than(const MPNumber *x, const MPNumber *y)
 {
     return mp_compare(x, y) > 0;
 }
 
-
 bool
 mp_is_less_equal(const MPNumber *x, const MPNumber *y)
 {
@@ -378,7 +364,6 @@ mp_ln(const MPNumber *x, MPNumber *z)
         mpfr_abs(mpc_imagref(z->num), mpc_imagref(z->num), MPFR_RNDN);
 }
 
-
 void
 mp_logarithm(int64_t n, const MPNumber *x, MPNumber *z)
 {
@@ -406,21 +391,18 @@ mp_multiply(const MPNumber *x, const MPNumber *y, MPNumber *z)
     mpc_mul(z->num, x->num, y->num, MPC_RNDNN);
 }
 
-
 void
 mp_multiply_integer(const MPNumber *x, int64_t y, MPNumber *z)
 {
     mpc_mul_si(z->num, x->num, (long) y, MPC_RNDNN);
 }
 
-
 void
 mp_invert_sign(const MPNumber *x, MPNumber *z)
 {
     mpc_neg(z->num, x->num, MPC_RNDNN);
 }
 
-
 void
 mp_reciprocal(const MPNumber *x, MPNumber *z)
 {
@@ -469,14 +451,12 @@ mp_root(const MPNumber *x, int64_t n, MPNumber *z)
     }
 }
 
-
 void
 mp_sqrt(const MPNumber *x, MPNumber *z)
 {
     mp_root(x, 2, z);
 }
 
-
 void
 mp_factorial(const MPNumber *x, MPNumber *z)
 {
@@ -516,7 +496,6 @@ mp_factorial(const MPNumber *x, MPNumber *z)
     }
 }
 
-
 void
 mp_modulus_divide(const MPNumber *x, const MPNumber *y, MPNumber *z)
 {
@@ -544,7 +523,6 @@ mp_modulus_divide(const MPNumber *x, const MPNumber *y, MPNumber *z)
     mp_clear(&t2);
 }
 
-
 void
 mp_modular_exponentiation(const MPNumber *x, const MPNumber *y, const MPNumber *p, MPNumber *z)
 {
@@ -588,7 +566,6 @@ mp_modular_exponentiation(const MPNumber *x, const MPNumber *y, const MPNumber *
     mp_clear(&tmp);
 }
 
-
 void
 mp_xpowy(const MPNumber *x, const MPNumber *y, MPNumber *z)
 {
@@ -616,7 +593,6 @@ mp_xpowy(const MPNumber *x, const MPNumber *y, MPNumber *z)
     mpc_pow(z->num, x->num, y->num, MPC_RNDNN);
 }
 
-
 void
 mp_xpowy_integer(const MPNumber *x, int64_t n, MPNumber *z)
 {
@@ -665,6 +641,205 @@ mp_zeta(const MPNumber *x, MPNumber *z)
     mp_clear(&one);
 }
 
+/***********************************************************************/
+/**                          FACTORIZATION                            **/
+/***********************************************************************/
+
+
+/**
+ * mp_is_pprime uses the Miller-Rabin primality test to decide
+ * whether or not a number is probable prime.
+ * For special values of @n and @rounds it can be deterministic,
+ * but in general the probability of declaring @n as prime if it
+ * is not is at most 4^(-@rounds).
+ * @n has to be odd.
+ * Returns TRUE if @n is probable prime and FALSE otherwise.
+ */
+static bool
+mp_is_pprime(MPNumber *n, uint64_t rounds)
+{
+    MPNumber tmp = mp_new();
+    MPNumber two = mp_new_from_integer(2);
+    uint64_t l = 0;
+    bool is_pprime = TRUE;
+
+    /* Write t := n-1 = 2^l * q with q odd */
+    MPNumber q = mp_new();
+    MPNumber t = mp_new();
+    mp_add_integer(n, -1, &t);
+    mp_set_from_mp(&t, &q);
+    do
+    {
+        mp_divide_integer(&q, 2, &q);
+        mp_modulus_divide(&q, &two, &tmp);
+        l++;
+    } while (mp_is_zero(&tmp));
+
+    /* @rounds Miller-Rabin tests to bases a = 2,3,...,@rounds+1 */
+    MPNumber one = mp_new_from_integer(1);
+    MPNumber a = mp_new_from_integer(1);
+    MPNumber b = mp_new();
+
+    for (uint64_t i = 1; (i < mp_to_integer(&t)) && (i <= rounds+1); i++)
+    {
+        mp_add_integer(&a, 1, &a);
+        mp_modular_exponentiation(&a, &q, n, &b);
+        if (mp_compare(&one, &b) == 0 || mp_compare(&t, &b) == 0)
+        {
+            continue;
+        }
+
+        bool is_witness = FALSE;
+        for (int i = 1; i < l; i++)
+        {
+            mp_modular_exponentiation(&b, &two, n, &b);
+            if (mp_compare(&b, &t) == 0)
+            {
+                is_witness = TRUE;
+                break;
+            }
+        }
+
+        if (!is_witness)
+        {
+            is_pprime = FALSE;
+            break;
+        }
+    }
+
+    mp_clear(&t);
+    mp_clear(&q);
+    mp_clear(&a);
+    mp_clear(&b);
+    mp_clear(&one);
+    mp_clear(&two);
+    mp_clear(&tmp);
+
+    return is_pprime;
+}
+
+/**
+ * Sets z = gcd(a,b) where gcd(a,b) is the greatest common divisor of @a and @b.
+ */
+static void
+mp_gcd (const MPNumber *a, const MPNumber *b, MPNumber *z)
+{
+    MPNumber null = mp_new_from_integer(0);
+    MPNumber t1 = mp_new();
+    MPNumber t2 = mp_new();
+
+    mp_set_from_mp(a, z);
+    mp_set_from_mp(b, &t2);
+
+    while (mp_compare(&t2, &null) != 0)
+    {
+        mp_set_from_mp(&t2, &t1);
+        mp_modulus_divide(z, &t2, &t2);
+        mp_set_from_mp(&t1, z);
+    }
+
+    mp_clear(&null);
+    mp_clear(&t1);
+    mp_clear(&t2);
+}
+
+/**
+ * mp_pollard_rho searches a divisor of @n using Pollard's rho algorithm.
+ * @i is the start value of the pseudorandom sequence which is generated
+ * by the polynomial x^2+1 mod n.
+ *
+ * Returns TRUE if a divisor was found and stores the divisor in z.
+ * Returns FALSE otherwise.
+ */
+static bool
+mp_pollard_rho (const MPNumber *n, uint64_t i, MPNumber *z)
+{
+    MPNumber one = mp_new_from_integer(1);
+    MPNumber two = mp_new_from_integer(2);
+    MPNumber x = mp_new_from_integer(i);
+    MPNumber y = mp_new_from_integer(2);
+    MPNumber d = mp_new_from_integer(1);
+
+    while (mp_compare(&d, &one) == 0)
+    {
+        mp_modular_exponentiation(&x, &two, n, &x);
+        mp_add(&x, &one, &x);
+
+        mp_modular_exponentiation(&y, &two, n, &y);
+        mp_add(&y, &one, &y);
+
+        mp_modular_exponentiation(&y, &two, n, &y);
+        mp_add(&y, &one, &y);
+
+        mp_subtract(&x, &y,z);
+        mp_abs(z, z);
+        mp_gcd(z, n, &d);
+    }
+
+    if (mp_compare(&d, n) == 0)
+    {
+        mp_clear(&one);
+        mp_clear(&two);
+        mp_clear(&x);
+        mp_clear(&y);
+        mp_clear(&d);
+
+        return FALSE;
+    }
+    else
+    {
+        mp_set_from_mp(&d, z);
+
+        mp_clear(&one);
+        mp_clear(&two);
+        mp_clear(&x);
+        mp_clear(&y);
+        mp_clear(&d);
+
+        return TRUE;
+    }
+}
+
+/**
+ * find_big_prime_factor acts as driver function for mp_pollard_rho which
+ * is run as long as a prime factor is found.
+ * On success sets @z to a prime factor of @n.
+ */
+static void
+find_big_prime_factor (const MPNumber *n, MPNumber *z)
+{
+    MPNumber tmp = mp_new();
+    uint64_t i = 2;
+
+    while (TRUE)
+    {
+        while (mp_pollard_rho (n, i, &tmp) == FALSE)
+        {
+            i++;
+        }
+
+        if (!mp_is_pprime(&tmp, 50))
+        {
+            mp_divide(n, &tmp, &tmp);
+        }
+        else
+            break;
+    }
+
+    mp_set_from_mp(&tmp, z);
+    mp_clear(&tmp);
+}
+
+/**
+ * mp_factorize tries to factorize the value of @x.
+ * If @x < 2^64 it calls mp_factorize_unit64 which deals in integers
+ * and should be fast enough for most values.
+ * If @x > 2^64 the approach to find factors of @x is as follows:
+ *   - Try to divide @x by prime numbers 2,3,5,7,.. up to min(2^13, sqrt(x))
+ *   - Use Pollard rho to find prime factors > 2^13.
+ * Returns a pointer to a GList with all prime factors of @x which needs to
+ * be freed.
+ */
 GList*
 mp_factorize(const MPNumber *x)
 {
@@ -711,8 +886,7 @@ mp_factorize(const MPNumber *x)
         return list;
     }
 
-    MPNumber divisor = mp_new();
-    mp_set_from_integer(2, &divisor);
+    MPNumber divisor = mp_new_from_integer(2);
     while (TRUE)
     {
         mp_divide(&value, &divisor, &tmp);
@@ -728,10 +902,13 @@ mp_factorize(const MPNumber *x)
             break;
     }
 
-    MPNumber root = mp_new();
     mp_set_from_integer(3, &divisor);
+
+    MPNumber root = mp_new();
     mp_sqrt(&value, &root);
-    while (mp_is_less_equal(&divisor, &root))
+    uint64_t max_trial_division = (uint64_t) (1 << 10);
+    uint64_t iter = 0;
+    while (mp_is_less_equal(&divisor, &root) && (iter++ < max_trial_division))
     {
         mp_divide(&value, &divisor, &tmp);
         if (mp_is_integer(&tmp))
@@ -745,8 +922,22 @@ mp_factorize(const MPNumber *x)
         }
         else
         {
-            mp_add_integer(&divisor, 2, &tmp);
-            mp_set_from_mp(&tmp, &divisor);
+            mp_add_integer(&divisor, 2, &divisor);
+        }
+    }
+
+    while (!mp_is_pprime(&value, 50))
+    {
+        find_big_prime_factor (&value, &divisor);
+
+        mp_divide(&value, &divisor, &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));
+            mpc_init2(factor->num, PRECISION);
         }
     }
 
@@ -765,7 +956,6 @@ mp_factorize(const MPNumber *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);
diff --git a/src/mp.h b/src/mp.h
index 3e03089..8c527e3 100644
--- a/src/mp.h
+++ b/src/mp.h
@@ -69,6 +69,8 @@ void        mperr(const char *format, ...) __attribute__((format(printf, 1, 2)))
 /* Returns initialized MPNumber object */
 MPNumber    mp_new(void);
 
+MPNumber    mp_new_from_integer(uint64_t x);
+
 MPNumber*   mp_new_ptr(void);
 
 void        mp_clear(MPNumber *z);