X-Git-Url: https://oss.titaniummirror.com/gitweb?a=blobdiff_plain;f=gmp%2Fmpz%2Fnextprime.c;fp=gmp%2Fmpz%2Fnextprime.c;h=0fd8da57e91ea87a104e5097f335d183b4487821;hb=6fed43773c9b0ce596dca5686f37ac3fc0fa11c0;hp=0000000000000000000000000000000000000000;hpb=27b11d56b743098deb193d510b337ba22dc52e5c;p=msp430-gcc.git

diff --git a/gmp/mpz/nextprime.c b/gmp/mpz/nextprime.c
new file mode 100644
index 00000000..0fd8da57
--- /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;
+}