First commit - moving from LuaForge to SourceForge

1 files changed, 135 insertions, 0 deletions

diff --git a/src/fftw3/kernel/primes.c b/src/fftw3/kernel/primes.c
new file mode 100644
index 0000000..608e51c
--- /dev/null
+++ b/src/fftw3/kernel/primes.c
@@ -0,0 +1,135 @@
+/*
+ * Copyright (c) 2003 Matteo Frigo
+ * Copyright (c) 2003 Massachusetts Institute of Technology
+ *
+ * This program is free software; you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation; either version 2 of the License, or
+ * (at your option) any later version.
+ *
+ * This program 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 General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program; if not, write to the Free Software
+ * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
+ *
+ */
+
+/* $Id: primes.c,v 1.1 2008/10/17 06:11:29 scuri Exp $ */
+
+#include "ifftw.h"
+
+/***************************************************************************/
+
+/* Rader's algorithm requires lots of modular arithmetic, and if we
+   aren't careful we can have errors due to integer overflows. */
+
+#ifdef SAFE_MULMOD
+
+#  include <limits.h>
+
+/* compute (x * y) mod p, but watch out for integer overflows; we must
+   have x, y >= 0, p > 0.  This routine is slow. */
+int X(safe_mulmod)(int x, int y, int p)
+{
+     if (y == 0 || x <= INT_MAX / y)
+	  return((x * y) % p);
+     else {
+	  int y2 = y/2;
+	  return((X(safe_mulmod)(x, y2, p) +
+		  X(safe_mulmod)(x, y - y2, p)) % p);
+     }
+}
+#endif /* safe_mulmod ('long long' unavailable) */
+
+/***************************************************************************/
+
+/* Compute n^m mod p, where m >= 0 and p > 0.  If we really cared, we
+   could make this tail-recursive. */
+int X(power_mod)(int n, int m, int p)
+{
+     A(p > 0);
+     if (m == 0)
+	  return 1;
+     else if (m % 2 == 0) {
+	  int x = X(power_mod)(n, m / 2, p);
+	  return MULMOD(x, x, p);
+     }
+     else
+	  return MULMOD(n, X(power_mod)(n, m - 1, p), p);
+}
+
+/* the following two routines were contributed by Greg Dionne. */
+static int get_prime_factors(int n, int *primef)
+{
+     int i;
+     int size = 0;
+
+     primef[size++] = 2;
+     do
+	  n >>= 1;
+     while ((n & 1) == 0);
+
+     if (n == 1)
+	  return size;
+
+     for (i = 3; i * i <= n; i += 2)
+	  if (!(n % i)) {
+	       primef[size++] = i;
+	       do
+		    n /= i;
+	       while (!(n % i));
+	  }
+     if (n == 1)
+	  return size;
+     primef[size++] = n;
+     return size;
+}
+
+int X(find_generator)(int p)
+{
+    int n, i, size;
+    int primef[16];     /* smallest number = 32589158477190044730 > 2^64 */
+    int pm1 = p - 1;
+
+    if (p == 2)
+	 return 1;
+
+    size = get_prime_factors(pm1, primef);
+    n = 2;
+    for (i = 0; i < size; i++)
+        if (X(power_mod)(n, pm1 / primef[i], p) == 1) {
+            i = -1;
+            n++;
+        }
+    return n;
+}
+
+/* Return first prime divisor of n  (It would be at best slightly faster to
+   search a static table of primes; there are 6542 primes < 2^16.)  */
+int X(first_divisor)(int n)
+{
+     int i;
+     if (n <= 1)
+	  return n;
+     if (n % 2 == 0)
+	  return 2;
+     for (i = 3; i*i <= n; i += 2)
+	  if (n % i == 0)
+	       return i;
+     return n;
+}
+
+int X(is_prime)(int n)
+{
+     return(n > 1 && X(first_divisor)(n) == n);
+}
+
+int X(next_prime)(int n)
+{
+     while (!X(is_prime)(n)) ++n;
+     return n;
+}