Imported gcc-4.4.3

[msp430-gcc.git] / gmp / mpz / nextprime.c

diff --git a/gmp/mpz/nextprime.c b/gmp/mpz/nextprime.c
new file mode 100644 (file)
index 0000000..0fd8da5
--- /dev/null
+++ b/gmp/mpz/nextprime.c
@@ -0,0 +1,120 @@
+/* mpz_nextprime(p,t) - compute the next prime > t and store that in p.
+
+Copyright 1999, 2000, 2001, 2008, 2009 Free Software Foundation, Inc.
+
+Contributed to the GNU project by Niels Möller and Torbjörn Granlund.
+
+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 3 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.  If not, see http://www.gnu.org/licenses/.  */
+
+#include "gmp.h"
+#include "gmp-impl.h"
+#include "longlong.h"
+
+static const unsigned char primegap[] =
+{
+  2,2,4,2,4,2,4,6,2,6,4,2,4,6,6,2,6,4,2,6,4,6,8,4,2,4,2,4,14,4,6,
+  2,10,2,6,6,4,6,6,2,10,2,4,2,12,12,4,2,4,6,2,10,6,6,6,2,6,4,2,10,14,4,2,
+  4,14,6,10,2,4,6,8,6,6,4,6,8,4,8,10,2,10,2,6,4,6,8,4,2,4,12,8,4,8,4,6,
+  12,2,18,6,10,6,6,2,6,10,6,6,2,6,6,4,2,12,10,2,4,6,6,2,12,4,6,8,10,8,10,8,
+  6,6,4,8,6,4,8,4,14,10,12,2,10,2,4,2,10,14,4,2,4,14,4,2,4,20,4,8,10,8,4,6,
+  6,14,4,6,6,8,6,12
+};
+
+#define NUMBER_OF_PRIMES 167
+
+void
+mpz_nextprime (mpz_ptr p, mpz_srcptr n)
+{
+  unsigned short *moduli;
+  unsigned long difference;
+  int i;
+  unsigned prime_limit;
+  unsigned long prime;
+  int cnt;
+  mp_size_t pn;
+  unsigned long nbits;
+  unsigned incr;
+  TMP_SDECL;
+
+  /* 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;
+
+  pn = SIZ(p);
+  count_leading_zeros (cnt, PTR(p)[pn - 1]);
+  nbits = pn * GMP_NUMB_BITS - (cnt - GMP_NAIL_BITS);
+  if (nbits / 2 >= NUMBER_OF_PRIMES)
+    prime_limit = NUMBER_OF_PRIMES - 1;
+  else
+    prime_limit = nbits / 2;
+
+  TMP_SMARK;
+
+  /* Compute residues modulo small odd primes */
+  moduli = TMP_SALLOC_TYPE (prime_limit * sizeof moduli[0], unsigned short);
+
+  for (;;)
+    {
+      /* FIXME: Compute lazily? */
+      prime = 3;
+      for (i = 0; i < prime_limit; i++)
+       {
+         moduli[i] = mpz_fdiv_ui (p, prime);
+         prime += primegap[i];
+       }
+
+#define INCR_LIMIT 0x10000     /* deep science */
+
+      for (difference = incr = 0; incr < INCR_LIMIT; difference += 2)
+       {
+         /* First check residues */
+         prime = 3;
+         for (i = 0; i < prime_limit; i++)
+           {
+             unsigned r;
+             /* FIXME: Reduce moduli + incr and store back, to allow for
+                division-free reductions.  Alternatively, table primes[]'s
+                inverses (mod 2^16).  */
+             r = (moduli[i] + incr) % prime;
+             prime += primegap[i];
+
+             if (r == 0)
+               goto next;
+           }
+
+         mpz_add_ui (p, p, difference);
+         difference = 0;
+
+         /* Miller-Rabin test */
+         if (mpz_millerrabin (p, 10))
+           goto done;
+       next:;
+         incr += 2;
+       }
+      mpz_add_ui (p, p, difference);
+      difference = 0;
+    }
+ done:
+  TMP_SFREE;
+}