1 /*
2 * Copyright (C) 1987-2008 Sun Microsystems, Inc. All Rights Reserved.
3 * Copyright (C) 2008-2011 Robert Ancell
4 *
5 * This program is free software: you can redistribute it and/or modify it under
6 * the terms of the GNU General Public License as published by the Free Software
7 * Foundation, either version 2 of the License, or (at your option) any later
8 * version. See http://www.gnu.org/copyleft/gpl.html the full text of the
9 * license.
10 */
11
12 #include <stdlib.h>
13 #include <stdio.h>
14 #include <stdint.h>
15 #include <math.h>
16 #include <errno.h>
17
18 #include "mp.h"
19
20 char *mp_error = NULL;
21
22 /* THIS ROUTINE IS CALLED WHEN AN ERROR CONDITION IS ENCOUNTERED, AND
23 * AFTER A MESSAGE HAS BEEN WRITTEN TO STDERR.
24 */
25 void
mperr(const char * format,...)26 mperr(const char *format, ...)
27 {
28 char text[1024];
29 va_list args;
30
31 va_start(args, format);
32 vsnprintf(text, 1024, format, args);
33 va_end(args);
34
35 if (mp_error)
36 free(mp_error);
37 mp_error = strdup(text);
38 }
39
40 const char *
mp_get_error()41 mp_get_error()
42 {
43 return mp_error;
44 }
45
mp_clear_error()46 void mp_clear_error()
47 {
48 if (mp_error)
49 free(mp_error);
50 mp_error = NULL;
51 }
52
53 MPNumber
mp_new(void)54 mp_new(void)
55 {
56 MPNumber z;
57 mpc_init2(z.num, PRECISION);
58 return z;
59 }
60
61 MPNumber
mp_new_from_unsigned_integer(unsigned long x)62 mp_new_from_unsigned_integer(unsigned long x)
63 {
64 MPNumber z;
65 mpc_init2(z.num, PRECISION);
66 mpc_set_ui(z.num, x, MPC_RNDNN);
67 return z;
68 }
69
70 MPNumber *
mp_new_ptr(void)71 mp_new_ptr(void)
72 {
73 MPNumber *z = malloc(sizeof(MPNumber));
74 mpc_init2(z->num, PRECISION);
75 return z;
76 }
77
78 void
mp_clear(MPNumber * z)79 mp_clear(MPNumber *z)
80 {
81 if (z != NULL)
82 mpc_clear(z->num);
83 }
84
85 void
mp_free(MPNumber * z)86 mp_free(MPNumber *z)
87 {
88 if (z != NULL)
89 {
90 mpc_clear(z->num);
91 free(z);
92 }
93 }
94
95 void
mp_get_eulers(MPNumber * z)96 mp_get_eulers(MPNumber *z)
97 {
98 /* e^1, since mpfr doesn't have a function to return e */
99 mpfr_set_ui(mpc_realref(z->num), 1, MPFR_RNDN);
100 mpfr_exp(mpc_realref(z->num), mpc_realref(z->num), MPFR_RNDN);
101 mpfr_set_zero(mpc_imagref(z->num), 0);
102 }
103
104 void
mp_get_i(MPNumber * z)105 mp_get_i(MPNumber *z)
106 {
107 mpc_set_si_si(z->num, 0, 1, MPC_RNDNN);
108 }
109
110 void
mp_abs(const MPNumber * x,MPNumber * z)111 mp_abs(const MPNumber *x, MPNumber *z)
112 {
113 mpfr_set_zero(mpc_imagref(z->num), MPFR_RNDN);
114 mpc_abs(mpc_realref(z->num), x->num, MPC_RNDNN);
115 }
116
117 void
mp_arg(const MPNumber * x,MPAngleUnit unit,MPNumber * z)118 mp_arg(const MPNumber *x, MPAngleUnit unit, MPNumber *z)
119 {
120 if (mp_is_zero(x))
121 {
122 /* Translators: Error display when attempting to take argument of zero */
123 mperr(_("Argument not defined for zero"));
124 mp_set_from_integer(0, z);
125 return;
126 }
127
128 mpfr_set_zero(mpc_imagref(z->num), MPFR_RNDN);
129 mpc_arg(mpc_realref(z->num), x->num, MPC_RNDNN);
130 convert_from_radians(z, unit, z);
131 // MPC returns -π for the argument of negative real numbers if
132 // their imaginary part is -0, we want +π for all real negative
133 // numbers
134 if (!mp_is_complex (x) && mp_is_negative (x))
135 mpfr_abs(mpc_realref(z->num), mpc_realref(z->num), MPFR_RNDN);
136 }
137
138 void
mp_conjugate(const MPNumber * x,MPNumber * z)139 mp_conjugate(const MPNumber *x, MPNumber *z)
140 {
141 mpc_conj(z->num, x->num, MPC_RNDNN);
142 }
143
144 void
mp_real_component(const MPNumber * x,MPNumber * z)145 mp_real_component(const MPNumber *x, MPNumber *z)
146 {
147 mpc_set_fr(z->num, mpc_realref(x->num), MPC_RNDNN);
148 }
149
150 void
mp_imaginary_component(const MPNumber * x,MPNumber * z)151 mp_imaginary_component(const MPNumber *x, MPNumber *z)
152 {
153 mpc_set_fr(z->num, mpc_imagref(x->num), MPC_RNDNN);
154 }
155
156 void
mp_add(const MPNumber * x,const MPNumber * y,MPNumber * z)157 mp_add(const MPNumber *x, const MPNumber *y, MPNumber *z)
158 {
159 mpc_add(z->num, x->num, y->num, MPC_RNDNN);
160 }
161
162 void
mp_add_integer(const MPNumber * x,long y,MPNumber * z)163 mp_add_integer(const MPNumber *x, long y, MPNumber *z)
164 {
165 mpc_add_si(z->num, x->num, y, MPC_RNDNN);
166 }
167
168 void
mp_subtract(const MPNumber * x,const MPNumber * y,MPNumber * z)169 mp_subtract(const MPNumber *x, const MPNumber *y, MPNumber *z)
170 {
171 mpc_sub(z->num, x->num, y->num, MPC_RNDNN);
172 }
173
174 void
mp_sgn(const MPNumber * x,MPNumber * z)175 mp_sgn(const MPNumber *x, MPNumber *z)
176 {
177 mpc_set_si(z->num, mpfr_sgn(mpc_realref(x->num)), MPC_RNDNN);
178 }
179
180 void
mp_integer_component(const MPNumber * x,MPNumber * z)181 mp_integer_component(const MPNumber *x, MPNumber *z)
182 {
183 mpfr_set_zero(mpc_imagref(z->num), MPFR_RNDN);
184 mpfr_trunc(mpc_realref(z->num), mpc_realref(x->num));
185 }
186
187 void
mp_fractional_component(const MPNumber * x,MPNumber * z)188 mp_fractional_component(const MPNumber *x, MPNumber *z)
189 {
190 mpfr_set_zero(mpc_imagref(z->num), MPFR_RNDN);
191 mpfr_frac(mpc_realref(z->num), mpc_realref(x->num), MPFR_RNDN);
192 }
193
194 void
mp_fractional_part(const MPNumber * x,MPNumber * z)195 mp_fractional_part(const MPNumber *x, MPNumber *z)
196 {
197 MPNumber f = mp_new();
198 mp_floor(x, &f);
199 mp_subtract(x, &f, z);
200 mp_clear(&f);
201 }
202
203 void
mp_floor(const MPNumber * x,MPNumber * z)204 mp_floor(const MPNumber *x, MPNumber *z)
205 {
206 mpfr_set_zero(mpc_imagref(z->num), MPFR_RNDN);
207 mpfr_floor(mpc_realref(z->num), mpc_realref(x->num));
208 }
209
210 void
mp_ceiling(const MPNumber * x,MPNumber * z)211 mp_ceiling(const MPNumber *x, MPNumber *z)
212 {
213 mpfr_set_zero(mpc_imagref(z->num), MPFR_RNDN);
214 mpfr_ceil(mpc_realref(z->num), mpc_realref(x->num));
215 }
216
217 void
mp_round(const MPNumber * x,MPNumber * z)218 mp_round(const MPNumber *x, MPNumber *z)
219 {
220 mpfr_set_zero(mpc_imagref(z->num), MPFR_RNDN);
221 mpfr_round(mpc_realref(z->num), mpc_realref(x->num));
222 }
223
224 int
mp_compare(const MPNumber * x,const MPNumber * y)225 mp_compare(const MPNumber *x, const MPNumber *y)
226 {
227 return mpfr_cmp(mpc_realref(x->num), mpc_realref(y->num));
228 }
229
230 void
mp_divide(const MPNumber * x,const MPNumber * y,MPNumber * z)231 mp_divide(const MPNumber *x, const MPNumber *y, MPNumber *z)
232 {
233 if (mp_is_zero(y))
234 {
235 /* Translators: Error displayed attempted to divide by zero */
236 mperr(_("Division by zero is undefined"));
237 mp_set_from_integer(0, z);
238 return;
239 }
240 mpc_div(z->num, x->num, y->num, MPC_RNDNN);
241 }
242
243 void
mp_divide_integer(const MPNumber * x,long y,MPNumber * z)244 mp_divide_integer(const MPNumber *x, long y, MPNumber *z)
245 {
246 MPNumber t1 = mp_new();
247
248 mp_set_from_integer(y, &t1);
249 mp_divide(x, &t1, z);
250 mp_clear(&t1);
251 }
252
253 bool
mp_is_integer(const MPNumber * x)254 mp_is_integer(const MPNumber *x)
255 {
256 if (mp_is_complex(x))
257 return false;
258
259 return mpfr_integer_p(mpc_realref(x->num)) != 0;
260 }
261
262 bool
mp_is_positive_integer(const MPNumber * x)263 mp_is_positive_integer(const MPNumber *x)
264 {
265 if (mp_is_complex(x))
266 return false;
267 else
268 return mpfr_sgn(mpc_realref(x->num)) >= 0 && mp_is_integer(x);
269 }
270
271 bool
mp_is_natural(const MPNumber * x)272 mp_is_natural(const MPNumber *x)
273 {
274 if (mp_is_complex(x))
275 return false;
276 else
277 return mpfr_sgn(mpc_realref(x->num)) > 0 && mp_is_integer(x);
278 }
279
280 bool
mp_is_complex(const MPNumber * x)281 mp_is_complex(const MPNumber *x)
282 {
283 return !mpfr_zero_p(mpc_imagref(x->num));
284 }
285
286 bool
mp_is_equal(const MPNumber * x,const MPNumber * y)287 mp_is_equal(const MPNumber *x, const MPNumber *y)
288 {
289 int res = mpc_cmp(x->num, y->num);
290 return MPC_INEX_RE(res) == 0 && MPC_INEX_IM(res) == 0;
291 }
292
293 void
mp_epowy(const MPNumber * x,MPNumber * z)294 mp_epowy(const MPNumber *x, MPNumber *z)
295 {
296 mpc_exp(z->num, x->num, MPC_RNDNN);
297 }
298
299 bool
mp_is_zero(const MPNumber * x)300 mp_is_zero (const MPNumber *x)
301 {
302 int res = mpc_cmp_si_si(x->num, 0, 0);
303 return MPC_INEX_RE(res) == 0 && MPC_INEX_IM(res) == 0;
304 }
305
306 bool
mp_is_negative(const MPNumber * x)307 mp_is_negative(const MPNumber *x)
308 {
309 return mpfr_sgn(mpc_realref(x->num)) < 0;
310 }
311
312 bool
mp_is_greater_equal(const MPNumber * x,const MPNumber * y)313 mp_is_greater_equal(const MPNumber *x, const MPNumber *y)
314 {
315 return mp_compare(x, y) >= 0;
316 }
317
318 bool
mp_is_greater_than(const MPNumber * x,const MPNumber * y)319 mp_is_greater_than(const MPNumber *x, const MPNumber *y)
320 {
321 return mp_compare(x, y) > 0;
322 }
323
324 bool
mp_is_less_equal(const MPNumber * x,const MPNumber * y)325 mp_is_less_equal(const MPNumber *x, const MPNumber *y)
326 {
327 return mp_compare(x, y) <= 0;
328 }
329
330 bool
mp_is_less_than(const MPNumber * x,const MPNumber * y)331 mp_is_less_than(const MPNumber *x, const MPNumber *y)
332 {
333 return mp_compare(x, y) < 0;
334 }
335
336 void
mp_ln(const MPNumber * x,MPNumber * z)337 mp_ln(const MPNumber *x, MPNumber *z)
338 {
339 /* ln(0) undefined */
340 if (mp_is_zero(x))
341 {
342 /* Translators: Error displayed when attempting to take logarithm of zero */
343 mperr(_("Logarithm of zero is undefined"));
344 mp_set_from_integer(0, z);
345 return;
346 }
347
348 /* ln(-x) complex */
349 /* FIXME: Make complex numbers optional */
350 /*if (mp_is_negative(x)) {
351 // Translators: Error displayed attempted to take logarithm of negative value
352 mperr(_("Logarithm of negative values is undefined"));
353 mp_set_from_integer(0, z);
354 return;
355 }*/
356
357 mpc_log(z->num, x->num, MPC_RNDNN);
358 // MPC returns -π for the imaginary part of the log of
359 // negative real numbers if their imaginary part is -0, we want +π
360 if (!mp_is_complex (x) && mp_is_negative (x))
361 mpfr_abs(mpc_imagref(z->num), mpc_imagref(z->num), MPFR_RNDN);
362 }
363
364 void
mp_logarithm(long n,const MPNumber * x,MPNumber * z)365 mp_logarithm(long n, const MPNumber *x, MPNumber *z)
366 {
367 /* log(0) undefined */
368 if (mp_is_zero(x))
369 {
370 /* Translators: Error displayed when attempting to take logarithm of zero */
371 mperr(_("Logarithm of zero is undefined"));
372 mp_set_from_integer(0, z);
373 return;
374 }
375
376 /* logn(x) = ln(x) / ln(n) */
377 MPNumber t1 = mp_new();
378 mp_set_from_integer(n, &t1);
379 mp_ln(&t1, &t1);
380 mp_ln(x, z);
381 mp_divide(z, &t1, z);
382 mp_clear(&t1);
383 }
384
385 void
mp_multiply(const MPNumber * x,const MPNumber * y,MPNumber * z)386 mp_multiply(const MPNumber *x, const MPNumber *y, MPNumber *z)
387 {
388 mpc_mul(z->num, x->num, y->num, MPC_RNDNN);
389 }
390
391 void
mp_multiply_integer(const MPNumber * x,long y,MPNumber * z)392 mp_multiply_integer(const MPNumber *x, long y, MPNumber *z)
393 {
394 mpc_mul_si(z->num, x->num, y, MPC_RNDNN);
395 }
396
397 void
mp_invert_sign(const MPNumber * x,MPNumber * z)398 mp_invert_sign(const MPNumber *x, MPNumber *z)
399 {
400 mpc_neg(z->num, x->num, MPC_RNDNN);
401 }
402
403 void
mp_reciprocal(const MPNumber * x,MPNumber * z)404 mp_reciprocal(const MPNumber *x, MPNumber *z)
405 {
406 mpc_t temp;
407 mpc_init2(temp, PRECISION);
408 mpc_set_ui(temp, 1, MPC_RNDNN);
409 mpc_fr_div(z->num, mpc_realref(temp), x->num, MPC_RNDNN);
410 mpc_clear(temp);
411 }
412
413 void
mp_root(const MPNumber * x,long n,MPNumber * z)414 mp_root(const MPNumber *x, long n, MPNumber *z)
415 {
416 unsigned long p;
417
418 if (n < 0)
419 {
420 mpc_ui_div(z->num, 1, x->num, MPC_RNDNN);
421
422 if (n == LONG_MIN)
423 p = (unsigned long) LONG_MAX + 1;
424 else
425 p = (unsigned long) -n;
426 }
427 else if (n > 0)
428 {
429 mp_set_from_mp(x, z);
430 p = n;
431 }
432 else
433 { /* Translators: Error displayed when attempting to take zeroth root */
434 mperr(_("The zeroth root of a number is undefined"));
435 mp_set_from_integer(0, z);
436 return;
437 }
438 if (!mp_is_complex(x) && (!mp_is_negative(x) || (p & 1) == 1))
439 {
440 mpfr_rootn_ui(mpc_realref(z->num), mpc_realref(z->num), p, MPFR_RNDN);
441 mpfr_set_zero(mpc_imagref(z->num), MPFR_RNDN);
442 }
443 else
444 {
445 mpfr_t tmp;
446 mpfr_init2(tmp, PRECISION);
447 mpfr_set_ui(tmp, p, MPFR_RNDN);
448 mpfr_ui_div(tmp, 1, tmp, MPFR_RNDN);
449 mpc_pow_fr(z->num, z->num, tmp, MPC_RNDNN);
450 mpfr_clear(tmp);
451 }
452 }
453
454 void
mp_sqrt(const MPNumber * x,MPNumber * z)455 mp_sqrt(const MPNumber *x, MPNumber *z)
456 {
457 mp_root(x, 2, z);
458 }
459
460 void
mp_factorial(const MPNumber * x,MPNumber * z)461 mp_factorial(const MPNumber *x, MPNumber *z)
462 {
463 /* 0! == 1 */
464 if (mp_is_zero(x))
465 {
466 mpc_set_si(z->num, 1, MPC_RNDNN);
467 return;
468 }
469 if (!mp_is_natural(x))
470 {
471 /* Factorial Not defined for Complex or for negative numbers */
472 if (mp_is_negative(x) || mp_is_complex(x))
473 { /* Translators: Error displayed when attempted take the factorial of a negative or complex number */
474 mperr(_("Factorial is only defined for non-negative real numbers"));
475 mp_set_from_integer(0, z);
476 return;
477 }
478 MPNumber tmp = mp_new();
479 mpfr_t tmp2;
480 mpfr_init2(tmp2, PRECISION);
481 mp_set_from_integer(1, &tmp);
482 mp_add(&tmp, x, &tmp);
483
484 /* Factorial(x) = Gamma(x+1) - This is the formula used to calculate Factorial of positive real numbers.*/
485 mpfr_gamma(tmp2, mpc_realref(tmp.num), MPFR_RNDN);
486 mpc_set_fr(z->num, tmp2, MPC_RNDNN);
487 mp_clear(&tmp);
488 mpfr_clear(tmp2);
489 }
490 else
491 {
492 /* Convert to integer - if couldn't be converted then the factorial would be too big anyway */
493 unsigned long value = mp_to_unsigned_integer(x);
494 mpfr_fac_ui(mpc_realref(z->num), value, MPFR_RNDN);
495 mpfr_set_zero(mpc_imagref(z->num), MPFR_RNDN);
496 }
497 }
498
499 void
mp_modulus_divide(const MPNumber * x,const MPNumber * y,MPNumber * z)500 mp_modulus_divide(const MPNumber *x, const MPNumber *y, MPNumber *z)
501 {
502 if (!mp_is_integer(x) || !mp_is_integer(y))
503 { /* Translators: Error displayed when attemping to do a modulus division on non-integer numbers */
504 mperr(_("Modulus division is only defined for integers"));
505 mp_set_from_integer(0, z);
506 return;
507 }
508
509 MPNumber t1 = mp_new();
510 MPNumber t2 = mp_new();
511
512 mp_divide(x, y, &t1);
513 mp_floor(&t1, &t1);
514 mp_multiply(&t1, y, &t2);
515 mp_subtract(x, &t2, z);
516
517 mp_set_from_integer(0, &t1);
518 if ((mp_compare(y, &t1) < 0 && mp_compare(z, &t1) > 0) ||
519 (mp_compare(y, &t1) > 0 && mp_compare(z, &t1) < 0))
520 mp_add(z, y, z);
521
522 mp_clear(&t1);
523 mp_clear(&t2);
524 }
525
526 void
mp_modular_exponentiation(const MPNumber * x,const MPNumber * y,const MPNumber * p,MPNumber * z)527 mp_modular_exponentiation(const MPNumber *x, const MPNumber *y, const MPNumber *p, MPNumber *z)
528 {
529 MPNumber base_value = mp_new();
530 MPNumber exp_value = mp_new();
531 MPNumber ans = mp_new();
532 MPNumber two = mp_new();
533 MPNumber tmp = mp_new();
534
535 mp_set_from_integer(1, &ans);
536 mp_set_from_integer(2, &two);
537 mp_abs(y, &exp_value);
538
539 if (mp_is_negative(y))
540 mp_reciprocal(x, &base_value);
541 else
542 mp_set_from_mp(x, &base_value);
543
544 while (!mp_is_zero(&exp_value))
545 {
546 mp_modulus_divide(&exp_value, &two, &tmp);
547
548 bool is_even = mp_is_zero(&tmp);
549 if (!is_even)
550 {
551 mp_multiply(&ans, &base_value, &ans);
552 mp_modulus_divide(&ans, p, &ans);
553 }
554 mp_multiply(&base_value, &base_value, &base_value);
555 mp_modulus_divide(&base_value, p, &base_value);
556 mp_divide_integer(&exp_value, 2, &exp_value);
557 mp_floor(&exp_value, &exp_value);
558 }
559
560 mp_modulus_divide(&ans, p, z);
561
562 mp_clear(&base_value);
563 mp_clear(&exp_value);
564 mp_clear(&ans);
565 mp_clear(&two);
566 mp_clear(&tmp);
567 }
568
569 void
mp_xpowy(const MPNumber * x,const MPNumber * y,MPNumber * z)570 mp_xpowy(const MPNumber *x, const MPNumber *y, MPNumber *z)
571 {
572 /* 0^-n invalid */
573 if (mp_is_zero(x) && mp_is_negative(y))
574 { /* Translators: Error displayed when attempted to raise 0 to a negative exponent */
575 mperr(_("The power of zero is undefined for a negative exponent"));
576 mp_set_from_integer(0, z);
577 return;
578 }
579
580 if (!mp_is_complex(x) && !mp_is_complex(y) && !mp_is_integer(y))
581 {
582 MPNumber reciprocal = mp_new();
583 mp_reciprocal(y, &reciprocal);
584 if (mp_is_integer(&reciprocal))
585 {
586 mp_root(x, mp_to_integer(&reciprocal), z);
587 mp_clear(&reciprocal);
588 return;
589 }
590 mp_clear(&reciprocal);
591 }
592
593 mpc_pow(z->num, x->num, y->num, MPC_RNDNN);
594 }
595
596 void
mp_xpowy_integer(const MPNumber * x,long n,MPNumber * z)597 mp_xpowy_integer(const MPNumber *x, long n, MPNumber *z)
598 {
599 /* 0^-n invalid */
600 if (mp_is_zero(x) && n < 0)
601 { /* Translators: Error displayed when attempted to raise 0 to a negative re_exponent */
602 mperr(_("The power of zero is undefined for a negative exponent"));
603 mp_set_from_integer(0, z);
604 return;
605 }
606
607 mpc_pow_si(z->num, x->num, n, MPC_RNDNN);
608 }
609
610 void
mp_erf(const MPNumber * x,MPNumber * z)611 mp_erf(const MPNumber *x, MPNumber *z)
612 {
613 if (mp_is_complex(x))
614 { /* Translators: Error displayed when error function (erf) value is undefined */
615 mperr(_("The error function is only defined for real numbers"));
616 mp_set_from_integer(0, z);
617 return;
618 }
619
620 mpfr_set_zero(mpc_imagref(z->num), MPFR_RNDN);
621 mpfr_erf(mpc_realref(z->num), mpc_realref(x->num), MPFR_RNDN);
622 }
623
624 void
mp_zeta(const MPNumber * x,MPNumber * z)625 mp_zeta(const MPNumber *x, MPNumber *z)
626 {
627 MPNumber one = mp_new();
628
629 mp_set_from_integer(1, &one);
630 if (mp_is_complex(x) || mp_compare(x, &one) == 0)
631 { /* Translators: Error displayed when zeta function value is undefined */
632 mperr(_("The Riemann zeta function is only defined for real numbers ≠1"));
633 mp_set_from_integer(0, z);
634 mp_clear(&one);
635 return;
636 }
637
638 mpfr_set_zero(mpc_imagref(z->num), MPFR_RNDN);
639 mpfr_zeta(mpc_realref(z->num), mpc_realref(x->num), MPFR_RNDN);
640
641 mp_clear(&one);
642 }
643
644 /***********************************************************************/
645 /** FACTORIZATION **/
646 /***********************************************************************/
647
648
649 /**
650 * mp_is_pprime uses the Miller-Rabin primality test to decide
651 * whether or not a number is probable prime.
652 * For special values of @n and @rounds it can be deterministic,
653 * but in general the probability of declaring @n as prime if it
654 * is not is at most 4^(-@rounds).
655 * @n has to be odd.
656 * Returns TRUE if @n is probable prime and FALSE otherwise.
657 */
658 static bool
mp_is_pprime(MPNumber * n,unsigned long rounds)659 mp_is_pprime(MPNumber *n, unsigned long rounds)
660 {
661 MPNumber tmp = mp_new();
662 MPNumber two = mp_new_from_unsigned_integer(2);
663 unsigned long l = 0;
664 bool is_pprime = TRUE;
665
666 /* Write t := n-1 = 2^l * q with q odd */
667 MPNumber q = mp_new();
668 MPNumber t = mp_new();
669 mp_add_integer(n, -1, &t);
670 mp_set_from_mp(&t, &q);
671 do
672 {
673 mp_divide_integer(&q, 2, &q);
674 mp_modulus_divide(&q, &two, &tmp);
675 l++;
676 } while (mp_is_zero(&tmp));
677
678 /* @rounds Miller-Rabin tests to bases a = 2,3,...,@rounds+1 */
679 MPNumber one = mp_new_from_unsigned_integer(1);
680 MPNumber a = mp_new_from_unsigned_integer(1);
681 MPNumber b = mp_new();
682
683 for (unsigned long i = 1; (i < mp_to_integer(&t)) && (i <= rounds+1); i++)
684 {
685 mp_add_integer(&a, 1, &a);
686 mp_modular_exponentiation(&a, &q, n, &b);
687 if (mp_compare(&one, &b) == 0 || mp_compare(&t, &b) == 0)
688 {
689 continue;
690 }
691
692 bool is_witness = FALSE;
693 for (long j = 1; j < l; j++)
694 {
695 mp_modular_exponentiation(&b, &two, n, &b);
696 if (mp_compare(&b, &t) == 0)
697 {
698 is_witness = TRUE;
699 break;
700 }
701 }
702
703 if (!is_witness)
704 {
705 is_pprime = FALSE;
706 break;
707 }
708 }
709
710 mp_clear(&t);
711 mp_clear(&q);
712 mp_clear(&a);
713 mp_clear(&b);
714 mp_clear(&one);
715 mp_clear(&two);
716 mp_clear(&tmp);
717
718 return is_pprime;
719 }
720
721 /**
722 * Sets z = gcd(a,b) where gcd(a,b) is the greatest common divisor of @a and @b.
723 */
724 static void
mp_gcd(const MPNumber * a,const MPNumber * b,MPNumber * z)725 mp_gcd (const MPNumber *a, const MPNumber *b, MPNumber *z)
726 {
727 MPNumber null = mp_new_from_unsigned_integer(0);
728 MPNumber t1 = mp_new();
729 MPNumber t2 = mp_new();
730
731 mp_set_from_mp(a, z);
732 mp_set_from_mp(b, &t2);
733
734 while (mp_compare(&t2, &null) != 0)
735 {
736 mp_set_from_mp(&t2, &t1);
737 mp_modulus_divide(z, &t2, &t2);
738 mp_set_from_mp(&t1, z);
739 }
740
741 mp_clear(&null);
742 mp_clear(&t1);
743 mp_clear(&t2);
744 }
745
746 /**
747 * mp_pollard_rho searches a divisor of @n using Pollard's rho algorithm.
748 * @i is the start value of the pseudorandom sequence which is generated
749 * by the polynomial x^2+1 mod n.
750 *
751 * Returns TRUE if a divisor was found and stores the divisor in z.
752 * Returns FALSE otherwise.
753 */
754 static bool
mp_pollard_rho(const MPNumber * n,unsigned long i,MPNumber * z)755 mp_pollard_rho (const MPNumber *n, unsigned long i, MPNumber *z)
756 {
757 MPNumber one = mp_new_from_unsigned_integer(1);
758 MPNumber two = mp_new_from_unsigned_integer(2);
759 MPNumber x = mp_new_from_unsigned_integer(i);
760 MPNumber y = mp_new_from_unsigned_integer(2);
761 MPNumber d = mp_new_from_unsigned_integer(1);
762
763 while (mp_compare(&d, &one) == 0)
764 {
765 mp_modular_exponentiation(&x, &two, n, &x);
766 mp_add(&x, &one, &x);
767
768 mp_modular_exponentiation(&y, &two, n, &y);
769 mp_add(&y, &one, &y);
770
771 mp_modular_exponentiation(&y, &two, n, &y);
772 mp_add(&y, &one, &y);
773
774 mp_subtract(&x, &y,z);
775 mp_abs(z, z);
776 mp_gcd(z, n, &d);
777 }
778
779 if (mp_compare(&d, n) == 0)
780 {
781 mp_clear(&one);
782 mp_clear(&two);
783 mp_clear(&x);
784 mp_clear(&y);
785 mp_clear(&d);
786
787 return FALSE;
788 }
789 else
790 {
791 mp_set_from_mp(&d, z);
792
793 mp_clear(&one);
794 mp_clear(&two);
795 mp_clear(&x);
796 mp_clear(&y);
797 mp_clear(&d);
798
799 return TRUE;
800 }
801 }
802
803 /**
804 * find_big_prime_factor acts as driver function for mp_pollard_rho which
805 * is run as long as a prime factor is found.
806 * On success sets @z to a prime factor of @n.
807 */
808 static void
find_big_prime_factor(const MPNumber * n,MPNumber * z)809 find_big_prime_factor (const MPNumber *n, MPNumber *z)
810 {
811 MPNumber tmp = mp_new();
812 unsigned long i = 2;
813
814 while (TRUE)
815 {
816 while (mp_pollard_rho (n, i, &tmp) == FALSE)
817 {
818 i++;
819 }
820
821 if (!mp_is_pprime(&tmp, 50))
822 {
823 mp_divide(n, &tmp, &tmp);
824 }
825 else
826 break;
827 }
828
829 mp_set_from_mp(&tmp, z);
830 mp_clear(&tmp);
831 }
832
833 /**
834 * mp_factorize tries to factorize the value of @x.
835 * If @x < 2^64 it calls mp_factorize_unit64 which deals in integers
836 * and should be fast enough for most values.
837 * If @x > 2^64 the approach to find factors of @x is as follows:
838 * - Try to divide @x by prime numbers 2,3,5,7,.. up to min(2^13, sqrt(x))
839 * - Use Pollard rho to find prime factors > 2^13.
840 * Returns a pointer to a GList with all prime factors of @x which needs to
841 * be freed.
842 */
843 GList*
mp_factorize(const MPNumber * x)844 mp_factorize(const MPNumber *x)
845 {
846 GList *list = NULL;
847 MPNumber *factor = g_slice_alloc0(sizeof(MPNumber));
848 mpc_init2(factor->num, PRECISION);
849
850 MPNumber value = mp_new();
851 mp_abs(x, &value);
852
853 if (mp_is_zero(&value))
854 {
855 mp_set_from_mp(&value, factor);
856 list = g_list_append(list, factor);
857 mp_clear(&value);
858 return list;
859 }
860
861 MPNumber tmp = mp_new();
862 mp_set_from_integer(1, &tmp);
863 if (mp_is_equal(&value, &tmp))
864 {
865 mp_set_from_mp(x, factor);
866 list = g_list_append(list, factor);
867 mp_clear(&value);
868 mp_clear(&tmp);
869 return list;
870 }
871
872 /* If value < 2^64-1, call for factorize_uint64 function which deals in integers */
873 uint64_t num = 1;
874 num = num << 63;
875 num += (num - 1);
876 MPNumber int_max = mp_new();
877 mp_set_from_unsigned_integer(num, &int_max);
878 if (mp_is_less_equal(x, &int_max))
879 {
880 list = mp_factorize_unit64(mp_to_unsigned_integer(&value));
881 if (mp_is_negative(x))
882 mp_invert_sign(list->data, list->data);
883 mp_clear(&value);
884 mp_clear(&tmp);
885 mp_clear(&int_max);
886 return list;
887 }
888
889 MPNumber divisor = mp_new_from_unsigned_integer(2);
890 while (TRUE)
891 {
892 mp_divide(&value, &divisor, &tmp);
893 if (mp_is_integer(&tmp))
894 {
895 mp_set_from_mp(&tmp, &value);
896 mp_set_from_mp(&divisor, factor);
897 list = g_list_append(list, factor);
898 factor = g_slice_alloc0(sizeof(MPNumber));
899 mpc_init2(factor->num, PRECISION);
900 }
901 else
902 break;
903 }
904
905 mp_set_from_integer(3, &divisor);
906
907 MPNumber root = mp_new();
908 mp_sqrt(&value, &root);
909 uint64_t max_trial_division = (uint64_t) (1 << 10);
910 uint64_t iter = 0;
911 while (mp_is_less_equal(&divisor, &root) && (iter++ < max_trial_division))
912 {
913 mp_divide(&value, &divisor, &tmp);
914 if (mp_is_integer(&tmp))
915 {
916 mp_set_from_mp(&tmp, &value);
917 mp_sqrt(&value, &root);
918 mp_set_from_mp(&divisor, factor);
919 list = g_list_append(list, factor);
920 factor = g_slice_alloc0(sizeof(MPNumber));
921 mpc_init2(factor->num, PRECISION);
922 }
923 else
924 {
925 mp_add_integer(&divisor, 2, &divisor);
926 }
927 }
928
929 while (!mp_is_pprime(&value, 50))
930 {
931 find_big_prime_factor (&value, &divisor);
932
933 mp_divide(&value, &divisor, &tmp);
934 if (mp_is_integer(&tmp))
935 {
936 mp_set_from_mp(&tmp, &value);
937 mp_set_from_mp(&divisor, factor);
938 list = g_list_append(list, factor);
939 factor = g_slice_alloc0(sizeof(MPNumber));
940 mpc_init2(factor->num, PRECISION);
941 }
942 }
943
944 mp_set_from_integer(1, &tmp);
945 if (mp_is_greater_than(&value, &tmp))
946 {
947 mp_set_from_mp(&value, factor);
948 list = g_list_append(list, factor);
949 }
950 else
951 {
952 mpc_clear(factor->num);
953 g_slice_free(MPNumber, factor);
954 }
955
956 if (mp_is_negative(x))
957 mp_invert_sign(list->data, list->data);
958
959 mp_clear(&value);
960 mp_clear(&tmp);
961 mp_clear(&divisor);
962 mp_clear(&root);
963
964 return list;
965 }
966
967 GList*
mp_factorize_unit64(uint64_t n)968 mp_factorize_unit64(uint64_t n)
969 {
970 GList *list = NULL;
971 MPNumber *factor = g_slice_alloc0(sizeof(MPNumber));
972 mpc_init2(factor->num, PRECISION);
973
974 MPNumber tmp = mp_new();
975 mp_set_from_unsigned_integer(2, &tmp);
976 while (n % 2 == 0)
977 {
978 n /= 2;
979 mp_set_from_mp(&tmp, factor);
980 list = g_list_append(list, factor);
981 factor = g_slice_alloc0(sizeof(MPNumber));
982 mpc_init2(factor->num, PRECISION);
983 }
984
985 for (uint64_t divisor = 3; divisor <= n / divisor; divisor +=2)
986 {
987 while (n % divisor == 0)
988 {
989 n /= divisor;
990 mp_set_from_unsigned_integer(divisor, factor);
991 list = g_list_append(list, factor);
992 factor = g_slice_alloc0(sizeof(MPNumber));
993 mpc_init2(factor->num, PRECISION);
994 }
995 }
996
997 if (n > 1)
998 {
999 mp_set_from_unsigned_integer(n, factor);
1000 list = g_list_append(list, factor);
1001 }
1002 else
1003 {
1004 mpc_clear(factor->num);
1005 g_slice_free(MPNumber, factor);
1006 }
1007 mp_clear(&tmp);
1008
1009 return list;
1010 }
1011