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