diff options
Diffstat (limited to 'src/mp.c')
| -rw-r--r-- | src/mp.c | 309 |
1 files changed, 250 insertions, 59 deletions
@@ -37,14 +37,12 @@ mperr(const char *format, ...) mp_error = strdup(text); } - const char * mp_get_error() { return mp_error; } - void mp_clear_error() { if (mp_error) @@ -52,7 +50,6 @@ void mp_clear_error() mp_error = NULL; } - MPNumber mp_new(void) { @@ -61,6 +58,14 @@ mp_new(void) return z; } +MPNumber +mp_new_from_unsigned_integer(ulong x) +{ + MPNumber z; + mpc_init2(z.num, PRECISION); + mpc_set_ui(z.num, x, MPC_RNDNN); + return z; +} MPNumber * mp_new_ptr(void) @@ -70,7 +75,6 @@ mp_new_ptr(void) return z; } - void mp_clear(MPNumber *z) { @@ -78,7 +82,6 @@ mp_clear(MPNumber *z) mpc_clear(z->num); } - void mp_free(MPNumber *z) { @@ -89,7 +92,6 @@ mp_free(MPNumber *z) } } - void mp_get_eulers(MPNumber *z) { @@ -99,14 +101,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 +114,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 +135,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,9 +159,8 @@ 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) +mp_add_integer(const MPNumber *x, long y, MPNumber *z) { mpc_add_si(z->num, x->num, y, MPC_RNDNN); } @@ -196,7 +191,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 +200,6 @@ mp_fractional_part(const MPNumber *x, MPNumber *z) mp_clear(&f); } - void mp_floor(const MPNumber *x, MPNumber *z) { @@ -214,7 +207,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 +214,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) { @@ -250,7 +241,7 @@ mp_divide(const MPNumber *x, const MPNumber *y, MPNumber *z) } void -mp_divide_integer(const MPNumber *x, int64_t y, MPNumber *z) +mp_divide_integer(const MPNumber *x, long y, MPNumber *z) { MPNumber t1 = mp_new(); @@ -268,7 +259,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 +268,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 +277,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 +290,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 +309,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,9 +361,8 @@ 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) +mp_logarithm(long n, const MPNumber *x, MPNumber *z) { /* log(0) undefined */ if (mp_is_zero(x)) @@ -406,39 +388,41 @@ 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) +mp_multiply_integer(const MPNumber *x, long y, MPNumber *z) { - mpc_mul_si(z->num, x->num, (long) y, MPC_RNDNN); + mpc_mul_si(z->num, x->num, 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) { - mpc_set_si(z->num, 1, MPC_RNDNN); - mpc_fr_div(z->num, mpc_realref(z->num), x->num, MPC_RNDNN); + mpc_t temp; + mpc_init2(temp, PRECISION); + mpc_set_ui(temp, 1, MPC_RNDNN); + mpc_fr_div(z->num, mpc_realref(temp), x->num, MPC_RNDNN); + mpc_clear(temp); } void -mp_root(const MPNumber *x, int64_t n, MPNumber *z) +mp_root(const MPNumber *x, long n, MPNumber *z) { - uint64_t p; + ulong p; if (n < 0) { - if (n == INT64_MIN) - p = (uint64_t) INT64_MAX + 1; + mpc_ui_div(z->num, 1, x->num, MPC_RNDNN); + + if (n == LONG_MIN) + p = (ulong) LONG_MAX + 1; else - p = -n; + p = (ulong) -n; } else if (n > 0) { @@ -453,28 +437,26 @@ mp_root(const MPNumber *x, int64_t n, MPNumber *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_rootn_ui(mpc_realref(z->num), mpc_realref(z->num), p, MPFR_RNDN); mpfr_set_zero(mpc_imagref(z->num), MPFR_RNDN); } else { mpfr_t tmp; mpfr_init2(tmp, PRECISION); - mpfr_set_ui(tmp, (uint64_t) p, MPFR_RNDN); + mpfr_set_ui(tmp, 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) { mp_root(x, 2, z); } - void mp_factorial(const MPNumber *x, MPNumber *z) { @@ -508,13 +490,12 @@ mp_factorial(const MPNumber *x, MPNumber *z) else { /* Convert to integer - if couldn't be converted then the factorial would be too big anyway */ - int64_t value = mp_to_integer(x); + ulong value = mp_to_unsigned_integer(x); mpfr_fac_ui(mpc_realref(z->num), value, MPFR_RNDN); mpfr_set_zero(mpc_imagref(z->num), MPFR_RNDN); } } - void mp_modulus_divide(const MPNumber *x, const MPNumber *y, MPNumber *z) { @@ -542,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) { @@ -586,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) { @@ -614,9 +593,8 @@ 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) +mp_xpowy_integer(const MPNumber *x, long n, MPNumber *z) { /* 0^-n invalid */ if (mp_is_zero(x) && n < 0) @@ -626,7 +604,7 @@ mp_xpowy_integer(const MPNumber *x, int64_t n, MPNumber *z) return; } - mpc_pow_si(z->num, x->num, (long) n, MPC_RNDNN); + mpc_pow_si(z->num, x->num, n, MPC_RNDNN); } void @@ -663,6 +641,204 @@ 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, ulong rounds) +{ + MPNumber tmp = mp_new(); + MPNumber two = mp_new_from_unsigned_integer(2); + ulong 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_unsigned_integer(1); + MPNumber a = mp_new_from_unsigned_integer(1); + MPNumber b = mp_new(); + + for (ulong 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 (long j = 1; j < l; j++) + { + 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_unsigned_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, ulong i, MPNumber *z) +{ + MPNumber one = mp_new_from_unsigned_integer(1); + MPNumber two = mp_new_from_unsigned_integer(2); + MPNumber x = mp_new_from_unsigned_integer(i); + MPNumber y = mp_new_from_unsigned_integer(2); + MPNumber d = mp_new_from_unsigned_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(); + ulong 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) { @@ -709,8 +885,7 @@ mp_factorize(const MPNumber *x) return list; } - MPNumber divisor = mp_new(); - mp_set_from_integer(2, &divisor); + MPNumber divisor = mp_new_from_unsigned_integer(2); while (TRUE) { mp_divide(&value, &divisor, &tmp); @@ -726,10 +901,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)) @@ -743,8 +921,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); } } @@ -763,7 +955,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); |