1 // Fast integer multiplication using FFT over the complex numbers, modified
2 // to work with real numbers only in the implementation.
3 // [Donald Ervin Knuth: The Art of Computer Programming, Vol. II:
4 // Seminumerical Algorithms, second edition. Section 4.3.3, p. 290-294.]
5 // Bruno Haible 6.5.1996, 24.-25.8.1996, 27.-30.8.1996
7 // FFT in the complex domain has the drawback that it needs careful round-off
8 // error analysis. But for CPUs with good floating-point performance it might
9 // nevertheless be better than FFT mod Z/mZ.
11 // The usual FFT(2^n) computes the values of a polynomial p(z) mod (z^(2^n)-1)
12 // at the (2^n)-th roots of unity. For our purposes, we start out with
13 // polynomials with real coefficients. So the values at z_j = exp(2 pi i j/N)
14 // and z_-j = exp(- 2 pi i j/N) will be complex conjugates of each other,
15 // which means that there exists a polynomial r(z) of degree <= 1 with real
16 // coefficients such that p(z_j) = r(z_j) and p(z_-j) = r(z_-j). This
17 // implies that r(z) is the remainder of p(z) divided by
18 // (z - z_j) (z - z_-j) = (z^2 - 2 cos(2 pi j/N) z + 1).
20 // Based on this insight, we replace the usual n FFT steps
22 // (z^(2^n)-1) = prod(j=0..2^(n-m)-1, (z^(2^m) - exp(2 pi j/2^(n-m))))
25 // (z^(2^n)-1) = (z^(2^m)-1) * prod(j=1..2^(n-m)-1, factor_m[j](z))
27 // factor_m[j](z) = prod(k mod 2^n with k == j or k == -j mod 2^(n-m+1),
28 // (z - exp(2 pi i k/N)) )
29 // = prod(k=0..2^(m-1)-1,
30 // (z - exp(2 pi i (j+2^(n-m+1)k)/N))
31 // (z - exp(- 2 pi i (j+2^(n-m+1)k)/N)) )
32 // = (z^(2^(m-1)) - exp(2 pi i j 2^(m-1)/N))
33 // (z^(2^(m-1)) - exp(- 2 pi i j 2^(m-1)/N))
34 // = (z^(2^(m-1)) - exp(2 pi i j/2^(n-m+1)))
35 // (z^(2^(m-1)) - exp(- 2 pi i j/2^(n-m+1)))
36 // = (z^(2^m) - 2 cos(2 pi j/2^(n-m+1)) z^(2^(m-1)) + 1).
37 // The factors and the input are real polynomials, hence all intermediate
38 // and final remainders will be real as well.
40 // The usual FFT algorithm
41 // Input: polynomial p in x[0..2^n-1].
42 // for l = n-1..0: // step m=l+1 -> m=l
43 // for s in {0,..,2^(n-1-l)-1}:
44 // exp := bit_reverse(n-1-l,s)*2^l,
45 // // chinese remainder algorithm for (z^(2^(l+1)) - w^(2*exp)) =
46 // // = (z^(2^l) - w^exp) * (z^(2^l) - w^(exp+2^(n-1))).
47 // for t in {0,..,2^l-1}:
48 // i1 := s*2^(l+1) + t, i2 := s*2^(l+1) + 2^l + t,
49 // replace (x[i1],x[i2]) by (x[i1] + w^exp*x[i2], x[i1] - w^exp*x[i2])
52 // for j in {0..2^(n-m)-1}:
53 // p(z) mod (z^(2^m) - exp(2 pi i j/2^(n-m)))
54 // in x[bit_reverse(n-m,j)*2^m .. bit_reverse(n-m,j)*2^m+2^m-1].
55 // Output: p(z_j) in x[bit_reverse(n,j)].
56 // is thus replaced by the algorithm
57 // Input: polynomial p in x[0..2^n-1].
58 // for l = n-1..1: // step m=l+1 -> m=l
60 // // chinese remainder algorithm for (z^(2^(l+1)) - 1) =
61 // // = (z^(2^l) - 1) * (z^(2^l) + 1).
62 // for t in {0,..,2^l-1}:
63 // i1 := t, i2 := 2^l + t,
64 // replace (x[i1],x[i2]) by (x[i1] + x[i2], x[i1] - x[i2])
65 // for s in {1,..,2^(n-1-l)-1}:
66 // exp := shuffle(n-1-l,s)*2^(l-1),
67 // // chinese remainder algorithm for
68 // // (z^(2^(l+1)) - 2 cos(2 pi 2*exp/N) z^(2^l) + 1) =
69 // // = (z^(2^l) - 2 cos(2 pi exp/N) z^(2^(l-1)) + 1)
70 // // * (z^(2^l) - 2 cos(2 pi (2^(n-1)-exp)/N) z^(2^(l-1)) + 1)
71 // gam := 2 cos(2 pi exp/N),
72 // for t in {0,..,2^(l-1)-1}:
73 // i1 := s*2^(l+1) + t, i2 := s*2^(l+1) + 2^(l-1) + t,
74 // i3 := s*2^(l+1) + 2^l + t, i4 := s*2^(l+1) + 2^l + 2^(l-1) + t,
75 // replace (x[i1],x[i2],x[i3],x[i4]) by
76 // (x[i1]-x[i3] - x[i4]*gam, x[i3]*gam + x[i2]+x[i4]*(gam^2-1),
77 // x[i1]-x[i3] + x[i4]*gam, - x[i3]*gam + x[i2]+x[i4]*(gam^2-1))
80 // p(z) mod (z^(2^m) - 1) in x[0..2^m-1],
81 // for j in {1,..,2^(n-m)-1}:
82 // p(z) mod (z^(2^m) - 2 cos(2 pi j/2^(n-m+1)) z^(2^(m-1)) + 1)
83 // in x[invshuffle(n-m,j)*2^m .. invshuffle(n-m,j)*2^m+2^m-1].
84 // Output: p(z) mod (z^2 - 1) in x[0],x[1],
85 // p(z) mod (z^2 - 2 cos(2 pi j/2^n) z + 1) (0 < j < 2^(n-1))
86 // in x[2*invshuffle(n-1,j)],x[2*invshuffle(n-1,j)+1].
88 // The shuffle function is defined like this:
89 // shuffle(n,j) defined for n >= 0, 0 < j < 2^n, yields 0 < shuffle(n,j) < 2^n.
90 // Definition by splitting off the least significant bit:
92 // n > 0: shuffle(n,1) = 2^(n-1),
93 // n > 0, 0 < j < 2^(n-1): shuffle(n,2*j) = shuffle(n-1,j),
94 // n > 0, 0 < j < 2^(n-1): shuffle(n,2*j+1) = 2^n - shuffle(n-1,j).
95 // Its inverse function is defined like this:
96 // invshuffle(n,j) defined for n >= 0, 0 < j < 2^n, 0 < invshuffle(n,j) < 2^n.
97 // Definition by splitting off the most significant bit:
99 // n > 0, 0 < j < 2^(n-1): invshuffle(n,j) = invshuffle(n-1,j)*2,
100 // n > 0, j = 2^(n-1): invshuffle(n,j) = 1,
101 // n > 0, 2^(n-1) < j < 2^n: invshuffle(n,j) = invshuffle(n-1,2^n-j)*2+1.
102 // Note that shuffle(n,.) and invshuffle(n,.) are _not_ the same permutation
105 // It is important to have precise round-off error estimates. Although
106 // (by laws of statistics) in the average the actual round-off error makes up
107 // only half of the bits provided for round-off protection, we cannot rely
108 // on this average behaviour, but have to produce correct results.
110 // Knuth's formula (42), p. 294, says:
111 // If we want to multiply l-bit words using an FFT(2^k), our floating point
112 // numbers shall have m >= 2(k+l) + k + log_2 k + 3.5 mantissa bits.
114 // Here is a more careful analysis, using absolute error estimates.
116 // 1. We assume floating point numbers with radix 2, with the properties:
117 // (i) Multiplication with 2^n and 2^-n is exact.
118 // (ii) Negation x -> -x is exact.
119 // (iii) Addition: When adding x and y, with |x| <= 2^a, |y| <= 2^a,
120 // the result |x+y| <= 2^(a+1) has an error <= e*2^(a+1-m).
121 // (iv) Multiplication: When multiplying x and y, with |x| <= 2^a,
122 // |y| <= 2^b, the result |x*y| <= 2^(a+b) has an error <= e*2^(a+b-m).
123 // (v) Division: When dividing x by y, with |x| <= 2^a,
124 // 2^b <= |y| <= 2^(b+1), the result |x/y| <= 2^(a-b) has an error
126 // Here e = 1 for a truncating arithmetic, but e = 1/2 for a rounding
127 // arithmetic like IEEE single and double floats.
128 // 2. Let's introduce some notation: err(x) means |x'-x| where x is the
129 // exact mathematical value and x' is its representation in the machine.
130 // 3. From 1. we get for real numbers x,y:
131 // (i) err(2^n*x) = 2^n * err(x),
132 // (ii) err(-x) = err(x),
133 // (iii) |x| <= 2^a, |y| <= 2^a, then
134 // err(x+y) <= err(x) + err(y) + e*2^(a+1-m),
135 // [or .... ............... + e*2^(a-m) if |x+y| <= 2^a].
136 // (iv) |x| <= 2^a, |y| <= 2^b, then
137 // err(x*y) <= 2^a * err(y) + 2^b * err(x) + e*2^(a+b-m).
138 // (v) |x| <= 2^a, 2^b <= |y| <= 2^(b+1), then
139 // err(x/y) <= 2^(-b) * err(x) + 2^(a-2b) * err(y) + e*2^(a-b-m).
140 // 4. Our complex arithmetic will be based on the formulas:
141 // (i) 2^n*(x+iy) = (2^n*x)+i(2^n*y)
142 // (ii) -(x+iy) = (-x)+i(-y)
143 // (iii) (x+iy)+(u+iv) = (x+u)+i(y+v)
144 // (iv) (x+iy)*(u+iv) = (x*u-y*v)+i(x*v+y*u)
145 // The notation err(z) means |z'-z|, as above, with |.| being the usual
146 // absolute value on complex numbers (_not_ the L^1 norm).
147 // 5. From 3. and 4. we get for complex numbers x,y:
148 // (i) err(2^n*x) = 2^n * err(x),
149 // (ii) err(-x) = err(x),
150 // (iii) |x| <= 2^a, |y| <= 2^a, then
151 // err(x+y) <= err(x) + err(y) + e*2^(a+3/2-m),
152 // (iv) |x| <= 2^a, |y| <= 2^b, then
153 // err(x*y) <= 2^a * err(y) + 2^b * err(x) + 3*e*2^(a+b+1/2-m).
154 // 6. We start out with precomputed roots of unity:
155 // |exp(2 pi i/2^n)| <= 1,
156 // err(exp(2 pi i/2^n)) <= e*2^(1/2-m), (even err(..)=0 for n=0,1,2),
157 // and compute exp(2 pi i * j/2^k) according to the binary digits of j.
158 // This way, each root of unity will be a product of at most k precomputed
159 // roots. If (j mod 2^(k-2)) has n bits, then exp(2 pi i * j/2^k) will
160 // be computed using n factors, i.e. n-1 complex multiplications, and by
162 // err(exp(2 pi i * j/2^k)) <= n*e*2^(1/2-m) + max(n-1,0)*3*e*2^(1/2-m)
163 // = max(4*n-3,0)*e*2^(1/2-m).
164 // Hence the maximum roots-of-unity error is (set n=k-2)
165 // err(w^j) <= (4*k-11)*e*2^(1/2-m),
166 // and the average roots-of-unity error is (set n=(k-2)/2)
167 // < 2*(k-2)*e*2^(1/2-m).
168 // 7. We use precomputed values of gam = 2*cos(2*pi*j/2^k) (0 < j < 2^(k-2)),
169 // 12*2^-k <= |gam| <= 2,
170 // err(gam) <= (4*k-11)*e*2^(3/2-m),
172 // err(gam-1) <= (4*k-11)*e*2^(3/2-m),
173 // err(gam+1) <= (4*k-9)*e*2^(3/2-m),
174 // err(gam^2-1) <= (20*k-51)*e*2^(3/2-m),
175 // err(1/gam) <= (4*k-10)*e*(2*k-9/2-m).
176 // 8. Now we start the FFT.
177 // Before the first step, x_i are integral, |x_i| < 2^l and err(x_i) = 0.
178 // After the first butterfly, which replaces (x(i1),x(i2)) by
179 // (x(i1) + x(i2), x(i1) - x(i2)), we have |x_i| < 2^(l+1) and err(x_i) = 0.
180 // Then, for each of the remaining k-2 steps, a butterfly replaces
181 // (x[i1],x[i2],x[i3],x[i4]) by
182 // (x[i1]-x[i3] - x[i4]*gam, x[i3]*gam + x[i2]+x[i4]*(gam^2-1),
183 // x[i1]-x[i3] + x[i4]*gam, - x[i3]*gam + x[i2]+x[i4]*(gam^2-1)).
184 // Thus, after n steps we have |x_i| < 2^(l+3n) and err(x_i) <= E_i where
189 // This is just too complicated. We use fftc's round-off error estimates.
190 // As a consequence, sometimes the algorithm bails out, saying
191 // "FFT problem: checksum error"!!!
195 #error "real fft implemented only for intDsize==32"
199 #include "cln/floatparam.h"
200 #include "cln/exception.h"
203 #if (long_double_mant_bits > double_mant_bits) && (defined(__i386__) || defined(__m68k__) || (defined(__sparc__) && 0))
204 // Only these CPUs have fast "long double"s in hardware.
205 // On SPARC, "long double"s are emulated in software and don't work.
206 typedef long double fftr_real;
207 #define fftr_real_mant_bits long_double_mant_bits
208 #define fftr_real_rounds long_double_rounds
210 typedef double fftr_real;
211 #define fftr_real_mant_bits double_mant_bits
212 #define fftr_real_rounds double_rounds
215 // We need complex numbers for precomputing the roots of unity.
216 typedef struct fftr_complex {
221 // For every integer exp, we need to precompute gam = 2 cos(2 pi exp/N).
223 fftr_complex c; // exp(2 pi i exp/N)
225 fftr_real d; // gam = 2 cos(2 pi exp/N)
226 fftr_real e; // gam^2-1
227 fftr_real f; // 1/gam
231 static const fftr_complex fftr_roots_of_1 [32+1] =
232 // roots_of_1[n] is a (2^n)th root of unity in C.
233 // Also roots_of_1[n-1] = roots_of_1[n]^2.
234 // For simplicity we choose roots_of_1[n] = exp(2 pi i/2^n).
236 #if (fftr_real_mant_bits == double_mant_bits)
237 // These values have 64 bit precision.
241 { 0.7071067811865475244, 0.7071067811865475244 },
242 { 0.9238795325112867561, 0.38268343236508977172 },
243 { 0.9807852804032304491, 0.19509032201612826784 },
244 { 0.99518472667219688623, 0.098017140329560601996 },
245 { 0.9987954562051723927, 0.049067674327418014254 },
246 { 0.9996988186962042201, 0.024541228522912288032 },
247 { 0.99992470183914454094, 0.0122715382857199260795 },
248 { 0.99998117528260114264, 0.0061358846491544753597 },
249 { 0.9999952938095761715, 0.00306795676296597627 },
250 { 0.99999882345170190993, 0.0015339801862847656123 },
251 { 0.99999970586288221914, 7.6699031874270452695e-4 },
252 { 0.99999992646571785114, 3.8349518757139558907e-4 },
253 { 0.9999999816164292938, 1.9174759731070330744e-4 },
254 { 0.9999999954041073129, 9.5873799095977345874e-5 },
255 { 0.99999999885102682754, 4.793689960306688455e-5 },
256 { 0.99999999971275670683, 2.3968449808418218729e-5 },
257 { 0.9999999999281891767, 1.1984224905069706422e-5 },
258 { 0.99999999998204729416, 5.9921124526424278428e-6 },
259 { 0.99999999999551182357, 2.9960562263346607504e-6 },
260 { 0.99999999999887795586, 1.4980281131690112288e-6 },
261 { 0.999999999999719489, 7.4901405658471572114e-7 },
262 { 0.99999999999992987223, 3.7450702829238412391e-7 },
263 { 0.99999999999998246807, 1.8725351414619534487e-7 },
264 { 0.999999999999995617, 9.36267570730980828e-8 },
265 { 0.99999999999999890425, 4.6813378536549092695e-8 },
266 { 0.9999999999999997261, 2.340668926827455276e-8 },
267 { 0.99999999999999993153, 1.1703344634137277181e-8 },
268 { 0.99999999999999998287, 5.8516723170686386908e-9 },
269 { 0.9999999999999999957, 2.925836158534319358e-9 },
270 { 0.9999999999999999989, 1.4629180792671596806e-9 }
271 #else (fftr_real_mant_bits > double_mant_bits)
272 // These values have 128 bit precision.
276 { 0.707106781186547524400844362104849039284L, 0.707106781186547524400844362104849039284L },
277 { 0.923879532511286756128183189396788286823L, 0.38268343236508977172845998403039886676L },
278 { 0.980785280403230449126182236134239036975L, 0.195090322016128267848284868477022240928L },
279 { 0.995184726672196886244836953109479921574L, 0.098017140329560601994195563888641845861L },
280 { 0.998795456205172392714771604759100694444L, 0.0490676743274180142549549769426826583147L },
281 { 0.99969881869620422011576564966617219685L, 0.0245412285229122880317345294592829250654L },
282 { 0.99992470183914454092164649119638322435L, 0.01227153828571992607940826195100321214037L },
283 { 0.999981175282601142656990437728567716173L, 0.00613588464915447535964023459037258091705L },
284 { 0.999995293809576171511580125700119899554L, 0.00306795676296597627014536549091984251894L },
285 { 0.99999882345170190992902571017152601905L, 0.001533980186284765612303697150264079079954L },
286 { 0.999999705862882219160228217738765677117L, 7.66990318742704526938568357948576643142e-4L },
287 { 0.99999992646571785114473148070738785695L, 3.83495187571395589072461681181381263396e-4L },
288 { 0.999999981616429293808346915402909714504L, 1.91747597310703307439909561989000933469e-4L },
289 { 0.99999999540410731289097193313960614896L, 9.58737990959773458705172109764763511872e-5L },
290 { 0.9999999988510268275626733077945541084L, 4.79368996030668845490039904946588727468e-5L },
291 { 0.99999999971275670684941397221864177609L, 2.39684498084182187291865771650218200947e-5L },
292 { 0.999999999928189176709775095883850490262L, 1.198422490506970642152156159698898480473e-5L },
293 { 0.99999999998204729417728262414778410738L, 5.99211245264242784287971180889086172999e-6L },
294 { 0.99999999999551182354431058417299732444L, 2.99605622633466075045481280835705981183e-6L },
295 { 0.999999999998877955886077016551752536504L, 1.49802811316901122885427884615536112069e-6L },
296 { 0.999999999999719488971519214794719584451L, 7.49014056584715721130498566730655637157e-7L },
297 { 0.99999999999992987224287980123972873676L, 3.74507028292384123903169179084633177398e-7L },
298 { 0.99999999999998246806071995015624773673L, 1.8725351414619534486882457659356361712e-7L },
299 { 0.999999999999995617015179987529456656217L, 9.3626757073098082799067286680885620193e-8L },
300 { 0.999999999999998904253794996881763834182L, 4.68133785365490926951155181385400969594e-8L },
301 { 0.99999999999999972606344874922040343793L, 2.34066892682745527595054934190348440379e-8L },
302 { 0.999999999999999931515862187305098514444L, 1.170334463413727718124621350323810379807e-8L },
303 { 0.999999999999999982878965546826274482047L, 5.8516723170686386908097901008341396944e-9L },
304 { 0.999999999999999995719741386706568611352L, 2.92583615853431935792823046906895590202e-9L },
305 { 0.999999999999999998929935346676642152265L, 1.46291807926715968052953216186596371037e-9L }
309 // Define this for (cheap) consistency checks.
312 static fftr_real fftr_pow2_table[64] = // table of powers of 2
371 144115188075855872.0,
372 288230376151711744.0,
373 576460752303423488.0,
374 1152921504606846976.0,
375 2305843009213693952.0,
376 4611686018427387904.0,
377 9223372036854775808.0
380 // For a constant expression n (0 <= n < 128), returns 2^n of type fftr_real.
381 #define fftr_pow2(n) \
382 (((n) & 64 ? (fftr_real)18446744073709551616.0 : (fftr_real)1.0) \
383 * ((n) & 32 ? (fftr_real)4294967296.0 : (fftr_real)1.0) \
384 * ((n) & 16 ? (fftr_real)65536.0 : (fftr_real)1.0) \
385 * ((n) & 8 ? (fftr_real)256.0 : (fftr_real)1.0) \
386 * ((n) & 4 ? (fftr_real)16.0 : (fftr_real)1.0) \
387 * ((n) & 2 ? (fftr_real)4.0 : (fftr_real)1.0) \
388 * ((n) & 1 ? (fftr_real)2.0 : (fftr_real)1.0) \
392 static inline void mul (const fftr_complex& a, const fftr_complex& b, fftr_complex& r)
394 var fftr_real r_re = a.re * b.re - a.im * b.im;
395 var fftr_real r_im = a.re * b.im + a.im * b.re;
400 static uintC shuffle (uintL n, uintC x)
404 // Invariant: y + v*shuffle(n,x).
408 return y + (v << (n-1));
414 } while (!(--n == 0));
415 throw runtime_exception();
419 static uintC invshuffle (uintL n, uintC x)
423 // Invariant: y + v*invshuffle(n,x).
425 if (x == ((uintC)1 << (n-1)))
427 else if (x > ((uintC)1 << (n-1))) {
428 x = ((uintC)1 << n) - x;
432 } while (!(--n == 0));
433 throw runtime_exception();
437 // Compute a real convolution using FFT: z[0..N-1] := x[0..N-1] * y[0..N-1].
438 static void fftr_convolution (const uintL n, const uintC N, // N = 2^n
439 fftr_real * x, // N numbers
440 fftr_real * y, // N numbers
441 fftr_real * z // N numbers result
445 var fftr_cosinus* const w = cl_alloc_array(fftr_cosinus,N>>2);
447 // Initialize w[exp].c to exp(2 pi i exp/N).
448 w[0].c = fftr_roots_of_1[0];
451 for (j = n-3; j>=0; j--) {
452 var fftr_complex r_j = fftr_roots_of_1[n-j];
454 for (i = (2<<j); i < N>>2; i += (2<<j))
455 mul(w[i].c,r_j, w[i+(1<<j)].c);
458 // Initialize w[exp] to 2 cos(2 pi exp/N).
459 for (i = 0; i < N>>2; i++) {
460 var fftr_real gam = w[i].c.re * (fftr_real)2.0;
462 w[i].gam.e = (gam - (fftr_real)1.0) * (gam + (fftr_real)1.0);
463 w[i].gam.f = (fftr_real)1.0 / gam;
465 var bool squaring = (x == y);
466 // Do an FFT of length N on x.
469 for (l = n-1; l > 0; l--) {
471 var const uintC tmax = (uintC)1 << l;
472 for (var uintC t = 0; t < tmax; t++) {
474 var uintC i2 = i1 + tmax;
475 // replace (x[i1],x[i2]) by
476 // (x[i1] + x[i2], x[i1] - x[i2])
483 var const uintC smax = (uintC)1 << (n-1-l);
484 var const uintC tmax = (uintC)1 << (l-1);
485 for (var uintC s = 1; s < smax; s++) {
486 var uintC exp = shuffle(n-1-l,s) << (l-1);
487 for (var uintC t = 0; t < tmax; t++) {
488 var uintC i1 = (s << (l+1)) + t;
489 var uintC i2 = i1 + tmax;
490 var uintC i3 = i2 + tmax;
491 var uintC i4 = i3 + tmax;
492 // replace (x[i1],x[i2],x[i3],x[i4]) by
493 // (x[i1]-x[i3] - x[i4]*gam,
494 // x[i3]*gam + x[i2]+x[i4]*(gam^2-1),
495 // x[i1]-x[i3] + x[i4]*gam,
496 // - x[i3]*gam + x[i2]+x[i4]*(gam^2-1))
501 sum13 = x[i1] - x[i3];
502 tmp3 = x[i3]*w[exp].gam.d;
503 tmp4 = x[i4]*w[exp].gam.d;
504 sum24 = x[i2]+x[i4]*w[exp].gam.e;
505 x[i1] = sum13 - tmp4;
506 x[i3] = sum13 + tmp4;
507 x[i2] = sum24 + tmp3;
508 x[i4] = sum24 - tmp3;
513 // replace (x[0],x[1]) by (x[0]+x[1], x[0]-x[1])
520 // Do an FFT of length N on y.
523 for (l = n-1; l > 0; l--) {
525 var const uintC tmax = (uintC)1 << l;
526 for (var uintC t = 0; t < tmax; t++) {
528 var uintC i2 = i1 + tmax;
529 // replace (y[i1],y[i2]) by
530 // (y[i1] + y[i2], y[i1] - y[i2])
537 var const uintC smax = (uintC)1 << (n-1-l);
538 var const uintC tmax = (uintC)1 << (l-1);
539 for (var uintC s = 1; s < smax; s++) {
540 var uintC exp = shuffle(n-1-l,s) << (l-1);
541 for (var uintC t = 0; t < tmax; t++) {
542 var uintC i1 = (s << (l+1)) + t;
543 var uintC i2 = i1 + tmax;
544 var uintC i3 = i2 + tmax;
545 var uintC i4 = i3 + tmax;
546 // replace (y[i1],y[i2],y[i3],y[i4]) by
547 // (y[i1]-y[i3] - y[i4]*gam,
548 // y[i3]*gam + y[i2]+y[i4]*(gam^2-1),
549 // y[i1]-y[i3] + y[i4]*gam,
550 // - y[i3]*gam + y[i2]+y[i4]*(gam^2-1))
555 sum13 = y[i1] - y[i3];
556 tmp3 = y[i3]*w[exp].gam.d;
557 tmp4 = y[i4]*w[exp].gam.d;
558 sum24 = y[i2]+y[i4]*w[exp].gam.e;
559 y[i1] = sum13 - tmp4;
560 y[i3] = sum13 + tmp4;
561 y[i2] = sum24 + tmp3;
562 y[i4] = sum24 - tmp3;
567 // replace (y[0],y[1]) by (y[0]+y[1], y[0]-y[1])
574 // Multiply the transformed vectors into z.
576 // Multiplication mod (z-1).
578 // Multiplication mod (z+1).
580 #if 0 // This needs w[1..2^(n-1)-1].
581 var const uintC smax = (uintC)1 << (n-1);
582 for (var uintC s = 1; s < smax; s++) {
583 // Multiplication mod (z^2 - 2 cos(2 pi j/2^n) z + 1)
584 // with j = shuffle(n-1,s).
585 var uintC exp = shuffle(n-1,s);
586 var fftr_real tmp0 = x[2*s] * y[2*s];
587 var fftr_real tmp1 = x[2*s] * y[2*s+1] + x[2*s+1] * y[2*s];
588 var fftr_real tmp2 = x[2*s+1] * y[2*s+1];
589 z[2*s] = tmp0 - tmp2;
590 z[2*s+1] = tmp1 + tmp2*w[exp].gam.d;
592 #else // This only needs w[1..2^(n-2)-1].
593 // Multiplication mod (z^2+1).
595 var fftr_real tmp0 = x[2] * y[2];
596 var fftr_real tmp1 = x[2] * y[3] + x[3] * y[2];
597 var fftr_real tmp2 = x[3] * y[3];
601 var const uintC smax = (uintC)1 << (n-2);
602 for (var uintC s = 1; s < smax; s++) {
603 var uintC exp = shuffle(n-2,s);
604 // Multiplication mod (z^2 - 2 cos(2 pi j/2^n) z + 1)
605 // with j = shuffle(n-1,2*s) = shuffle(n-2,s).
606 var fftr_real gam = w[exp].gam.d;
608 var fftr_real tmp0 = x[4*s] * y[4*s];
609 var fftr_real tmp1 = x[4*s] * y[4*s+1] + x[4*s+1] * y[4*s];
610 var fftr_real tmp2 = x[4*s+1] * y[4*s+1];
611 z[4*s] = tmp0 - tmp2;
612 z[4*s+1] = tmp1 + tmp2 * gam;
614 // Multiplication mod (z^2 - 2 cos(2 pi j/2^n) z + 1)
615 // with j = shuffle(n-1,2*s+1) = 2^(n-1) - shuffle(n-2,s).
617 var fftr_real tmp0 = x[4*s+2] * y[4*s+2];
618 var fftr_real tmp1 = x[4*s+2] * y[4*s+3] + x[4*s+3] * y[4*s+2];
619 var fftr_real tmp2 = x[4*s+3] * y[4*s+3];
620 z[4*s+2] = tmp0 - tmp2;
621 z[4*s+3] = tmp1 - tmp2 * gam;
626 // Undo an FFT of length N on z.
630 // replace (z[0],z[1]) by ((z[0]+z[1])/2, (z[0]-z[1])/2)
633 z[1] = (z[0] - tmp) * (fftr_real)0.5;
634 z[0] = (z[0] + tmp) * (fftr_real)0.5;
636 for (l = 1; l < n; l++) {
638 var const uintC tmax = (uintC)1 << l;
639 for (var uintC t = 0; t < tmax; t++) {
641 var uintC i2 = i1 + tmax;
642 // replace (z[i1],z[i2]) by
643 // ((z[i1]+z[i2])/2, (z[i1]-z[i2])/2)
644 // Do the division by 2 later.
651 var const uintC smax = (uintC)1 << (n-1-l);
652 var const uintC tmax = (uintC)1 << (l-1);
653 for (var uintC s = 1; s < smax; s++) {
654 var uintC exp = shuffle(n-1-l,s) << (l-1);
655 for (var uintC t = 0; t < tmax; t++) {
656 var uintC i1 = (s << (l+1)) + t;
657 var uintC i2 = i1 + tmax;
658 var uintC i3 = i2 + tmax;
659 var uintC i4 = i3 + tmax;
660 // replace (z[i1],z[i2],z[i3],z[i4]) by
661 // ((z[i1]+z[i3]+(z[i2]-z[i4])/gam)/2,
662 // (z[i2]+z[i4]-(z[i3]-z[i1])/gam*(gam^2-1))/2,
663 // (z[i2]-z[i4])/(gam*2),
664 // (z[i3]-z[i1])/(gam*2))
665 // Do the division by 2 later.
670 sum13 = z[i1] + z[i3];
671 sum24 = z[i2] + z[i4];
672 tmp4 = (z[i3] - z[i1]) * w[exp].gam.f;
673 tmp3 = (z[i2] - z[i4]) * w[exp].gam.f;
674 z[i1] = sum13 + tmp3;
675 z[i2] = sum24 - tmp4*w[exp].gam.e;
681 // Do all divisions by 2 now.
683 var fftr_real f = (fftr_real)2.0 / (fftr_real)N; // 2^-(n-1)
684 for (i = 0; i < N; i++)
690 // For a given k >= 2, the maximum l is determined by
691 // 2*l < m - 1/2 - 2*k - log_2(12*k-15) - (1 if e=1.0, 0 if e=0.5).
692 // This is a decreasing function of k.
694 (int)((fftr_real_mant_bits \
696 - ((k)<=2 ? 4 : (k)<=3 ? 5 : (k)<=5 ? 6 : (k)<=8 ? 7 : (k)<=16 ? 8 : (k)<=31 ? 9 : 10) \
697 - (fftr_real_rounds == rounds_to_nearest ? 0 : 1)) \
699 static int max_l_table[32+1] =
700 { 0, 0, max_l(2), max_l(3), max_l(4), max_l(5), max_l(6),
701 max_l(7), max_l(8), max_l(9), max_l(10), max_l(11), max_l(12), max_l(13),
702 max_l(14), max_l(15), max_l(16), max_l(17), max_l(18), max_l(19), max_l(20),
703 max_l(21), max_l(22), max_l(23), max_l(24), max_l(25), max_l(26), max_l(27),
704 max_l(28), max_l(29), max_l(30), max_l(31), max_l(32)
707 // Split len uintD's below sourceptr into chunks of l bits, thus filling
708 // N real numbers at x.
709 static void fill_factor (uintC N, fftr_real* x, uintL l,
710 const uintD* sourceptr, uintC len)
713 if (max_l(2) > intDsize && l > intDsize) {
715 if (max_l(2) > 64 && l > 64) {
716 throw runtime_exception("FFT problem: l > 64 not supported by pow2_table");
718 var fftr_real carry = 0;
719 var sintL carrybits = 0; // number of bits in carry (>=0, <l)
722 var uintD digit = lsprefnext(sourceptr);
723 if (carrybits+intDsize >= l) {
724 x[i] = carry + (fftr_real)(digit & bitm(l-carrybits)) * fftr_pow2_table[carrybits];
726 carry = (l-carrybits == intDsize ? (fftr_real)0 : (fftr_real)(digit >> (l-carrybits)));
727 carrybits = carrybits+intDsize-l;
729 carry = carry + (fftr_real)digit * fftr_pow2_table[carrybits];
730 carrybits = carrybits+intDsize;
739 throw runtime_exception();
740 } else if (max_l(2) >= intDsize && l == intDsize) {
743 throw runtime_exception();
744 for (i = 0; i < len; i++) {
745 var uintD digit = lsprefnext(sourceptr);
746 x[i] = (fftr_real)digit;
750 var const uintD l_mask = bit(l)-1;
752 var sintL carrybits = 0; // number of bits in carry (>=0, <intDsize)
753 for (i = 0; i < N; i++) {
754 if (carrybits >= l) {
755 x[i] = (fftr_real)(carry & l_mask);
762 var uintD digit = lsprefnext(sourceptr);
763 x[i] = (fftr_real)((carry | (digit << carrybits)) & l_mask);
764 carry = digit >> (l-carrybits);
765 carrybits = intDsize - (l-carrybits);
768 while (carrybits > 0) {
770 throw runtime_exception();
771 x[i] = (fftr_real)(carry & l_mask);
777 throw runtime_exception();
783 // Given a not too large floating point number, round it to the nearest integer.
784 static inline fftr_real fftr_fround (fftr_real x)
787 #if (fftr_real_rounds == rounds_to_nearest)
788 (x + (fftr_pow2(fftr_real_mant_bits-1)+fftr_pow2(fftr_real_mant_bits-2)))
789 - (fftr_pow2(fftr_real_mant_bits-1)+fftr_pow2(fftr_real_mant_bits-2));
790 #elif (fftr_real_rounds == rounds_to_infinity)
791 (fftr_pow2(fftr_real_mant_bits-1)+fftr_pow2(fftr_real_mant_bits-2))
792 - ((fftr_pow2(fftr_real_mant_bits-1)+fftr_pow2(fftr_real_mant_bits-2))
793 - (x + (fftr_real)0.5));
794 #else // rounds_to_zero, rounds_to_minus_infinity
795 ((x + (fftr_real)0.5)
796 + (fftr_pow2(fftr_real_mant_bits-1)+fftr_pow2(fftr_real_mant_bits-2)))
797 - (fftr_pow2(fftr_real_mant_bits-1)+fftr_pow2(fftr_real_mant_bits-2));
801 // Given a not too large floating point number, round it down.
802 static inline fftr_real fftr_ffloor (fftr_real x)
804 #if (fftr_real_rounds == rounds_to_nearest)
806 (x + (fftr_pow2(fftr_real_mant_bits-1)+fftr_pow2(fftr_real_mant_bits-2)))
807 - (fftr_pow2(fftr_real_mant_bits-1)+fftr_pow2(fftr_real_mant_bits-2));
811 return y - (fftr_real)1.0;
812 #elif (fftr_real_rounds == rounds_to_infinity)
814 (fftr_pow2(fftr_real_mant_bits-1)+fftr_pow2(fftr_real_mant_bits-2))
815 - ((fftr_pow2(fftr_real_mant_bits-1)+fftr_pow2(fftr_real_mant_bits-2))
817 #else // rounds_to_zero, rounds_to_minus_infinity
819 (x + (fftr_pow2(fftr_real_mant_bits-1)+fftr_pow2(fftr_real_mant_bits-2)))
820 - (fftr_pow2(fftr_real_mant_bits-1)+fftr_pow2(fftr_real_mant_bits-2));
824 // Combine the N real numbers at z into uintD's below destptr.
825 // The z[i] are known to be approximately integers >= 0, < N*2^(2*l).
826 // Assumes room for floor(N*l/intDsize)+(1+ceiling((n+2*l)/intDsize)) uintD's
827 // below destptr. Fills len digits and returns (destptr lspop len).
828 static uintD* unfill_product (uintL n, uintC N, // N = 2^n
829 const fftr_real * z, uintL l,
833 if (n + 2*l <= intDsize) {
834 // 2-digit carry is sufficient, l < intDsize
835 var uintD carry0 = 0;
836 var uintD carry1 = 0;
837 var uintL shift = 0; // shift next digit before adding it to the carry, >=0, <intDsize
838 for (i = 0; i < N; i++) {
839 // Fetch z[i] and round it to the nearest uintD.
840 var uintD digit = (uintD)(z[i] + (fftr_real)0.5);
842 carry1 += digit >> (intDsize-shift);
843 digit = digit << shift;
845 if ((carry0 += digit) < digit)
848 if (shift >= intDsize) {
849 lsprefnext(destptr) = carry0;
855 lsprefnext(destptr) = carry0;
856 lsprefnext(destptr) = carry1;
857 } else if (n + 2*l <= 2*intDsize) {
858 // 3-digit carry is sufficient, l < intDsize
860 var uintDD carry0 = 0;
861 var uintD carry1 = 0;
862 var uintL shift = 0; // shift next digit before adding it to the carry, >=0, <intDsize
863 for (i = 0; i < N; i++) {
864 // Fetch z[i] and round it to the nearest uintDD.
865 var uintDD digit = (uintDD)(z[i] + (fftr_real)0.5);
867 carry1 += (uintD)(digit >> (2*intDsize-shift));
868 digit = digit << shift;
870 if ((carry0 += digit) < digit)
873 if (shift >= intDsize) {
874 lsprefnext(destptr) = lowD(carry0);
875 carry0 = highlowDD(carry1,highD(carry0));
880 lsprefnext(destptr) = lowD(carry0);
881 lsprefnext(destptr) = highD(carry0);
882 lsprefnext(destptr) = carry1;
884 var uintD carry0 = 0;
885 var uintD carry1 = 0;
886 var uintD carry2 = 0;
887 var uintL shift = 0; // shift next digit before adding it to the carry, >=0, <intDsize
888 for (i = 0; i < N; i++) {
889 // Fetch z[i] and round it to the nearest uintDD.
890 var fftr_real zi = z[i] + (fftr_real)0.5;
891 var const fftr_real pow2_intDsize = fftr_pow2(intDsize);
892 var const fftr_real invpow2_intDsize = (fftr_real)1.0/pow2_intDsize;
893 var uintD digit1 = (uintD)(zi * invpow2_intDsize);
894 var uintD digit0 = (uintD)(zi - (fftr_real)digit1 * pow2_intDsize);
896 carry2 += digit1 >> (intDsize-shift);
897 digit1 = (digit1 << shift) | (digit0 >> (intDsize-shift));
898 digit0 = digit0 << shift;
900 if ((carry0 += digit0) < digit0)
901 if ((carry1 += 1) == 0)
903 if ((carry1 += digit1) < digit1)
906 if (shift >= intDsize) {
907 lsprefnext(destptr) = carry0;
914 lsprefnext(destptr) = carry0;
915 lsprefnext(destptr) = carry1;
916 lsprefnext(destptr) = carry2;
919 // 1-digit+1-float carry is sufficient
920 var uintD carry0 = 0;
921 var fftr_real carry1 = 0;
922 var uintL shift = 0; // shift next digit before adding it to the carry, >=0, <intDsize
923 for (i = 0; i < N; i++) {
924 // Fetch z[i] and round it to the nearest integer.
925 var fftr_real digit = fftr_fround(z[i]);
927 if (!(digit >= (fftr_real)0
928 && z[i] > digit - (fftr_real)0.5
929 && z[i] < digit + (fftr_real)0.5))
930 throw runtime_exception();
933 digit = digit * fftr_pow2_table[shift];
934 var fftr_real digit1 = fftr_ffloor(digit*((fftr_real)1.0/fftr_pow2(intDsize)));
935 var uintD digit0 = (uintD)(digit - digit1*fftr_pow2(intDsize));
937 if ((carry0 += digit0) < digit0)
938 carry1 += (fftr_real)1.0;
940 while (shift >= intDsize) {
941 lsprefnext(destptr) = carry0;
942 var fftr_real tmp = fftr_ffloor(carry1*((fftr_real)1.0/fftr_pow2(intDsize)));
943 carry0 = (uintD)(carry1 - tmp*fftr_pow2(intDsize));
948 if (carry0 > 0 || carry1 > (fftr_real)0.0) {
949 lsprefnext(destptr) = carry0;
950 while (carry1 > (fftr_real)0.0) {
951 var fftr_real tmp = fftr_ffloor(carry1*((fftr_real)1.0/fftr_pow2(intDsize)));
952 lsprefnext(destptr) = (uintD)(carry1 - tmp*fftr_pow2(intDsize));
960 static inline void mulu_fftr_nocheck (const uintD* sourceptr1, uintC len1,
961 const uintD* sourceptr2, uintC len2,
963 // Es ist 2 <= len1 <= len2.
965 // We have to find parameters l and k such that
966 // ceiling(len1*intDsize/l) + ceiling(len2*intDsize/l) - 1 <= 2^k,
967 // and (12*k-15)*e*2^(2*l+2*k+1/2-m) < 1/2.
968 // Try primarily to minimize k. Minimizing l buys you nothing.
970 // Computing k: If len1 and len2 differ much, we'll split source2 -
971 // hence for the moment just substitute len1 for len2.
973 // First approximation of k: A necessary condition for
974 // 2*ceiling(len1*intDsize/l) - 1 <= 2^k
975 // is 2*len1*intDsize/l_max - 1 <= 2^k.
977 var const int l = max_l(2);
978 var uintC lhs = 2*ceiling(len1*intDsize,l) - 1; // >=1
982 integerlengthC(lhs-1, k=); // k>=2
984 // Try whether this k is ok or whether we have to increase k.
986 if (k >= sizeof(max_l_table)/sizeof(max_l_table[0])
987 || max_l_table[k] <= 0) {
988 throw runtime_exception("FFT problem: numbers too big, floating point precision not sufficient");
990 if (2*ceiling(len1*intDsize,max_l_table[k])-1 <= ((uintC)1 << k))
993 // We could try to reduce l, keeping the same k. But why should we?
994 // Calculate the number of pieces in which source2 will have to be
995 // split. Each of the pieces must satisfy
996 // ceiling(len1*intDsize/l) + ceiling(len2*intDsize/l) - 1 <= 2^k,
998 // Try once with k, once with k+1. Compare them.
1000 var uintC remaining_k = ((uintC)1 << k) + 1 - ceiling(len1*intDsize,max_l_table[k]);
1001 var uintC max_piecelen_k = floor(remaining_k*max_l_table[k],intDsize);
1002 var uintC numpieces_k = ceiling(len2,max_piecelen_k);
1003 var uintC remaining_k1 = ((uintC)1 << (k+1)) + 1 - ceiling(len1*intDsize,max_l_table[k+1]);
1004 var uintC max_piecelen_k1 = floor(remaining_k1*max_l_table[k+1],intDsize);
1005 var uintC numpieces_k1 = ceiling(len2,max_piecelen_k1);
1006 if (numpieces_k <= 2*numpieces_k1) {
1008 len2p = max_piecelen_k;
1012 len2p = max_piecelen_k1;
1015 var const uintL l = max_l_table[k];
1016 var const uintL n = k;
1017 var const uintC N = (uintC)1 << n;
1019 var fftr_real* const x = cl_alloc_array(fftr_real,N);
1020 var fftr_real* const y = cl_alloc_array(fftr_real,N);
1022 var fftr_real* const z = cl_alloc_array(fftr_real,N);
1024 var fftr_real* const z = x; // put z in place of x - saves memory
1026 var uintD* const tmpprod1 = cl_alloc_array(uintD,len1+1);
1027 var uintC tmpprod_len = floor(l<<n,intDsize)+6;
1028 var uintD* const tmpprod = cl_alloc_array(uintD,tmpprod_len);
1029 var uintC destlen = len1+len2;
1030 clear_loop_lsp(destptr,destlen);
1036 var uintD* tmpptr = arrayLSDptr(tmpprod1,len1+1);
1037 mulu_loop_lsp(lspref(sourceptr2,0),sourceptr1,tmpptr,len1);
1038 if (addto_loop_lsp(tmpptr,destptr,len1+1))
1039 if (inc_loop_lsp(destptr lspop (len1+1),destlen-(len1+1)))
1040 throw runtime_exception();
1042 var bool squaring = ((sourceptr1 == sourceptr2) && (len1 == len2p));
1044 fill_factor(N,x,l,sourceptr1,len1);
1047 fill_factor(N,y,l,sourceptr2,len2p);
1050 fftr_convolution(n,N, &x[0], &y[0], &z[0]);
1052 fftr_convolution(n,N, &x[0], &x[0], &z[0]);
1056 var fftr_real re_lo_limit = (fftr_real)(-0.5);
1057 var fftr_real re_hi_limit = (fftr_real)N * fftr_pow2_table[l] * fftr_pow2_table[l] + (fftr_real)0.5;
1058 for (var uintC i = 0; i < N; i++)
1059 if (!(z[i] > re_lo_limit
1060 && z[i] < re_hi_limit))
1061 throw runtime_exception();
1064 var uintD* tmpLSDptr = arrayLSDptr(tmpprod,tmpprod_len);
1065 var uintD* tmpMSDptr = unfill_product(n,N,z,l,tmpLSDptr);
1067 #if CL_DS_BIG_ENDIAN_P
1068 tmpLSDptr - tmpMSDptr;
1070 tmpMSDptr - tmpLSDptr;
1072 if (tmplen > tmpprod_len)
1073 throw runtime_exception();
1074 // Add result to destptr[-destlen..-1]:
1075 if (tmplen > destlen) {
1076 if (test_loop_msp(tmpMSDptr,tmplen-destlen))
1077 throw runtime_exception();
1080 if (addto_loop_lsp(tmpLSDptr,destptr,tmplen))
1081 if (inc_loop_lsp(destptr lspop tmplen,destlen-tmplen))
1082 throw runtime_exception();
1085 destptr = destptr lspop len2p;
1087 sourceptr2 = sourceptr2 lspop len2p;
1095 // Compute a checksum: number mod (2^intDsize-1).
1096 static uintD compute_checksum (const uintD* sourceptr, uintC len)
1098 var uintD tmp = ~(uintD)0; // -1-(sum mod 2^intDsize-1), always >0
1100 var uintD digit = lsprefnext(sourceptr);
1102 tmp -= digit; // subtract digit
1104 tmp -= digit+1; // subtract digit-(2^intDsize-1)
1105 } while (--len > 0);
1109 // Multiply two checksums modulo (2^intDsize-1).
1110 static inline uintD multiply_checksum (uintD checksum1, uintD checksum2)
1115 var uintDD cksum = muluD(checksum1,checksum2);
1116 cksum_hi = highD(cksum); checksum = lowD(cksum);
1118 muluD(checksum1,checksum2, cksum_hi =, checksum =);
1120 if ((checksum += cksum_hi) + 1 <= cksum_hi)
1127 static void mulu_fftr (const uintD* sourceptr1, uintC len1,
1128 const uintD* sourceptr2, uintC len2,
1131 // Compute checksums of the arguments and multiply them.
1132 var uintD checksum1 = compute_checksum(sourceptr1,len1);
1133 var uintD checksum2 = compute_checksum(sourceptr2,len2);
1134 var uintD checksum = multiply_checksum(checksum1,checksum2);
1135 mulu_fftr_nocheck(sourceptr1,len1,sourceptr2,len2,destptr);
1136 if (!(checksum == compute_checksum(destptr,len1+len2))) {
1137 throw runtime_exception("FFT problem: checksum error.");