summaryrefslogtreecommitdiff
path: root/rts/gmp/mpz/nextprime.c
diff options
context:
space:
mode:
Diffstat (limited to 'rts/gmp/mpz/nextprime.c')
-rw-r--r--rts/gmp/mpz/nextprime.c120
1 files changed, 120 insertions, 0 deletions
diff --git a/rts/gmp/mpz/nextprime.c b/rts/gmp/mpz/nextprime.c
new file mode 100644
index 0000000000..f024dd1206
--- /dev/null
+++ b/rts/gmp/mpz/nextprime.c
@@ -0,0 +1,120 @@
+/* mpz_nextprime(p,t) - compute the next prime > t and store that in p.
+
+Copyright (C) 1999, 2000 Free Software Foundation, Inc.
+
+This file is part of the GNU MP Library.
+
+The GNU MP Library is free software; you can redistribute it and/or modify
+it under the terms of the GNU Lesser General Public License as published by
+the Free Software Foundation; either version 2.1 of the License, or (at your
+option) any later version.
+
+The GNU MP Library is distributed in the hope that it will be useful, but
+WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
+or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
+License for more details.
+
+You should have received a copy of the GNU Lesser General Public License
+along with the GNU MP Library; see the file COPYING.LIB. If not, write to
+the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
+MA 02111-1307, USA. */
+
+#include "gmp.h"
+#include "gmp-impl.h"
+
+void
+#if __STDC__
+mpz_nextprime (mpz_ptr p, mpz_srcptr t)
+#else
+mpz_nextprime (p, t)
+ mpz_ptr p;
+ mpz_srcptr t;
+#endif
+{
+ mpz_add_ui (p, t, 1L);
+ while (! mpz_probab_prime_p (p, 5))
+ mpz_add_ui (p, p, 1L);
+}
+
+#if 0
+/* This code is not yet tested. Will be enabled in 3.1. */
+
+status unsigned short primes[] =
+{
+3,5,7,11,13,17,19,23,29,31,37,41,43,47,53,59,61,67,71,73,79,83,89,97,
+101,103,107,109,113,127,131,137,139,149,151,157,163,167,173,179,181,
+191,193,197,199,211,223,227,229,233,239,241,251,257,263,269,271,277,
+281,283,293,307,311,313,317,331,337,347,349,353,359,367,373,379,383,
+389,397,401,409,419,421,431,433,439,443,449,457,461,463,467,479,487,
+491,499,503,509,521,523,541,547,557,563,569,571,577,587,593,599,601,
+607,613,617,619,631,641,643,647,653,659,661,673,677,683,691,701,709,
+719,727,733,739,743,751,757,761,769,773,787,797,809,811,821,823,827,
+829,839,853,857,859,863,877,881,883,887,907,911,919,929,937,941,947,
+953,967,971,977,983,991,997
+};
+
+#define NUMBER_OF_PRIMES 167
+
+void
+#if __STDC__
+mpz_nextprime (mpz_ptr p, mpz_srcptr n)
+#else
+mpz_nextprime (p, n)
+ mpz_ptr p;
+ mpz_srcptr n;
+#endif
+{
+ mpz_t tmp;
+ unsigned short *moduli;
+ unsigned long difference;
+ int i;
+ int composite;
+
+ /* First handle tiny numbers */
+ if (mpz_cmp_ui (n, 2) < 0)
+ {
+ mpz_set_ui (p, 2);
+ return;
+ }
+ mpz_add_ui (p, n, 1);
+ mpz_setbit (p, 0);
+
+ if (mpz_cmp_ui (p, 7) <= 0)
+ return;
+
+ prime_limit = NUMBER_OF_PRIMES - 1;
+ if (mpz_cmp_ui (p, primes[prime_limit]) <= 0)
+ /* Just use first three entries (3,5,7) of table for small numbers */
+ prime_limit = 3;
+ if (prime_limit)
+ {
+ /* Compute residues modulo small odd primes */
+ moduli = (unsigned short *) TMP_ALLOC (prime_limit * sizeof moduli[0]);
+ for (i = 0; i < prime_limit; i++)
+ moduli[i] = mpz_fdiv_ui (p, primes[i]);
+ }
+ for (difference = 0; ; difference += 2)
+ {
+ composite = 0;
+
+ /* First check residues */
+ for (i = 0; i < prime_limit; i++)
+ {
+ int acc, pr;
+ composite |= (moduli[i] == 0);
+ acc = moduli[i] + 2;
+ pr = primes[i];
+ moduli[i] = acc >= pr ? acc - pr : acc;
+ }
+ if (composite)
+ continue;
+
+ mpz_add_ui (p, p, difference);
+ difference = 0;
+
+ /* Miller-Rabin test */
+ if (mpz_millerrabin (p, 2))
+ break;
+ }
+}
+#endif