← Back to branch summary

~diresu/blender/blender-command-port

« back to all changes in this revision

Viewing changes to extern/fftw/kernel/primes.c

  • Committer: theeth
  • Date: 2008-10-14 16:52:04 UTC
  • Revision ID: vcs-imports@canonical.com-20081014165204-r32w2gm6s0osvdhn
copy back trunk

Show diffs side-by-side

added added

removed removed

Lines of Context:
extern/fftw/kernel/primes.c
 
1
/*
 
2
 * Copyright (c) 2003, 2006 Matteo Frigo
 
3
 * Copyright (c) 2003, 2006 Massachusetts Institute of Technology
 
4
 *
 
5
 * This program is free software; you can redistribute it and/or modify
 
6
 * it under the terms of the GNU General Public License as published by
 
7
 * the Free Software Foundation; either version 2 of the License, or
 
8
 * (at your option) any later version.
 
9
 *
 
10
 * This program is distributed in the hope that it will be useful,
 
11
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 
12
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 
13
 * GNU General Public License for more details.
 
14
 *
 
15
 * You should have received a copy of the GNU General Public License
 
16
 * along with this program; if not, write to the Free Software
 
17
 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
 
18
 *
 
19
 */
 
20
 
 
21
/* $Id: primes.c,v 1.25 2006-02-05 23:19:55 athena Exp $ */
 
22
 
 
23
#include "ifftw.h"
 
24
 
 
25
/***************************************************************************/
 
26
 
 
27
/* Rader's algorithm requires lots of modular arithmetic, and if we
 
28
   aren't careful we can have errors due to integer overflows. */
 
29
 
 
30
/* Compute (x * y) mod p, but watch out for integer overflows; we must
 
31
   have 0 <= {x, y} < p.
 
32
 
 
33
   If overflow is common, this routine is somewhat slower than
 
34
   e.g. using 'long long' arithmetic.  However, it has the advantage
 
35
   of working when INT is 64 bits, and is also faster when overflow is
 
36
   rare.  FFTW calls this via the MULMOD macro, which further
 
37
   optimizes for the case of small integers. 
 
38
*/
 
39
 
 
40
#define ADD_MOD(x, y, p) ((x) >= (p) - (y)) ? ((x) + ((y) - (p))) : ((x) + (y))
 
41
 
 
42
INT X(safe_mulmod)(INT x, INT y, INT p)
 
43
{
 
44
     INT r;
 
45
 
 
46
     if (y > x) 
 
47
          return X(safe_mulmod)(y, x, p);
 
48
 
 
49
     A(0 <= y && x < p);
 
50
 
 
51
     r = 0;
 
52
     while (y) {
 
53
          r = ADD_MOD(r, x*(y&1), p); y >>= 1;
 
54
          x = ADD_MOD(x, x, p);
 
55
     }
 
56
 
 
57
     return r;
 
58
}
 
59
 
 
60
/***************************************************************************/
 
61
 
 
62
/* Compute n^m mod p, where m >= 0 and p > 0.  If we really cared, we
 
63
   could make this tail-recursive. */
 
64
 
 
65
INT X(power_mod)(INT n, INT m, INT p)
 
66
{
 
67
     A(p > 0);
 
68
     if (m == 0)
 
69
          return 1;
 
70
     else if (m % 2 == 0) {
 
71
          INT x = X(power_mod)(n, m / 2, p);
 
72
          return MULMOD(x, x, p);
 
73
     }
 
74
     else
 
75
          return MULMOD(n, X(power_mod)(n, m - 1, p), p);
 
76
}
 
77
 
 
78
/* the following two routines were contributed by Greg Dionne. */
 
79
static INT get_prime_factors(INT n, INT *primef)
 
80
{
 
81
     INT i;
 
82
     INT size = 0;
 
83
 
 
84
     A(n % 2 == 0); /* this routine is designed only for even n */
 
85
     primef[size++] = (INT)2;
 
86
     do
 
87
          n >>= 1;
 
88
     while ((n & 1) == 0);
 
89
 
 
90
     if (n == 1)
 
91
          return size;
 
92
 
 
93
     for (i = 3; i * i <= n; i += 2)
 
94
          if (!(n % i)) {
 
95
               primef[size++] = i;
 
96
               do
 
97
                    n /= i;
 
98
               while (!(n % i));
 
99
          }
 
100
     if (n == 1)
 
101
          return size;
 
102
     primef[size++] = n;
 
103
     return size;
 
104
}
 
105
 
 
106
INT X(find_generator)(INT p)
 
107
{
 
108
    INT n, i, size;
 
109
    INT primef[16];     /* smallest number = 32589158477190044730 > 2^64 */
 
110
    INT pm1 = p - 1;
 
111
 
 
112
    if (p == 2)
 
113
         return 1;
 
114
 
 
115
    size = get_prime_factors(pm1, primef);
 
116
    n = 2;
 
117
    for (i = 0; i < size; i++)
 
118
        if (X(power_mod)(n, pm1 / primef[i], p) == 1) {
 
119
            i = -1;
 
120
            n++;
 
121
        }
 
122
    return n;
 
123
}
 
124
 
 
125
/* Return first prime divisor of n  (It would be at best slightly faster to
 
126
   search a static table of primes; there are 6542 primes < 2^16.)  */
 
127
INT X(first_divisor)(INT n)
 
128
{
 
129
     INT i;
 
130
     if (n <= 1)
 
131
          return n;
 
132
     if (n % 2 == 0)
 
133
          return 2;
 
134
     for (i = 3; i*i <= n; i += 2)
 
135
          if (n % i == 0)
 
136
               return i;
 
137
     return n;
 
138
}
 
139
 
 
140
int X(is_prime)(INT n)
 
141
{
 
142
     return(n > 1 && X(first_divisor)(n) == n);
 
143
}
 
144
 
 
145
INT X(next_prime)(INT n)
 
146
{
 
147
     while (!X(is_prime)(n)) ++n;
 
148
     return n;
 
149
}
 
150
 
 
151
int X(factors_into)(INT n, const INT *primes)
 
152
{
 
153
     for (; *primes != 0; ++primes) 
 
154
          while ((n % *primes) == 0) 
 
155
               n /= *primes;
 
156
     return (n == 1);
 
157
}
 
158
 
 
159
/* integer square root.  Return floor(sqrt(N)) */
 
160
INT X(isqrt)(INT n)
 
161
{
 
162
     INT guess, iguess;
 
163
 
 
164
     A(n >= 0);
 
165
     if (n == 0) return 0;
 
166
 
 
167
     guess = n; iguess = 1;
 
168
 
 
169
     do {
 
170
          guess = (guess + iguess) / 2;
 
171
          iguess = n / guess;
 
172
     } while (guess > iguess);
 
173
 
 
174
     return guess;
 
175
}
 
176
 
 
177
static INT isqrt_maybe(INT n)
 
178
{
 
179
     INT guess = X(isqrt)(n);
 
180
     return guess * guess == n ? guess : 0;
 
181
}
 
182
 
 
183
#define divides(a, b) (((b) % (a)) == 0)
 
184
INT X(choose_radix)(INT r, INT n)
 
185
{
 
186
     if (r > 0) {
 
187
          if (divides(r, n)) return r;
 
188
          return 0;
 
189
     } else if (r == 0) {
 
190
          return X(first_divisor)(n);
 
191
     } else {
 
192
          /* r is negative.  If n = (-r) * q^2, take q as the radix */
 
193
          r = 0 - r;
 
194
          return (n > r && divides(r, n)) ? isqrt_maybe(n / r) : 0;
 
195
     }
 
196
}