Chris@10: /*
Chris@10:  * Copyright (c) 2003, 2007-11 Matteo Frigo
Chris@10:  * Copyright (c) 2003, 2007-11 Massachusetts Institute of Technology
Chris@10:  *
Chris@10:  * This program is free software; you can redistribute it and/or modify
Chris@10:  * it under the terms of the GNU General Public License as published by
Chris@10:  * the Free Software Foundation; either version 2 of the License, or
Chris@10:  * (at your option) any later version.
Chris@10:  *
Chris@10:  * This program is distributed in the hope that it will be useful,
Chris@10:  * but WITHOUT ANY WARRANTY; without even the implied warranty of
Chris@10:  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
Chris@10:  * GNU General Public License for more details.
Chris@10:  *
Chris@10:  * You should have received a copy of the GNU General Public License
Chris@10:  * along with this program; if not, write to the Free Software
Chris@10:  * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA
Chris@10:  *
Chris@10:  */
Chris@10: 
Chris@10: 
Chris@10: #include "ifftw.h"
Chris@10: 
Chris@10: /***************************************************************************/
Chris@10: 
Chris@10: /* Rader's algorithm requires lots of modular arithmetic, and if we
Chris@10:    aren't careful we can have errors due to integer overflows. */
Chris@10: 
Chris@10: /* Compute (x * y) mod p, but watch out for integer overflows; we must
Chris@10:    have 0 <= {x, y} < p.
Chris@10: 
Chris@10:    If overflow is common, this routine is somewhat slower than
Chris@10:    e.g. using 'long long' arithmetic.  However, it has the advantage
Chris@10:    of working when INT is 64 bits, and is also faster when overflow is
Chris@10:    rare.  FFTW calls this via the MULMOD macro, which further
Chris@10:    optimizes for the case of small integers. 
Chris@10: */
Chris@10: 
Chris@10: #define ADD_MOD(x, y, p) ((x) >= (p) - (y)) ? ((x) + ((y) - (p))) : ((x) + (y))
Chris@10: 
Chris@10: INT X(safe_mulmod)(INT x, INT y, INT p)
Chris@10: {
Chris@10:      INT r;
Chris@10: 
Chris@10:      if (y > x) 
Chris@10: 	  return X(safe_mulmod)(y, x, p);
Chris@10: 
Chris@10:      A(0 <= y && x < p);
Chris@10: 
Chris@10:      r = 0;
Chris@10:      while (y) {
Chris@10: 	  r = ADD_MOD(r, x*(y&1), p); y >>= 1;
Chris@10: 	  x = ADD_MOD(x, x, p);
Chris@10:      }
Chris@10: 
Chris@10:      return r;
Chris@10: }
Chris@10: 
Chris@10: /***************************************************************************/
Chris@10: 
Chris@10: /* Compute n^m mod p, where m >= 0 and p > 0.  If we really cared, we
Chris@10:    could make this tail-recursive. */
Chris@10: 
Chris@10: INT X(power_mod)(INT n, INT m, INT p)
Chris@10: {
Chris@10:      A(p > 0);
Chris@10:      if (m == 0)
Chris@10: 	  return 1;
Chris@10:      else if (m % 2 == 0) {
Chris@10: 	  INT x = X(power_mod)(n, m / 2, p);
Chris@10: 	  return MULMOD(x, x, p);
Chris@10:      }
Chris@10:      else
Chris@10: 	  return MULMOD(n, X(power_mod)(n, m - 1, p), p);
Chris@10: }
Chris@10: 
Chris@10: /* the following two routines were contributed by Greg Dionne. */
Chris@10: static INT get_prime_factors(INT n, INT *primef)
Chris@10: {
Chris@10:      INT i;
Chris@10:      INT size = 0;
Chris@10: 
Chris@10:      A(n % 2 == 0); /* this routine is designed only for even n */
Chris@10:      primef[size++] = (INT)2;
Chris@10:      do
Chris@10: 	  n >>= 1;
Chris@10:      while ((n & 1) == 0);
Chris@10: 
Chris@10:      if (n == 1)
Chris@10: 	  return size;
Chris@10: 
Chris@10:      for (i = 3; i * i <= n; i += 2)
Chris@10: 	  if (!(n % i)) {
Chris@10: 	       primef[size++] = i;
Chris@10: 	       do
Chris@10: 		    n /= i;
Chris@10: 	       while (!(n % i));
Chris@10: 	  }
Chris@10:      if (n == 1)
Chris@10: 	  return size;
Chris@10:      primef[size++] = n;
Chris@10:      return size;
Chris@10: }
Chris@10: 
Chris@10: INT X(find_generator)(INT p)
Chris@10: {
Chris@10:     INT n, i, size;
Chris@10:     INT primef[16];     /* smallest number = 32589158477190044730 > 2^64 */
Chris@10:     INT pm1 = p - 1;
Chris@10: 
Chris@10:     if (p == 2)
Chris@10: 	 return 1;
Chris@10: 
Chris@10:     size = get_prime_factors(pm1, primef);
Chris@10:     n = 2;
Chris@10:     for (i = 0; i < size; i++)
Chris@10:         if (X(power_mod)(n, pm1 / primef[i], p) == 1) {
Chris@10:             i = -1;
Chris@10:             n++;
Chris@10:         }
Chris@10:     return n;
Chris@10: }
Chris@10: 
Chris@10: /* Return first prime divisor of n  (It would be at best slightly faster to
Chris@10:    search a static table of primes; there are 6542 primes < 2^16.)  */
Chris@10: INT X(first_divisor)(INT n)
Chris@10: {
Chris@10:      INT i;
Chris@10:      if (n <= 1)
Chris@10: 	  return n;
Chris@10:      if (n % 2 == 0)
Chris@10: 	  return 2;
Chris@10:      for (i = 3; i*i <= n; i += 2)
Chris@10: 	  if (n % i == 0)
Chris@10: 	       return i;
Chris@10:      return n;
Chris@10: }
Chris@10: 
Chris@10: int X(is_prime)(INT n)
Chris@10: {
Chris@10:      return(n > 1 && X(first_divisor)(n) == n);
Chris@10: }
Chris@10: 
Chris@10: INT X(next_prime)(INT n)
Chris@10: {
Chris@10:      while (!X(is_prime)(n)) ++n;
Chris@10:      return n;
Chris@10: }
Chris@10: 
Chris@10: int X(factors_into)(INT n, const INT *primes)
Chris@10: {
Chris@10:      for (; *primes != 0; ++primes) 
Chris@10: 	  while ((n % *primes) == 0) 
Chris@10: 	       n /= *primes;
Chris@10:      return (n == 1);
Chris@10: }
Chris@10: 
Chris@10: /* integer square root.  Return floor(sqrt(N)) */
Chris@10: INT X(isqrt)(INT n)
Chris@10: {
Chris@10:      INT guess, iguess;
Chris@10: 
Chris@10:      A(n >= 0);
Chris@10:      if (n == 0) return 0;
Chris@10: 
Chris@10:      guess = n; iguess = 1;
Chris@10: 
Chris@10:      do {
Chris@10:           guess = (guess + iguess) / 2;
Chris@10: 	  iguess = n / guess;
Chris@10:      } while (guess > iguess);
Chris@10: 
Chris@10:      return guess;
Chris@10: }
Chris@10: 
Chris@10: static INT isqrt_maybe(INT n)
Chris@10: {
Chris@10:      INT guess = X(isqrt)(n);
Chris@10:      return guess * guess == n ? guess : 0;
Chris@10: }
Chris@10: 
Chris@10: #define divides(a, b) (((b) % (a)) == 0)
Chris@10: INT X(choose_radix)(INT r, INT n)
Chris@10: {
Chris@10:      if (r > 0) {
Chris@10: 	  if (divides(r, n)) return r;
Chris@10: 	  return 0;
Chris@10:      } else if (r == 0) {
Chris@10: 	  return X(first_divisor)(n);
Chris@10:      } else {
Chris@10: 	  /* r is negative.  If n = (-r) * q^2, take q as the radix */
Chris@10: 	  r = 0 - r;
Chris@10: 	  return (n > r && divides(r, n)) ? isqrt_maybe(n / r) : 0;
Chris@10:      }
Chris@10: }
Chris@10: 
Chris@10: /* return A mod N, works for all A including A < 0 */
Chris@10: INT X(modulo)(INT a, INT n)
Chris@10: {
Chris@10:      A(n > 0);
Chris@10:      if (a >= 0)
Chris@10: 	  return a % n;
Chris@10:      else
Chris@10: 	  return (n - 1) - ((-(a + (INT)1)) % n);
Chris@10: }
Chris@10: 
Chris@10: /* TRUE if N factors into small primes */
Chris@10: int X(factors_into_small_primes)(INT n)
Chris@10: {
Chris@10:      static const INT primes[] = { 2, 3, 5, 0 };
Chris@10:      return X(factors_into)(n, primes);
Chris@10: }