cl_DS_mul_fftc.h

   1 // Fast integer multiplication using FFT over the complex numbers.
   2 // [Donald Ervin Knuth: The Art of Computer Programming, Vol. II:
   3 //  Seminumerical Algorithms, second edition. Section 4.3.3, p. 290-294.]
   4 // Bruno Haible 6.5.1996, 24.-25.8.1996
   5
   6 // FFT in the complex domain has the drawback that it needs careful round-off
   7 // error analysis. But for CPUs with good floating-point performance it might
   8 // nevertheless be better than FFT mod Z/mZ.
   9
  10 // It is important to have precise round-off error estimates. Although
  11 // (by laws of statistics) in the average the actual round-off error makes up
  12 // only half of the bits provided for round-off protection, we cannot rely
  13 // on this average behaviour, but have to produce correct results.
  14 //
  15 // Knuth's formula (42), p. 294, says:
  16 // If we want to multiply l-bit words using an FFT(2^k), our floating point
  17 // numbers shall have m >= 2(k+l) + k + log_2 k + 3.5 mantissa bits.
  18 //
  19 // Here is a more careful analysis, using absolute error estimates.
  20 //
  21 // 1. We assume floating point numbers with radix 2, with the properties:
  22 //    (i)   Multiplication with 2^n and 2^-n is exact.
  23 //    (ii)  Negation x -> -x is exact.
  24 //    (iii) Addition: When adding x and y, with |x| <= 2^a, |y| <= 2^a,
  25 //          the result |x+y| <= 2^(a+1) has an error <= e*2^(a+1-m).
  26 //    (iv)  Multiplication: When multiplying x and y, with |x| <= 2^a,
  27 //          |y| <= 2^b, the result |x*y| <= 2^(a+b) has an error <= e*2^(a+b-m).
  28 //    Here e = 1 for a truncating arithmetic, but e = 1/2 for a rounding
  29 //    arithmetic like IEEE single and double floats.
  30 // 2. Let's introduce some notation: err(x) means |x'-x| where x is the
  31 //    exact mathematical value and x' is its representation in the machine.
  32 // 3. From 1. we get for real numbers x,y:
  33 //    (i)   err(2^n*x) = 2^n * err(x),
  34 //    (ii)  err(-x) = err(x),
  35 //    (iii) |x| <= 2^a, |y| <= 2^a, then
  36 //          err(x+y) <= err(x) + err(y) + e*2^(a+1-m),
  37 //          [or ....    ............... + e*2^(a-m) if |x+y| <= 2^a].
  38 //    (iv)  |x| <= 2^a, |y| <= 2^b, then
  39 //          err(x*y) <= 2^a * err(y) + 2^b * err(x) + e*2^(a+b-m).
  40 // 4. Our complex arithmetic will be based on the formulas:
  41 //    (i)   2^n*(x+iy) = (2^n*x)+i(2^n*y)
  42 //    (ii)  -(x+iy) = (-x)+i(-y)
  43 //    (iii) (x+iy)+(u+iv) = (x+u)+i(y+v)
  44 //    (iv)  (x+iy)*(u+iv) = (x*u-y*v)+i(x*v+y*u)
  45 //    The notation err(z) means |z'-z|, as above, with |.| being the usual
  46 //    absolute value on complex numbers (_not_ the L^1 norm).
  47 // 5. From 3. and 4. we get for complex numbers x,y:
  48 //    (i)   err(2^n*x) = 2^n * err(x),
  49 //    (ii)  err(-x) = err(x),
  50 //    (iii) |x| <= 2^a, |y| <= 2^a, then
  51 //          err(x+y) <= err(x) + err(y) + e*2^(a+3/2-m),
  52 //    (iv)  |x| <= 2^a, |y| <= 2^b, then
  53 //          err(x*y) <= 2^a * err(y) + 2^b * err(x) + 3*e*2^(a+b+1/2-m).
  54 // 6. We start out with precomputed roots of unity:
  55 //          |exp(2 pi i/2^n)| <= 1,
  56 //          err(exp(2 pi i/2^n)) <= e*2^(1/2-m), (even err(..)=0 for n=0,1,2),
  57 //    and compute  exp(2 pi i * j/2^k)  according to the binary digits of j.
  58 //    This way, each root of unity will be a product of at most k precomputed
  59 //    roots. If (j mod 2^(k-2)) has n bits, then  exp(2 pi i * j/2^k)  will
  60 //    be computed using n factors, i.e. n-1 complex multiplications, and by
  61 //    5.iv. we'll have
  62 //          err(exp(2 pi i * j/2^k)) <= n*e*2^(1/2-m) + max(n-1,0)*3*e*2^(1/2-m)
  63 //                                    = max(4*n-3,0)*e*2^(1/2-m).
  64 //    Hence the maximum roots-of-unity error is (set n=k-2)
  65 //          err(w^j) <= (4*k-11)*e*2^(1/2-m),
  66 //    and the average roots-of-unity error is (set n=(k-2)/2)
  67 //          < 2*(k-2)*e*2^(1/2-m).
  68 // 7. Now we start the FFT.
  69 //    Before the first step, x_i are integral, |x_i| < 2^l and err(x_i) = 0.
  70 //    After the first butterfly, which replaces (x(i1),x(i2)) by
  71 //    (x(i1) + x(i2), x(i1) - x(i2)), we have |x_i| < 2^(l+1) and err(x_i) = 0.
  72 //    Then, for each of the remaining k-1 steps, a butterfly replaces
  73 //    (x(i1),x(i2)) by (x(i1) + w^exp*x(i2), x(i1) - w^exp*x(i2)). Thus,
  74 //    after n steps we have |x_i| < 2^(l+n) and err(x_i) <= E_i where
  75 //          E_0 = 0,
  76 //          E_1 = 0,
  77 //          E_(n+1) = (E_n + 2^(l+n)*(4*k-11)*e*2^(1/2-m) + 3*e*2^(l+n+1/2-m))
  78 //                    + E_n + e*2^(l+n+3/2-m)
  79 //                  = 2*E_n + (4*k-6)*e*2^(l+n+1/2-m)
  80 //    hence  E_n = (2^n-2)*(4*k-6)*e*2^(l+1/2-m).
  81 //    Setting n = k, we have proved that after the FFT ends, we have
  82 //          |x_i| < 2^(l+k) and err(x_i) <= (4*k-6)*e*2^(l+k+1/2-m).
  83 // 8. The same error analysis holds for the y_i and their FFT. After we
  84 //    multiply z_i := x_i * y_i, we have
  85 //          |z_i| < 2^(2*l+2*k) and err(z_i) <= (8*k-9)*e*2^(2*l+2*k+1/2-m).
  86 // 9. Then an inverse FFT on z_i is done, which is the same as an FFT
  87 //    followed by a permutation and a division by 2^k. After n steps of
  88 //    the FFT, we have |z_i| < 2^(2*l+2*k+n) and err(z_i) <= E_i where
  89 //          E_0 = (8*k-9)*e*2^(2*l+2*k+1/2-m),
  90 //          E_(n+1) = (E_n + 2^(2*l+2*k+n)*(4*k-11)*e*2^(1/2-m)
  91 //                         + 3*e*2^(2*l+2*k+n+1/2-m))
  92 //                    + E_n + e*2^(2*l+2*k+n+3/2-m)
  93 //                  = 2*E_n + (4*k-6)*e*2^(2*l+2*k+n+1/2-m)
  94 //    hence  E_n = 2^n*(8*k-9)*e*2^(2*l+2*k+1/2-m)
  95 //                 + (2^n-1)*(4*k-6)*e*2^(2*l+2*k+1/2-m).
  96 //    So, after the FFT, we have (set n=k) |z_i| < 2^(2*l+3*k) and
  97 //          err(z_i) <= (12*k-15)*e*2^(2*l+3*k+1/2-m).
  98 //    Permutation doesn't change the estimates. After division by 2^k, we get
  99 //    |z_i| < 2^(2*l+2*k) and
 100 //          err(z_i) <= (12*k-15)*e*2^(2*l+2*k+1/2-m).
 101 // 10. When converting the z_i back to integers, we know that z_i should be
 102 //     real, integral, and |z_i| < 2^(2*l+k). We can only guarantee that we
 103 //     can find the integral z_i from the floating-point computation if
 104 //          (12*k-15)*e*2^(2*l+2*k+1/2-m) < 1/2.
 105 // 11. Assuming e = 1/2 and m = 53 (typical values for IEEE double arithmetic),
 106 //     we get the constraint  2*l < m - 1/2 - 2*k - log_2(12*k-15).
 107 //          k = 2    l <= 22
 108 //          k = 3    l <= 21
 109 //          k = 4    l <= 19
 110 //          k = 5    l <= 18
 111 //          k = 6    l <= 17
 112 //          k = 7    l <= 16
 113 //          k = 8    l <= 15
 114 //          k = 9    l <= 13
 115 //          k = 10   l <= 12
 116 //          k = 11   l <= 11
 117 //          k = 12   l <= 10
 118 //          k = 13   l <= 9
 119 //          k = 14   l <= 8
 120 //          k = 15   l <= 7
 121 //          k = 16   l <= 6
 122 //          k = 17   l <= 5
 123 //          k = 18   l <= 4
 124 //          k = 19   l <= 3
 125 //          k = 20   l <= 2
 126 //     Assuming e = 1/2 and m = 64 ("long double" arithmetic on i387/i486/i586),
 127 //     we get the constraint  2*l < m - 1/2 - 2*k - log_2(12*k-15).
 128 //          k = 2    l <= 28
 129 //          k = 3    l <= 26
 130 //          k = 4    l <= 25
 131 //          k = 5    l <= 24
 132 //          k = 6    l <= 22
 133 //          k = 7    l <= 21
 134 //          k = 8    l <= 20
 135 //          k = 9    l <= 19
 136 //          k = 10   l <= 18
 137 //          k = 11   l <= 17
 138 //          k = 12   l <= 16
 139 //          k = 13   l <= 15
 140 //          k = 14   l <= 14
 141 //          k = 15   l <= 13
 142 //          k = 16   l <= 12
 143 //          k = 17   l <= 10
 144 //          k = 18   l <= 9
 145 //          k = 19   l <= 8
 146 //          k = 20   l <= 7
 147 //          k = 21   l <= 6
 148 //          k = 22   l <= 5
 149 //          k = 23   l <= 4
 150 //          k = 24   l <= 3
 151 //          k = 25   l <= 2
 152
 153
 154 #if !(intDsize==32)
 155 #error "complex fft implemented only for intDsize==32"
 156 #endif
 157
 158
 159 #include "cln/floatparam.h"
 160 #include "cln/exception.h"
 161
 162 #if defined(HAVE_LONGDOUBLE) && (long_double_mant_bits > double_mant_bits) && (defined(__i386__) || defined(__m68k__) || (defined(__sparc__) && 0))
 163 // Only these CPUs have fast "long double"s in hardware.
 164 // On SPARC, "long double"s are emulated in software and don't work.
 165 typedef long double fftc_real;
 166 #define fftc_real_mant_bits long_double_mant_bits
 167 #define fftc_real_rounds long_double_rounds
 168 #else
 169 typedef double fftc_real;
 170 #define fftc_real_mant_bits double_mant_bits
 171 #define fftc_real_rounds double_rounds
 172 #endif
 173
 174 typedef struct fftc_complex
 175   {
 176     fftc_real re;
 177     fftc_real im;
 178   }
 179   fftc_complex;
 180
 181 static const fftc_complex fftc_roots_of_1 [32+1] =
 182   // roots_of_1[n] is a (2^n)th root of unity in C.
 183   // Also roots_of_1[n-1] = roots_of_1[n]^2.
 184   // For simplicity we choose  roots_of_1[n] = exp(2 pi i/2^n).
 185   {
 186    #if (fftc_real_mant_bits == double_mant_bits)
 187     // These values have 64 bit precision.
 188     { 1.0,                    0.0 },
 189     { -1.0,                   0.0 },
 190     { 0.0,                    1.0 },
 191     { 0.7071067811865475244,  0.7071067811865475244 },
 192     { 0.9238795325112867561,  0.38268343236508977172 },
 193     { 0.9807852804032304491,  0.19509032201612826784 },
 194     { 0.99518472667219688623, 0.098017140329560601996 },
 195     { 0.9987954562051723927,  0.049067674327418014254 },
 196     { 0.9996988186962042201,  0.024541228522912288032 },
 197     { 0.99992470183914454094, 0.0122715382857199260795 },
 198     { 0.99998117528260114264, 0.0061358846491544753597 },
 199     { 0.9999952938095761715,  0.00306795676296597627 },
 200     { 0.99999882345170190993, 0.0015339801862847656123 },
 201     { 0.99999970586288221914, 7.6699031874270452695e-4 },
 202     { 0.99999992646571785114, 3.8349518757139558907e-4 },
 203     { 0.9999999816164292938,  1.9174759731070330744e-4 },
 204     { 0.9999999954041073129,  9.5873799095977345874e-5 },
 205     { 0.99999999885102682754, 4.793689960306688455e-5 },
 206     { 0.99999999971275670683, 2.3968449808418218729e-5 },
 207     { 0.9999999999281891767,  1.1984224905069706422e-5 },
 208     { 0.99999999998204729416, 5.9921124526424278428e-6 },
 209     { 0.99999999999551182357, 2.9960562263346607504e-6 },
 210     { 0.99999999999887795586, 1.4980281131690112288e-6 },
 211     { 0.999999999999719489,   7.4901405658471572114e-7 },
 212     { 0.99999999999992987223, 3.7450702829238412391e-7 },
 213     { 0.99999999999998246807, 1.8725351414619534487e-7 },
 214     { 0.999999999999995617,   9.36267570730980828e-8 },
 215     { 0.99999999999999890425, 4.6813378536549092695e-8 },
 216     { 0.9999999999999997261,  2.340668926827455276e-8 },
 217     { 0.99999999999999993153, 1.1703344634137277181e-8 },
 218     { 0.99999999999999998287, 5.8516723170686386908e-9 },
 219     { 0.9999999999999999957,  2.925836158534319358e-9 },
 220     { 0.9999999999999999989,  1.4629180792671596806e-9 }
 221    #else (fftc_real_mant_bits > double_mant_bits)
 222     // These values have 128 bit precision.
 223     { 1.0L,                                       0.0L },
 224     { -1.0L,                                      0.0L },
 225     { 0.0L,                                       1.0L },
 226     { 0.707106781186547524400844362104849039284L, 0.707106781186547524400844362104849039284L },
 227     { 0.923879532511286756128183189396788286823L, 0.38268343236508977172845998403039886676L },
 228     { 0.980785280403230449126182236134239036975L, 0.195090322016128267848284868477022240928L },
 229     { 0.995184726672196886244836953109479921574L, 0.098017140329560601994195563888641845861L },
 230     { 0.998795456205172392714771604759100694444L, 0.0490676743274180142549549769426826583147L },
 231     { 0.99969881869620422011576564966617219685L,  0.0245412285229122880317345294592829250654L },
 232     { 0.99992470183914454092164649119638322435L,  0.01227153828571992607940826195100321214037L },
 233     { 0.999981175282601142656990437728567716173L, 0.00613588464915447535964023459037258091705L },
 234     { 0.999995293809576171511580125700119899554L, 0.00306795676296597627014536549091984251894L },
 235     { 0.99999882345170190992902571017152601905L,  0.001533980186284765612303697150264079079954L },
 236     { 0.999999705862882219160228217738765677117L, 7.66990318742704526938568357948576643142e-4L },
 237     { 0.99999992646571785114473148070738785695L,  3.83495187571395589072461681181381263396e-4L },
 238     { 0.999999981616429293808346915402909714504L, 1.91747597310703307439909561989000933469e-4L },
 239     { 0.99999999540410731289097193313960614896L,  9.58737990959773458705172109764763511872e-5L },
 240     { 0.9999999988510268275626733077945541084L,   4.79368996030668845490039904946588727468e-5L },
 241     { 0.99999999971275670684941397221864177609L,  2.39684498084182187291865771650218200947e-5L },
 242     { 0.999999999928189176709775095883850490262L, 1.198422490506970642152156159698898480473e-5L },
 243     { 0.99999999998204729417728262414778410738L,  5.99211245264242784287971180889086172999e-6L },
 244     { 0.99999999999551182354431058417299732444L,  2.99605622633466075045481280835705981183e-6L },
 245     { 0.999999999998877955886077016551752536504L, 1.49802811316901122885427884615536112069e-6L },
 246     { 0.999999999999719488971519214794719584451L, 7.49014056584715721130498566730655637157e-7L },
 247     { 0.99999999999992987224287980123972873676L,  3.74507028292384123903169179084633177398e-7L },
 248     { 0.99999999999998246806071995015624773673L,  1.8725351414619534486882457659356361712e-7L },
 249     { 0.999999999999995617015179987529456656217L, 9.3626757073098082799067286680885620193e-8L },
 250     { 0.999999999999998904253794996881763834182L, 4.68133785365490926951155181385400969594e-8L },
 251     { 0.99999999999999972606344874922040343793L,  2.34066892682745527595054934190348440379e-8L },
 252     { 0.999999999999999931515862187305098514444L, 1.170334463413727718124621350323810379807e-8L },
 253     { 0.999999999999999982878965546826274482047L, 5.8516723170686386908097901008341396944e-9L },
 254     { 0.999999999999999995719741386706568611352L, 2.92583615853431935792823046906895590202e-9L },
 255     { 0.999999999999999998929935346676642152265L, 1.46291807926715968052953216186596371037e-9L }
 256    #endif
 257   };
 258
 259 // Define this for (cheap) consistency checks.
 260 //#define DEBUG_FFTC
 261
 262 // Define the algorithm of the backward FFT:
 263 // Either FORWARD (a normal FFT followed by a permutation)
 264 // or     RECIPROOT (an FFT with reciprocal root of unity)
 265 // or     CLEVER (an FFT with reciprocal root of unity but clever computation
 266 //                of the reciprocals).
 267 // Drawback of FORWARD: the permutation pass.
 268 // Drawback of RECIPROOT: need all the powers of the root, not only half of them.
 269 #define FORWARD   42
 270 #define RECIPROOT 43
 271 #define CLEVER    44
 272 #define FFTC_BACKWARD CLEVER
 273
 274 static fftc_real fftc_pow2_table[64] = // table of powers of 2
 275   {
 276                       1.0,
 277                       2.0,
 278                       4.0,
 279                       8.0,
 280                      16.0,
 281                      32.0,
 282                      64.0,
 283                     128.0,
 284                     256.0,
 285                     512.0,
 286                    1024.0,
 287                    2048.0,
 288                    4096.0,
 289                    8192.0,
 290                   16384.0,
 291                   32768.0,
 292                   65536.0,
 293                  131072.0,
 294                  262144.0,
 295                  524288.0,
 296                 1048576.0,
 297                 2097152.0,
 298                 4194304.0,
 299                 8388608.0,
 300                16777216.0,
 301                33554432.0,
 302                67108864.0,
 303               134217728.0,
 304               268435456.0,
 305               536870912.0,
 306              1073741824.0,
 307              2147483648.0,
 308              4294967296.0,
 309              8589934592.0,
 310             17179869184.0,
 311             34359738368.0,
 312             68719476736.0,
 313            137438953472.0,
 314            274877906944.0,
 315            549755813888.0,
 316           1099511627776.0,
 317           2199023255552.0,
 318           4398046511104.0,
 319           8796093022208.0,
 320          17592186044416.0,
 321          35184372088832.0,
 322          70368744177664.0,
 323         140737488355328.0,
 324         281474976710656.0,
 325         562949953421312.0,
 326        1125899906842624.0,
 327        2251799813685248.0,
 328        4503599627370496.0,
 329        9007199254740992.0,
 330       18014398509481984.0,
 331       36028797018963968.0,
 332       72057594037927936.0,
 333      144115188075855872.0,
 334      288230376151711744.0,
 335      576460752303423488.0,
 336     1152921504606846976.0,
 337     2305843009213693952.0,
 338     4611686018427387904.0,
 339     9223372036854775808.0
 340   };
 341
 342 // For a constant expression n (0 <= n < 128), returns 2^n of type fftc_real.
 343 #define fftc_pow2(n)  \
 344   (((n) & 64 ? (fftc_real)18446744073709551616.0 : (fftc_real)1.0)      \
 345    * ((n) & 32 ? (fftc_real)4294967296.0 : (fftc_real)1.0)              \
 346    * ((n) & 16 ? (fftc_real)65536.0 : (fftc_real)1.0)                   \
 347    * ((n) & 8 ? (fftc_real)256.0 : (fftc_real)1.0)                      \
 348    * ((n) & 4 ? (fftc_real)16.0 : (fftc_real)1.0)                       \
 349    * ((n) & 2 ? (fftc_real)4.0 : (fftc_real)1.0)                        \
 350    * ((n) & 1 ? (fftc_real)2.0 : (fftc_real)1.0)                        \
 351   )
 352
 353 // r := a + b
 354 static inline void add (const fftc_complex& a, const fftc_complex& b, fftc_complex& r)
 355 {
 356         r.re = a.re + b.re;
 357         r.im = a.im + b.im;
 358 }
 359
 360 // r := a - b
 361 static inline void sub (const fftc_complex& a, const fftc_complex& b, fftc_complex& r)
 362 {
 363         r.re = a.re - b.re;
 364         r.im = a.im - b.im;
 365 }
 366
 367 // r := a * b
 368 static inline void mul (fftc_real a, const fftc_complex& b, fftc_complex& r)
 369 {
 370         r.re = a * b.re;
 371         r.im = a * b.im;
 372 }
 373
 374 // r := a * b
 375 static inline void mul (const fftc_complex& a, const fftc_complex& b, fftc_complex& r)
 376 {
 377         var fftc_real r_re = a.re * b.re - a.im * b.im;
 378         var fftc_real r_im = a.re * b.im + a.im * b.re;
 379         r.re = r_re;
 380         r.im = r_im;
 381 }
 382
 383 #ifndef _BIT_REVERSE
 384 #define _BIT_REVERSE
 385 // Reverse an n-bit number x. n>0.
 386 static uintC bit_reverse (uintL n, uintC x)
 387 {
 388         var uintC y = 0;
 389         do {
 390                 y <<= 1;
 391                 y |= (x & 1);
 392                 x >>= 1;
 393         } while (!(--n == 0));
 394         return y;
 395 }
 396 #endif
 397
 398 // Compute a complex convolution using FFT: z[0..N-1] := x[0..N-1] * y[0..N-1].
 399 static void fftc_convolution (const uintL n, const uintC N, // N = 2^n
 400                               fftc_complex * x, // N numbers
 401                               fftc_complex * y, // N numbers
 402                               fftc_complex * z  // N numbers result
 403                              )
 404 {
 405         CL_ALLOCA_STACK;
 406         #if (FFTC_BACKWARD == RECIPROOT) || defined(DEBUG_FFTC)
 407         var fftc_complex* const w = cl_alloc_array(fftc_complex,N);
 408         #else
 409         var fftc_complex* const w = cl_alloc_array(fftc_complex,(N>>1)+1);
 410         #endif
 411         var uintC i;
 412         // Initialize w[i] to w^i, w a primitive N-th root of unity.
 413         w[0] = fftc_roots_of_1[0];
 414         #if (FFTC_BACKWARD == RECIPROOT) || defined(DEBUG_FFTC)
 415         {
 416                 var int j;
 417                 for (j = n-1; j>=0; j--) {
 418                         var fftc_complex r_j = fftc_roots_of_1[n-j];
 419                         w[1<<j] = r_j;
 420                         for (i = (2<<j); i < N; i += (2<<j))
 421                                 mul(w[i],r_j, w[i+(1<<j)]);
 422                 }
 423         }
 424         #else
 425         {
 426                 var int j;
 427                 for (j = n-1; j>=0; j--) {
 428                         var fftc_complex r_j = fftc_roots_of_1[n-j];
 429                         w[1<<j] = r_j;
 430                         for (i = (2<<j); i < N>>1; i += (2<<j))
 431                                 mul(w[i],r_j, w[i+(1<<j)]);
 432                 }
 433         }
 434         #endif
 435         #ifdef DEBUG_FFTC
 436         // Check that w is really a primitive N-th root of unity.
 437         {
 438                 var fftc_real epsilon;
 439                 epsilon = (fftc_real)(n>=3 ? 4*n-11 : 0);
 440                 epsilon = epsilon / fftc_pow2(fftc_real_mant_bits);
 441                 if (fftc_real_rounds == rounds_to_nearest)
 442                         epsilon = 0.5*epsilon;
 443                 epsilon = 1.414*epsilon;
 444                 // epsilon = (4*k-11)*e*2^(1/2-m).
 445                 var fftc_real part;
 446                 var fftc_complex& w_N = w[N>>1];
 447                 part = w_N.re - (fftc_real)(-1.0);
 448                 if (part > epsilon || part < -epsilon)
 449                         throw runtime_exception();
 450                 part = w_N.im;
 451                 if (part > epsilon || part < -epsilon)
 452                         throw runtime_exception();
 453         }
 454         {
 455                 var fftc_complex w_N;
 456                 mul(w[N-1],fftc_roots_of_1[n], w_N);
 457                 // Since there was one more multiplication,
 458                 // we have to replace n by (n+1) in the epsilon above.
 459                 var fftc_real epsilon;
 460                 epsilon = (fftc_real)(4*n-7);
 461                 epsilon = epsilon / fftc_pow2(fftc_real_mant_bits);
 462                 if (fftc_real_rounds == rounds_to_nearest)
 463                         epsilon = 0.5*epsilon;
 464                 epsilon = 1.414*epsilon;
 465                 // epsilon = (4*k-7)*e*2^(1/2-m).
 466                 var fftc_real part;
 467                 part = w_N.re - (fftc_real)1.0;
 468                 if (part > epsilon || part < -epsilon)
 469                         throw runtime_exception();
 470                 part = w_N.im;
 471                 if (part > epsilon || part < -epsilon)
 472                         throw runtime_exception();
 473         }
 474         #endif
 475         var bool squaring = (x == y);
 476         // Do an FFT of length N on x.
 477         {
 478                 var sintL l;
 479                 /* l = n-1 */ {
 480                         var const uintC tmax = N>>1; // tmax = 2^(n-1)
 481                         for (var uintC t = 0; t < tmax; t++) {
 482                                 var uintC i1 = t;
 483                                 var uintC i2 = i1 + tmax;
 484                                 // Butterfly: replace (x(i1),x(i2)) by
 485                                 // (x(i1) + x(i2), x(i1) - x(i2)).
 486                                 var fftc_complex tmp;
 487                                 tmp = x[i2];
 488                                 sub(x[i1],tmp, x[i2]);
 489                                 add(x[i1],tmp, x[i1]);
 490                         }
 491                 }
 492                 for (l = n-2; l>=0; l--) {
 493                         var const uintC smax = (uintC)1 << (n-1-l);
 494                         var const uintC tmax = (uintC)1 << l;
 495                         for (var uintC s = 0; s < smax; s++) {
 496                                 var uintC exp = bit_reverse(n-1-l,s) << l;
 497                                 for (var uintC t = 0; t < tmax; t++) {
 498                                         var uintC i1 = (s << (l+1)) + t;
 499                                         var uintC i2 = i1 + tmax;
 500                                         // Butterfly: replace (x(i1),x(i2)) by
 501                                         // (x(i1) + w^exp*x(i2), x(i1) - w^exp*x(i2)).
 502                                         var fftc_complex tmp;
 503                                         mul(x[i2],w[exp], tmp);
 504                                         sub(x[i1],tmp, x[i2]);
 505                                         add(x[i1],tmp, x[i1]);
 506                                 }
 507                         }
 508                 }
 509         }
 510         // Do an FFT of length N on y.
 511         if (!squaring) {
 512                 var sintL l;
 513                 /* l = n-1 */ {
 514                         var uintC const tmax = N>>1; // tmax = 2^(n-1)
 515                         for (var uintC t = 0; t < tmax; t++) {
 516                                 var uintC i1 = t;
 517                                 var uintC i2 = i1 + tmax;
 518                                 // Butterfly: replace (y(i1),y(i2)) by
 519                                 // (y(i1) + y(i2), y(i1) - y(i2)).
 520                                 var fftc_complex tmp;
 521                                 tmp = y[i2];
 522                                 sub(y[i1],tmp, y[i2]);
 523                                 add(y[i1],tmp, y[i1]);
 524                         }
 525                 }
 526                 for (l = n-2; l>=0; l--) {
 527                         var const uintC smax = (uintC)1 << (n-1-l);
 528                         var const uintC tmax = (uintC)1 << l;
 529                         for (var uintC s = 0; s < smax; s++) {
 530                                 var uintC exp = bit_reverse(n-1-l,s) << l;
 531                                 for (var uintC t = 0; t < tmax; t++) {
 532                                         var uintC i1 = (s << (l+1)) + t;
 533                                         var uintC i2 = i1 + tmax;
 534                                         // Butterfly: replace (y(i1),y(i2)) by
 535                                         // (y(i1) + w^exp*y(i2), y(i1) - w^exp*y(i2)).
 536                                         var fftc_complex tmp;
 537                                         mul(y[i2],w[exp], tmp);
 538                                         sub(y[i1],tmp, y[i2]);
 539                                         add(y[i1],tmp, y[i1]);
 540                                 }
 541                         }
 542                 }
 543         }
 544         // Multiply the transformed vectors into z.
 545         for (i = 0; i < N; i++)
 546                 mul(x[i],y[i], z[i]);
 547         // Undo an FFT of length N on z.
 548         {
 549                 var uintL l;
 550                 for (l = 0; l < n-1; l++) {
 551                         var const uintC smax = (uintC)1 << (n-1-l);
 552                         var const uintC tmax = (uintC)1 << l;
 553                         #if FFTC_BACKWARD != CLEVER
 554                         for (var uintC s = 0; s < smax; s++) {
 555                                 var uintC exp = bit_reverse(n-1-l,s) << l;
 556                                 #if FFTC_BACKWARD == RECIPROOT
 557                                 if (exp > 0)
 558                                         exp = N - exp; // negate exp (use w^-1 instead of w)
 559                                 #endif
 560                                 for (var uintC t = 0; t < tmax; t++) {
 561                                         var uintC i1 = (s << (l+1)) + t;
 562                                         var uintC i2 = i1 + tmax;
 563                                         // Inverse Butterfly: replace (z(i1),z(i2)) by
 564                                         // ((z(i1)+z(i2))/2, (z(i1)-z(i2))/(2*w^exp)).
 565                                         // Do the division by 2 later.
 566                                         var fftc_complex sum;
 567                                         var fftc_complex diff;
 568                                         add(z[i1],z[i2], sum);
 569                                         sub(z[i1],z[i2], diff);
 570                                         z[i1] = sum;
 571                                         mul(diff,w[exp], z[i2]);
 572                                 }
 573                         }
 574                         #else // FFTC_BACKWARD == CLEVER: clever handling of negative exponents
 575                         /* s = 0, exp = 0 */ {
 576                                 for (var uintC t = 0; t < tmax; t++) {
 577                                         var uintC i1 = t;
 578                                         var uintC i2 = i1 + tmax;
 579                                         // Inverse Butterfly: replace (z(i1),z(i2)) by
 580                                         // ((z(i1)+z(i2))/2, (z(i1)-z(i2))/(2*w^exp)),
 581                                         // with exp <-- 0.
 582                                         // Do the division by 2 later.
 583                                         var fftc_complex sum;
 584                                         var fftc_complex diff;
 585                                         add(z[i1],z[i2], sum);
 586                                         sub(z[i1],z[i2], diff);
 587                                         z[i1] = sum;
 588                                         z[i2] = diff;
 589                                 }
 590                         }
 591                         for (var uintC s = 1; s < smax; s++) {
 592                                 var uintC exp = bit_reverse(n-1-l,s) << l;
 593                                 exp = (N>>1) - exp; // negate exp (use w^-1 instead of w)
 594                                 for (var uintC t = 0; t < tmax; t++) {
 595                                         var uintC i1 = (s << (l+1)) + t;
 596                                         var uintC i2 = i1 + tmax;
 597                                         // Inverse Butterfly: replace (z(i1),z(i2)) by
 598                                         // ((z(i1)+z(i2))/2, (z(i1)-z(i2))/(2*w^exp)),
 599                                         // with exp <-- (N/2 - exp).
 600                                         // Do the division by 2 later.
 601                                         var fftc_complex sum;
 602                                         var fftc_complex diff;
 603                                         add(z[i1],z[i2], sum);
 604                                         sub(z[i2],z[i1], diff); // note that w^(N/2) = -1
 605                                         z[i1] = sum;
 606                                         mul(diff,w[exp], z[i2]);
 607                                 }
 608                         }
 609                         #endif
 610                 }
 611                 /* l = n-1 */ {
 612                         var const fftc_real pow2 = (fftc_real)1 / (fftc_real)N;
 613                         var const uintC tmax = N>>1; // tmax = 2^(n-1)
 614                         for (var uintC t = 0; t < tmax; t++) {
 615                                 var uintC i1 = t;
 616                                 var uintC i2 = i1 + tmax;
 617                                 // Inverse Butterfly: replace (z(i1),z(i2)) by
 618                                 // ((z(i1)+z(i2))/2, (z(i1)-z(i2))/2).
 619                                 // Do all the divisions by 2 now.
 620                                 var fftc_complex sum;
 621                                 var fftc_complex diff;
 622                                 add(z[i1],z[i2], sum);
 623                                 sub(z[i1],z[i2], diff);
 624                                 mul(pow2,sum, z[i1]);
 625                                 mul(pow2,diff, z[i2]);
 626                         }
 627                 }
 628         }
 629         #if FFTC_BACKWARD == FORWARD
 630         // Swap z[i] and z[N-i] for 0 < i < N/2.
 631         for (i = (N>>1)-1; i > 0; i--) {
 632                 var fftc_complex tmp = z[i];
 633                 z[i] = z[N-i];
 634                 z[N-i] = tmp;
 635         }
 636         #endif
 637 }
 638
 639 // For a given k >= 2, the maximum l is determined by
 640 // 2*l < m - 1/2 - 2*k - log_2(12*k-15) - (1 if e=1.0, 0 if e=0.5).
 641 // This is a decreasing function of k.
 642 #define max_l(k)  \
 643         (int)((fftc_real_mant_bits                      \
 644                - 2*(k)                                  \
 645                - ((k)<=2 ? 4 : (k)<=3 ? 5 : (k)<=5 ? 6 : (k)<=8 ? 7 : (k)<=16 ? 8 : (k)<=31 ? 9 : 10) \
 646                - (fftc_real_rounds == rounds_to_nearest ? 0 : 1)) \
 647               / 2)
 648 static int max_l_table[32+1] =
 649   {         0,         0,  max_l(2),  max_l(3),  max_l(4),  max_l(5),  max_l(6),
 650      max_l(7),  max_l(8),  max_l(9), max_l(10), max_l(11), max_l(12), max_l(13),
 651     max_l(14), max_l(15), max_l(16), max_l(17), max_l(18), max_l(19), max_l(20),
 652     max_l(21), max_l(22), max_l(23), max_l(24), max_l(25), max_l(26), max_l(27),
 653     max_l(28), max_l(29), max_l(30), max_l(31), max_l(32)
 654   };
 655
 656 // Split len uintD's below sourceptr into chunks of l bits, thus filling
 657 // N complex numbers at x.
 658 static void fill_factor (uintC N, fftc_complex* x, uintL l,
 659                          const uintD* sourceptr, uintC len)
 660 {
 661         var uintC i;
 662         if (max_l(2) > intDsize && l > intDsize) {
 663                 // l > intDsize
 664                 if (max_l(2) > 64 && l > 64) {
 665                         throw runtime_exception("FFT problem: l > 64 not supported by pow2_table");
 666                 }
 667                 var fftc_real carry = 0;
 668                 var sintL carrybits = 0; // number of bits in carry (>=0, <l)
 669                 i = 0;
 670                 while (len > 0) {
 671                         var uintD digit = lsprefnext(sourceptr);
 672                         if (carrybits+intDsize >= l) {
 673                                 x[i].re = carry + (fftc_real)(digit & bitm(l-carrybits)) * fftc_pow2_table[carrybits];
 674                                 x[i].im = (fftc_real)0;
 675                                 i++;
 676                                 carry = (l-carrybits == intDsize ? (fftc_real)0 : (fftc_real)(digit >> (l-carrybits)));
 677                                 carrybits = carrybits+intDsize-l;
 678                         } else {
 679                                 carry = carry + (fftc_real)digit * fftc_pow2_table[carrybits];
 680                                 carrybits = carrybits+intDsize;
 681                         }
 682                         len--;
 683                 }
 684                 if (carrybits > 0) {
 685                         x[i].re = carry;
 686                         x[i].im = (fftc_real)0;
 687                         i++;
 688                 }
 689                 if (i > N)
 690                         throw runtime_exception();
 691         } else if (max_l(2) >= intDsize && l == intDsize) {
 692                 // l = intDsize
 693                 if (len > N)
 694                         throw runtime_exception();
 695                 for (i = 0; i < len; i++) {
 696                         var uintD digit = lsprefnext(sourceptr);
 697                         x[i].re = (fftc_real)digit;
 698                         x[i].im = (fftc_real)0;
 699                 }
 700         } else {
 701                 // l < intDsize
 702                 var const uintD l_mask = bit(l)-1;
 703                 var uintD carry = 0;
 704                 var sintL carrybits = 0; // number of bits in carry (>=0, <intDsize)
 705                 for (i = 0; i < N; i++) {
 706                         if (carrybits >= l) {
 707                                 x[i].re = (fftc_real)(carry & l_mask);
 708                                 x[i].im = (fftc_real)0;
 709                                 carry >>= l;
 710                                 carrybits -= l;
 711                         } else {
 712                                 if (len == 0)
 713                                         break;
 714                                 len--;
 715                                 var uintD digit = lsprefnext(sourceptr);
 716                                 x[i].re = (fftc_real)((carry | (digit << carrybits)) & l_mask);
 717                                 x[i].im = (fftc_real)0;
 718                                 carry = digit >> (l-carrybits);
 719                                 carrybits = intDsize - (l-carrybits);
 720                         }
 721                 }
 722                 while (carrybits > 0) {
 723                         if (!(i < N))
 724                                 throw runtime_exception();
 725                         x[i].re = (fftc_real)(carry & l_mask);
 726                         x[i].im = (fftc_real)0;
 727                         carry >>= l;
 728                         carrybits -= l;
 729                         i++;
 730                 }
 731                 if (len > 0)
 732                         throw runtime_exception();
 733         }
 734         for ( ; i < N; i++) {
 735                 x[i].re = (fftc_real)0;
 736                 x[i].im = (fftc_real)0;
 737         }
 738 }
 739
 740 // Given a not too large floating point number, round it to the nearest integer.
 741 static inline fftc_real fftc_fround (fftc_real x)
 742 {
 743         return
 744           #if (fftc_real_rounds == rounds_to_nearest)
 745             (x + (fftc_pow2(fftc_real_mant_bits-1)+fftc_pow2(fftc_real_mant_bits-2)))
 746             - (fftc_pow2(fftc_real_mant_bits-1)+fftc_pow2(fftc_real_mant_bits-2));
 747           #elif (fftc_real_rounds == rounds_to_infinity)
 748             (fftc_pow2(fftc_real_mant_bits-1)+fftc_pow2(fftc_real_mant_bits-2))
 749             - ((fftc_pow2(fftc_real_mant_bits-1)+fftc_pow2(fftc_real_mant_bits-2))
 750                - (x + (fftc_real)0.5));
 751           #else // rounds_to_zero, rounds_to_minus_infinity
 752             ((x + (fftc_real)0.5)
 753              + (fftc_pow2(fftc_real_mant_bits-1)+fftc_pow2(fftc_real_mant_bits-2)))
 754             - (fftc_pow2(fftc_real_mant_bits-1)+fftc_pow2(fftc_real_mant_bits-2));
 755           #endif
 756 }
 757
 758 // Given a not too large floating point number, round it down.
 759 static inline fftc_real fftc_ffloor (fftc_real x)
 760 {
 761         #if (fftc_real_rounds == rounds_to_nearest)
 762           var fftc_real y =
 763             (x + (fftc_pow2(fftc_real_mant_bits-1)+fftc_pow2(fftc_real_mant_bits-2)))
 764             - (fftc_pow2(fftc_real_mant_bits-1)+fftc_pow2(fftc_real_mant_bits-2));
 765           if (y <= x)
 766                 return y;
 767           else
 768                 return y - (fftc_real)1.0;
 769         #elif (fftc_real_rounds == rounds_to_infinity)
 770           return
 771             (fftc_pow2(fftc_real_mant_bits-1)+fftc_pow2(fftc_real_mant_bits-2))
 772             - ((fftc_pow2(fftc_real_mant_bits-1)+fftc_pow2(fftc_real_mant_bits-2))
 773                - x);
 774         #else // rounds_to_zero, rounds_to_minus_infinity
 775           return
 776             (x + (fftc_pow2(fftc_real_mant_bits-1)+fftc_pow2(fftc_real_mant_bits-2)))
 777             - (fftc_pow2(fftc_real_mant_bits-1)+fftc_pow2(fftc_real_mant_bits-2));
 778         #endif
 779 }
 780
 781 // Combine the N complex numbers at z into uintD's below destptr.
 782 // The z[i] are known to be approximately integers >= 0, < N*2^(2*l).
 783 // Assumes room for floor(N*l/intDsize)+(1+ceiling((n+2*l)/intDsize)) uintD's
 784 // below destptr. Fills len digits and returns (destptr lspop len).
 785 static uintD* unfill_product (uintL n, uintC N, // N = 2^n
 786                               const fftc_complex * z, uintL l,
 787                               uintD* destptr)
 788 {
 789         var uintC i;
 790         if (n + 2*l <= intDsize) {
 791                 // 2-digit carry is sufficient, l < intDsize
 792                 var uintD carry0 = 0;
 793                 var uintD carry1 = 0;
 794                 var uintL shift = 0; // shift next digit before adding it to the carry, >=0, <intDsize
 795                 for (i = 0; i < N; i++) {
 796                         // Fetch z[i].re and round it to the nearest uintD.
 797                         var uintD digit = (uintD)(z[i].re + (fftc_real)0.5);
 798                         if (shift > 0) {
 799                                 carry1 += digit >> (intDsize-shift);
 800                                 digit = digit << shift;
 801                         }
 802                         if ((carry0 += digit) < digit)
 803                                 carry1 += 1;
 804                         shift += l;
 805                         if (shift >= intDsize) {
 806                                 lsprefnext(destptr) = carry0;
 807                                 carry0 = carry1;
 808                                 carry1 = 0;
 809                                 shift -= intDsize;
 810                         }
 811                 }
 812                 lsprefnext(destptr) = carry0;
 813                 lsprefnext(destptr) = carry1;
 814         } else if (n + 2*l <= 2*intDsize) {
 815                 // 3-digit carry is sufficient, l < intDsize
 816                 #if HAVE_DD
 817                 var uintDD carry0 = 0;
 818                 var uintD carry1 = 0;
 819                 var uintL shift = 0; // shift next digit before adding it to the carry, >=0, <intDsize
 820                 for (i = 0; i < N; i++) {
 821                         // Fetch z[i].re and round it to the nearest uintDD.
 822                         var uintDD digit = (uintDD)(z[i].re + (fftc_real)0.5);
 823                         if (shift > 0) {
 824                                 carry1 += (uintD)(digit >> (2*intDsize-shift));
 825                                 digit = digit << shift;
 826                         }
 827                         if ((carry0 += digit) < digit)
 828                                 carry1 += 1;
 829                         shift += l;
 830                         if (shift >= intDsize) {
 831                                 lsprefnext(destptr) = lowD(carry0);
 832                                 carry0 = highlowDD(carry1,highD(carry0));
 833                                 carry1 = 0;
 834                                 shift -= intDsize;
 835                         }
 836                 }
 837                 lsprefnext(destptr) = lowD(carry0);
 838                 lsprefnext(destptr) = highD(carry0);
 839                 lsprefnext(destptr) = carry1;
 840                 #else
 841                 var uintD carry0 = 0;
 842                 var uintD carry1 = 0;
 843                 var uintD carry2 = 0;
 844                 var uintL shift = 0; // shift next digit before adding it to the carry, >=0, <intDsize
 845                 for (i = 0; i < N; i++) {
 846                         // Fetch z[i].re and round it to the nearest uintDD.
 847                         var fftc_real zi = z[i].re + (fftc_real)0.5;
 848                         var const fftc_real pow2_intDsize = fftc_pow2(intDsize);
 849                         var const fftc_real invpow2_intDsize = (fftc_real)1.0/pow2_intDsize;
 850                         var uintD digit1 = (uintD)(zi * invpow2_intDsize);
 851                         var uintD digit0 = (uintD)(zi - (fftc_real)digit1 * pow2_intDsize);
 852                         if (shift > 0) {
 853                                 carry2 += digit1 >> (intDsize-shift);
 854                                 digit1 = (digit1 << shift) | (digit0 >> (intDsize-shift));
 855                                 digit0 = digit0 << shift;
 856                         }
 857                         if ((carry0 += digit0) < digit0)
 858                                 if ((carry1 += 1) == 0)
 859                                         carry2 += 1;
 860                         if ((carry1 += digit1) < digit1)
 861                                 carry2 += 1;
 862                         shift += l;
 863                         if (shift >= intDsize) {
 864                                 lsprefnext(destptr) = carry0;
 865                                 carry0 = carry1;
 866                                 carry1 = carry2;
 867                                 carry2 = 0;
 868                                 shift -= intDsize;
 869                         }
 870                 }
 871                 lsprefnext(destptr) = carry0;
 872                 lsprefnext(destptr) = carry1;
 873                 lsprefnext(destptr) = carry2;
 874                 #endif
 875         } else {
 876                 // 1-digit+1-float carry is sufficient
 877                 var uintD carry0 = 0;
 878                 var fftc_real carry1 = 0;
 879                 var uintL shift = 0; // shift next digit before adding it to the carry, >=0, <intDsize
 880                 for (i = 0; i < N; i++) {
 881                         // Fetch z[i].re and round it to the nearest integer.
 882                         var fftc_real digit = fftc_fround(z[i].re);
 883                         #ifdef DEBUG_FFTC
 884                         if (!(digit >= (fftc_real)0
 885                               && z[i].re > digit - (fftc_real)0.5
 886                               && z[i].re < digit + (fftc_real)0.5))
 887                                 throw runtime_exception();
 888                         #endif
 889                         if (shift > 0)
 890                                 digit = digit * fftc_pow2_table[shift];
 891                         var fftc_real digit1 = fftc_ffloor(digit*((fftc_real)1.0/fftc_pow2(intDsize)));
 892                         var uintD digit0 = (uintD)(digit - digit1*fftc_pow2(intDsize));
 893                         carry1 += digit1;
 894                         if ((carry0 += digit0) < digit0)
 895                                 carry1 += (fftc_real)1.0;
 896                         shift += l;
 897                         while (shift >= intDsize) {
 898                                 lsprefnext(destptr) = carry0;
 899                                 var fftc_real tmp = fftc_ffloor(carry1*((fftc_real)1.0/fftc_pow2(intDsize)));
 900                                 carry0 = (uintD)(carry1 - tmp*fftc_pow2(intDsize));
 901                                 carry1 = tmp;
 902                                 shift -= intDsize;
 903                         }
 904                 }
 905                 if (carry0 > 0 || carry1 > (fftc_real)0.0) {
 906                         lsprefnext(destptr) = carry0;
 907                         while (carry1 > (fftc_real)0.0) {
 908                                 var fftc_real tmp = fftc_ffloor(carry1*((fftc_real)1.0/fftc_pow2(intDsize)));
 909                                 lsprefnext(destptr) = (uintD)(carry1 - tmp*fftc_pow2(intDsize));
 910                                 carry1 = tmp;
 911                         }
 912                 }
 913         }
 914         return destptr;
 915 }
 916
 917 static inline void mulu_fftcomplex_nocheck (const uintD* sourceptr1, uintC len1,
 918                                             const uintD* sourceptr2, uintC len2,
 919                                             uintD* destptr)
 920 // Es ist 2 <= len1 <= len2.
 921 {
 922         // We have to find parameters l and k such that
 923         // ceiling(len1*intDsize/l) + ceiling(len2*intDsize/l) - 1 <= 2^k,
 924         // and (12*k-15)*e*2^(2*l+2*k+1/2-m) < 1/2.
 925         // Try primarily to minimize k. Minimizing l buys you nothing.
 926         var uintL k;
 927         // Computing k: If len1 and len2 differ much, we'll split source2 -
 928         // hence for the moment just substitute len1 for len2.
 929         //
 930         // First approximation of k: A necessary condition for
 931         // 2*ceiling(len1*intDsize/l) - 1 <= 2^k
 932         // is  2*len1*intDsize/l_max - 1 <= 2^k.
 933         {
 934                 var const int l = max_l(2);
 935                 var uintC lhs = 2*ceiling(len1*intDsize,l) - 1; // >=1
 936                 if (lhs < 3)
 937                         k = 2;
 938                 else
 939                         integerlengthC(lhs-1, k=); // k>=2
 940         }
 941         // Try whether this k is ok or whether we have to increase k.
 942         for ( ; ; k++) {
 943                 if (k >= sizeof(max_l_table)/sizeof(max_l_table[0])
 944                     || max_l_table[k] <= 0) {
 945                         throw runtime_exception("FFT problem: numbers too big, floating point precision not sufficient");
 946                 }
 947                 if (2*ceiling(len1*intDsize,max_l_table[k])-1 <= ((uintC)1 << k))
 948                         break;
 949         }
 950         // We could try to reduce l, keeping the same k. But why should we?
 951         // Calculate the number of pieces in which source2 will have to be
 952         // split. Each of the pieces must satisfy
 953         // ceiling(len1*intDsize/l) + ceiling(len2*intDsize/l) - 1 <= 2^k,
 954         var uintC len2p;
 955         // Try once with k, once with k+1. Compare them.
 956         {
 957                 var uintC remaining_k = ((uintC)1 << k) + 1 - ceiling(len1*intDsize,max_l_table[k]);
 958                 var uintC max_piecelen_k = floor(remaining_k*max_l_table[k],intDsize);
 959                 var uintC numpieces_k = ceiling(len2,max_piecelen_k);
 960                 var uintC remaining_k1 = ((uintC)1 << (k+1)) + 1 - ceiling(len1*intDsize,max_l_table[k+1]);
 961                 var uintC max_piecelen_k1 = floor(remaining_k1*max_l_table[k+1],intDsize);
 962                 var uintC numpieces_k1 = ceiling(len2,max_piecelen_k1);
 963                 if (numpieces_k <= 2*numpieces_k1) {
 964                         // keep k
 965                         len2p = max_piecelen_k;
 966                 } else {
 967                         // choose k+1
 968                         k = k+1;
 969                         len2p = max_piecelen_k1;
 970                 }
 971         }
 972         var const uintL l = max_l_table[k];
 973         var const uintL n = k;
 974         var const uintC N = (uintC)1 << n;
 975         CL_ALLOCA_STACK;
 976         var fftc_complex* const x = cl_alloc_array(fftc_complex,N);
 977         var fftc_complex* const y = cl_alloc_array(fftc_complex,N);
 978         #ifdef DEBUG_FFTC
 979         var fftc_complex* const z = cl_alloc_array(fftc_complex,N);
 980         #else
 981         var fftc_complex* const z = x; // put z in place of x - saves memory
 982         #endif
 983         var uintD* const tmpprod1 = cl_alloc_array(uintD,len1+1);
 984         var uintC tmpprod_len = floor(l<<n,intDsize)+6;
 985         var uintD* const tmpprod = cl_alloc_array(uintD,tmpprod_len);
 986         var uintC destlen = len1+len2;
 987         clear_loop_lsp(destptr,destlen);
 988         do {
 989                 if (len2p > len2)
 990                         len2p = len2;
 991                 if (len2p == 1) {
 992                         // cheap case
 993                         var uintD* tmpptr = arrayLSDptr(tmpprod1,len1+1);
 994                         mulu_loop_lsp(lspref(sourceptr2,0),sourceptr1,tmpptr,len1);
 995                         if (addto_loop_lsp(tmpptr,destptr,len1+1))
 996                                 if (inc_loop_lsp(destptr lspop (len1+1),destlen-(len1+1)))
 997                                         throw runtime_exception();
 998                 } else {
 999                         var bool squaring = ((sourceptr1 == sourceptr2) && (len1 == len2p));
1000                         // Fill factor x.
1001                         fill_factor(N,x,l,sourceptr1,len1);
1002                         // Fill factor y.
1003                         if (!squaring)
1004                                 fill_factor(N,y,l,sourceptr2,len2p);
1005                         // Multiply.
1006                         if (!squaring)
1007                                 fftc_convolution(n,N, &x[0], &y[0], &z[0]);
1008                         else
1009                                 fftc_convolution(n,N, &x[0], &x[0], &z[0]);
1010                         #ifdef DEBUG_FFTC
1011                         // Check result.
1012                         {
1013                                 var fftc_real re_lo_limit = (fftc_real)(-0.5);
1014                                 var fftc_real re_hi_limit = (fftc_real)N * fftc_pow2_table[l] * fftc_pow2_table[l] + (fftc_real)0.5;
1015                                 var fftc_real im_lo_limit = (fftc_real)(-0.5);
1016                                 var fftc_real im_hi_limit = (fftc_real)0.5;
1017                                 for (var uintC i = 0; i < N; i++) {
1018                                         if (!(z[i].im > im_lo_limit
1019                                               && z[i].im < im_hi_limit))
1020                                                 throw runtime_exception();
1021                                         if (!(z[i].re > re_lo_limit
1022                                               && z[i].re < re_hi_limit))
1023                                                 throw runtime_exception();
1024                                 }
1025                         }
1026                         #endif
1027                         var uintD* tmpLSDptr = arrayLSDptr(tmpprod,tmpprod_len);
1028                         var uintD* tmpMSDptr = unfill_product(n,N,z,l,tmpLSDptr);
1029                         var uintC tmplen =
1030                           #if CL_DS_BIG_ENDIAN_P
1031                             tmpLSDptr - tmpMSDptr;
1032                           #else
1033                             tmpMSDptr - tmpLSDptr;
1034                           #endif
1035                         if (tmplen > tmpprod_len)
1036                           throw runtime_exception();
1037                         // Add result to destptr[-destlen..-1]:
1038                         if (tmplen > destlen) {
1039                                 if (test_loop_msp(tmpMSDptr,tmplen-destlen))
1040                                         throw runtime_exception();
1041                                 tmplen = destlen;
1042                         }
1043                         if (addto_loop_lsp(tmpLSDptr,destptr,tmplen))
1044                                 if (inc_loop_lsp(destptr lspop tmplen,destlen-tmplen))
1045                                         throw runtime_exception();
1046                 }
1047                 // Decrement len2.
1048                 destptr = destptr lspop len2p;
1049                 destlen -= len2p;
1050                 sourceptr2 = sourceptr2 lspop len2p;
1051                 len2 -= len2p;
1052         } while (len2 > 0);
1053 }
1054
1055 #ifndef _CHECKSUM
1056 #define _CHECKSUM
1057
1058 // Compute a checksum: number mod (2^intDsize-1).
1059 static uintD compute_checksum (const uintD* sourceptr, uintC len)
1060 {
1061         var uintD tmp = ~(uintD)0; // -1-(sum mod 2^intDsize-1), always >0
1062         do {
1063                 var uintD digit = lsprefnext(sourceptr);
1064                 if (digit < tmp)
1065                         tmp -= digit; // subtract digit
1066                 else
1067                         tmp -= digit+1; // subtract digit-(2^intDsize-1)
1068         } while (--len > 0);
1069         return ~tmp;
1070 }
1071
1072 // Multiply two checksums modulo (2^intDsize-1).
1073 static inline uintD multiply_checksum (uintD checksum1, uintD checksum2)
1074 {
1075         var uintD checksum;
1076         var uintD cksum_hi;
1077         #if HAVE_DD
1078         var uintDD cksum = muluD(checksum1,checksum2);
1079         cksum_hi = highD(cksum); checksum = lowD(cksum);
1080         #else
1081         muluD(checksum1,checksum2, cksum_hi =, checksum =);
1082         #endif
1083         if ((checksum += cksum_hi) + 1 <= cksum_hi)
1084                 checksum += 1;
1085         return checksum;
1086 }
1087
1088 #endif // _CHECKSUM
1089
1090 static void mulu_fftcomplex (const uintD* sourceptr1, uintC len1,
1091                              const uintD* sourceptr2, uintC len2,
1092                              uintD* destptr)
1093 {
1094         // Compute checksums of the arguments and multiply them.
1095         var uintD checksum1 = compute_checksum(sourceptr1,len1);
1096         var uintD checksum2 = compute_checksum(sourceptr2,len2);
1097         var uintD checksum = multiply_checksum(checksum1,checksum2);
1098         mulu_fftcomplex_nocheck(sourceptr1,len1,sourceptr2,len2,destptr);
1099         if (!(checksum == compute_checksum(destptr,len1+len2))) {
1100                 throw runtime_exception("FFT problem: checksum error");
1101         }
1102 }