src/base/digitseq/cl_DS_mul_fftr.h

   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
   6
   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.
  10
  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).
  19 //
  20 // Based on this insight, we replace the usual n FFT steps
  21 //    for m = n...0:
  22 //      (z^(2^n)-1) = prod(j=0..2^(n-m)-1, (z^(2^m) - exp(2 pi j/2^(n-m))))
  23 // by
  24 //    for m = n...1:
  25 //      (z^(2^n)-1) = (z^(2^m)-1) * prod(j=1..2^(n-m)-1, factor_m[j](z))
  26 //      where
  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.
  39 //
  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])
  50 //    Invariant:
  51 //      for m = n..0:
  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
  59 //      for s in {0}:
  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))
  78 //    Invariant:
  79 //      for m = n..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].
  87 //
  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:
  91 //   n = 0: void.
  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:
  98 //   n = 0: void.
  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
 103 // for n>=4.
 104
 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.
 109 //
 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.
 113 //
 114 // Here is a more careful analysis, using absolute error estimates.
 115 //
 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
 125 //          <= e*2^(a-b-m).
 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
 161 //    5.iv. we'll have
 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),
 171 //    and
 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
 185 //          E_0 = 0,
 186 //          E_1 = 0,
 187 //          E_(n+1) = ...
 188 //
 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"!!!
 192
 193
 194 #if !(intDsize==32)
 195 #error "real fft implemented only for intDsize==32"
 196 #endif
 197
 198
 199 #include "cln/floatparam.h"
 200 #include "cln/io.h"
 201 #include "cln/abort.h"
 202
 203
 204 #if defined(HAVE_LONGDOUBLE) && (long_double_mant_bits > double_mant_bits) && (defined(__i386__) || defined(__m68k__) || (defined(__sparc__) && 0))
 205 // Only these CPUs have fast "long double"s in hardware.
 206 // On SPARC, "long double"s are emulated in software and don't work.
 207 typedef long double fftr_real;
 208 #define fftr_real_mant_bits long_double_mant_bits
 209 #define fftr_real_rounds long_double_rounds
 210 #else
 211 typedef double fftr_real;
 212 #define fftr_real_mant_bits double_mant_bits
 213 #define fftr_real_rounds double_rounds
 214 #endif
 215
 216 // We need complex numbers for precomputing the roots of unity.
 217 typedef struct fftr_complex {
 218         fftr_real re;
 219         fftr_real im;
 220 } fftr_complex;
 221
 222 // For every integer exp, we need to precompute  gam = 2 cos(2 pi exp/N).
 223 typedef union {
 224         fftr_complex c; // exp(2 pi i exp/N)
 225         struct {
 226                 fftr_real d; // gam = 2 cos(2 pi exp/N)
 227                 fftr_real e; // gam^2-1
 228                 fftr_real f; // 1/gam
 229         } gam;
 230 } fftr_cosinus;
 231
 232 static const fftr_complex fftr_roots_of_1 [32+1] =
 233   // roots_of_1[n] is a (2^n)th root of unity in C.
 234   // Also roots_of_1[n-1] = roots_of_1[n]^2.
 235   // For simplicity we choose  roots_of_1[n] = exp(2 pi i/2^n).
 236   {
 237    #if (fftr_real_mant_bits == double_mant_bits)
 238     // These values have 64 bit precision.
 239     { 1.0,                    0.0 },
 240     { -1.0,                   0.0 },
 241     { 0.0,                    1.0 },
 242     { 0.7071067811865475244,  0.7071067811865475244 },
 243     { 0.9238795325112867561,  0.38268343236508977172 },
 244     { 0.9807852804032304491,  0.19509032201612826784 },
 245     { 0.99518472667219688623, 0.098017140329560601996 },
 246     { 0.9987954562051723927,  0.049067674327418014254 },
 247     { 0.9996988186962042201,  0.024541228522912288032 },
 248     { 0.99992470183914454094, 0.0122715382857199260795 },
 249     { 0.99998117528260114264, 0.0061358846491544753597 },
 250     { 0.9999952938095761715,  0.00306795676296597627 },
 251     { 0.99999882345170190993, 0.0015339801862847656123 },
 252     { 0.99999970586288221914, 7.6699031874270452695e-4 },
 253     { 0.99999992646571785114, 3.8349518757139558907e-4 },
 254     { 0.9999999816164292938,  1.9174759731070330744e-4 },
 255     { 0.9999999954041073129,  9.5873799095977345874e-5 },
 256     { 0.99999999885102682754, 4.793689960306688455e-5 },
 257     { 0.99999999971275670683, 2.3968449808418218729e-5 },
 258     { 0.9999999999281891767,  1.1984224905069706422e-5 },
 259     { 0.99999999998204729416, 5.9921124526424278428e-6 },
 260     { 0.99999999999551182357, 2.9960562263346607504e-6 },
 261     { 0.99999999999887795586, 1.4980281131690112288e-6 },
 262     { 0.999999999999719489,   7.4901405658471572114e-7 },
 263     { 0.99999999999992987223, 3.7450702829238412391e-7 },
 264     { 0.99999999999998246807, 1.8725351414619534487e-7 },
 265     { 0.999999999999995617,   9.36267570730980828e-8 },
 266     { 0.99999999999999890425, 4.6813378536549092695e-8 },
 267     { 0.9999999999999997261,  2.340668926827455276e-8 },
 268     { 0.99999999999999993153, 1.1703344634137277181e-8 },
 269     { 0.99999999999999998287, 5.8516723170686386908e-9 },
 270     { 0.9999999999999999957,  2.925836158534319358e-9 },
 271     { 0.9999999999999999989,  1.4629180792671596806e-9 }
 272    #else (fftr_real_mant_bits > double_mant_bits)
 273     // These values have 128 bit precision.
 274     { 1.0L,                                       0.0L },
 275     { -1.0L,                                      0.0L },
 276     { 0.0L,                                       1.0L },
 277     { 0.707106781186547524400844362104849039284L, 0.707106781186547524400844362104849039284L },
 278     { 0.923879532511286756128183189396788286823L, 0.38268343236508977172845998403039886676L },
 279     { 0.980785280403230449126182236134239036975L, 0.195090322016128267848284868477022240928L },
 280     { 0.995184726672196886244836953109479921574L, 0.098017140329560601994195563888641845861L },
 281     { 0.998795456205172392714771604759100694444L, 0.0490676743274180142549549769426826583147L },
 282     { 0.99969881869620422011576564966617219685L,  0.0245412285229122880317345294592829250654L },
 283     { 0.99992470183914454092164649119638322435L,  0.01227153828571992607940826195100321214037L },
 284     { 0.999981175282601142656990437728567716173L, 0.00613588464915447535964023459037258091705L },
 285     { 0.999995293809576171511580125700119899554L, 0.00306795676296597627014536549091984251894L },
 286     { 0.99999882345170190992902571017152601905L,  0.001533980186284765612303697150264079079954L },
 287     { 0.999999705862882219160228217738765677117L, 7.66990318742704526938568357948576643142e-4L },
 288     { 0.99999992646571785114473148070738785695L,  3.83495187571395589072461681181381263396e-4L },
 289     { 0.999999981616429293808346915402909714504L, 1.91747597310703307439909561989000933469e-4L },
 290     { 0.99999999540410731289097193313960614896L,  9.58737990959773458705172109764763511872e-5L },
 291     { 0.9999999988510268275626733077945541084L,   4.79368996030668845490039904946588727468e-5L },
 292     { 0.99999999971275670684941397221864177609L,  2.39684498084182187291865771650218200947e-5L },
 293     { 0.999999999928189176709775095883850490262L, 1.198422490506970642152156159698898480473e-5L },
 294     { 0.99999999998204729417728262414778410738L,  5.99211245264242784287971180889086172999e-6L },
 295     { 0.99999999999551182354431058417299732444L,  2.99605622633466075045481280835705981183e-6L },
 296     { 0.999999999998877955886077016551752536504L, 1.49802811316901122885427884615536112069e-6L },
 297     { 0.999999999999719488971519214794719584451L, 7.49014056584715721130498566730655637157e-7L },
 298     { 0.99999999999992987224287980123972873676L,  3.74507028292384123903169179084633177398e-7L },
 299     { 0.99999999999998246806071995015624773673L,  1.8725351414619534486882457659356361712e-7L },
 300     { 0.999999999999995617015179987529456656217L, 9.3626757073098082799067286680885620193e-8L },
 301     { 0.999999999999998904253794996881763834182L, 4.68133785365490926951155181385400969594e-8L },
 302     { 0.99999999999999972606344874922040343793L,  2.34066892682745527595054934190348440379e-8L },
 303     { 0.999999999999999931515862187305098514444L, 1.170334463413727718124621350323810379807e-8L },
 304     { 0.999999999999999982878965546826274482047L, 5.8516723170686386908097901008341396944e-9L },
 305     { 0.999999999999999995719741386706568611352L, 2.92583615853431935792823046906895590202e-9L },
 306     { 0.999999999999999998929935346676642152265L, 1.46291807926715968052953216186596371037e-9L }
 307    #endif
 308   };
 309
 310 // Define this for (cheap) consistency checks.
 311 #define DEBUG_FFTR
 312
 313 static fftr_real fftr_pow2_table[64] = // table of powers of 2
 314   {
 315                       1.0,
 316                       2.0,
 317                       4.0,
 318                       8.0,
 319                      16.0,
 320                      32.0,
 321                      64.0,
 322                     128.0,
 323                     256.0,
 324                     512.0,
 325                    1024.0,
 326                    2048.0,
 327                    4096.0,
 328                    8192.0,
 329                   16384.0,
 330                   32768.0,
 331                   65536.0,
 332                  131072.0,
 333                  262144.0,
 334                  524288.0,
 335                 1048576.0,
 336                 2097152.0,
 337                 4194304.0,
 338                 8388608.0,
 339                16777216.0,
 340                33554432.0,
 341                67108864.0,
 342               134217728.0,
 343               268435456.0,
 344               536870912.0,
 345              1073741824.0,
 346              2147483648.0,
 347              4294967296.0,
 348              8589934592.0,
 349             17179869184.0,
 350             34359738368.0,
 351             68719476736.0,
 352            137438953472.0,
 353            274877906944.0,
 354            549755813888.0,
 355           1099511627776.0,
 356           2199023255552.0,
 357           4398046511104.0,
 358           8796093022208.0,
 359          17592186044416.0,
 360          35184372088832.0,
 361          70368744177664.0,
 362         140737488355328.0,
 363         281474976710656.0,
 364         562949953421312.0,
 365        1125899906842624.0,
 366        2251799813685248.0,
 367        4503599627370496.0,
 368        9007199254740992.0,
 369       18014398509481984.0,
 370       36028797018963968.0,
 371       72057594037927936.0,
 372      144115188075855872.0,
 373      288230376151711744.0,
 374      576460752303423488.0,
 375     1152921504606846976.0,
 376     2305843009213693952.0,
 377     4611686018427387904.0,
 378     9223372036854775808.0
 379   };
 380
 381 // For a constant expression n (0 <= n < 128), returns 2^n of type fftr_real.
 382 #define fftr_pow2(n)  \
 383   (((n) & 64 ? (fftr_real)18446744073709551616.0 : (fftr_real)1.0)      \
 384    * ((n) & 32 ? (fftr_real)4294967296.0 : (fftr_real)1.0)              \
 385    * ((n) & 16 ? (fftr_real)65536.0 : (fftr_real)1.0)                   \
 386    * ((n) & 8 ? (fftr_real)256.0 : (fftr_real)1.0)                      \
 387    * ((n) & 4 ? (fftr_real)16.0 : (fftr_real)1.0)                       \
 388    * ((n) & 2 ? (fftr_real)4.0 : (fftr_real)1.0)                        \
 389    * ((n) & 1 ? (fftr_real)2.0 : (fftr_real)1.0)                        \
 390   )
 391
 392 // r := a * b
 393 static inline void mul (const fftr_complex& a, const fftr_complex& b, fftr_complex& r)
 394 {
 395         var fftr_real r_re = a.re * b.re - a.im * b.im;
 396         var fftr_real r_im = a.re * b.im + a.im * b.re;
 397         r.re = r_re;
 398         r.im = r_im;
 399 }
 400
 401 static uintL shuffle (uintL n, uintL x)
 402 {
 403         var uintL y = 0;
 404         var sintL v = 1;
 405         // Invariant: y + v*shuffle(n,x).
 406         do {
 407                 if (x & 1)
 408                         if (x == 1)
 409                                 return y + (v << (n-1));
 410                         else {
 411                                 y = y + (v << n);
 412                                 v = -v;
 413                         }
 414                 x >>= 1;
 415         } while (!(--n == 0));
 416         cl_abort();
 417 }
 418
 419 #if 0 // unused
 420 static uintL invshuffle (uintL n, uintL x)
 421 {
 422         var uintL y = 0;
 423         var uintL v = 1;
 424         // Invariant: y + v*invshuffle(n,x).
 425         do {
 426                 if (x == ((uintL)1 << (n-1)))
 427                         return y + v;
 428                 else if (x > ((uintL)1 << (n-1))) {
 429                         x = ((uintL)1 << n) - x;
 430                         y = y+v;
 431                 }
 432                 v <<= 1;
 433         } while (!(--n == 0));
 434         cl_abort();
 435 }
 436 #endif
 437
 438 // Compute a real convolution using FFT: z[0..N-1] := x[0..N-1] * y[0..N-1].
 439 static void fftr_convolution (const uintL n, const uintL N, // N = 2^n
 440                               fftr_real * x, // N numbers
 441                               fftr_real * y, // N numbers
 442                               fftr_real * z  // N numbers result
 443                              )
 444 {
 445         CL_ALLOCA_STACK;
 446         var fftr_cosinus* const w = cl_alloc_array(fftr_cosinus,N>>2);
 447         var uintL i;
 448         // Initialize w[exp].c to exp(2 pi i exp/N).
 449         w[0].c = fftr_roots_of_1[0];
 450         {
 451                 var int j;
 452                 for (j = n-3; j>=0; j--) {
 453                         var fftr_complex r_j = fftr_roots_of_1[n-j];
 454                         w[1<<j].c = r_j;
 455                         for (i = (2<<j); i < N>>2; i += (2<<j))
 456                                 mul(w[i].c,r_j, w[i+(1<<j)].c);
 457                 }
 458         }
 459         // Initialize w[exp] to 2 cos(2 pi exp/N).
 460         for (i = 0; i < N>>2; i++) {
 461                 var fftr_real gam = w[i].c.re * (fftr_real)2.0;
 462                 w[i].gam.d = gam;
 463                 w[i].gam.e = (gam - (fftr_real)1.0) * (gam + (fftr_real)1.0);
 464                 w[i].gam.f = (fftr_real)1.0 / gam;
 465         }
 466         var bool squaring = (x == y);
 467         // Do an FFT of length N on x.
 468         {
 469                 var uintL l;
 470                 for (l = n-1; l > 0; l--) {
 471                         /* s = 0 */ {
 472                                 var const uintL tmax = (uintL)1 << l;
 473                                 for (var uintL t = 0; t < tmax; t++) {
 474                                         var uintL i1 = t;
 475                                         var uintL i2 = i1 + tmax;
 476                                         // replace (x[i1],x[i2]) by
 477                                         // (x[i1] + x[i2], x[i1] - x[i2])
 478                                         var fftr_real tmp;
 479                                         tmp = x[i2];
 480                                         x[i2] = x[i1] - tmp;
 481                                         x[i1] = x[i1] + tmp;
 482                                 }
 483                         }
 484                         var const uintL smax = (uintL)1 << (n-1-l);
 485                         var const uintL tmax = (uintL)1 << (l-1);
 486                         for (var uintL s = 1; s < smax; s++) {
 487                                 var uintL exp = shuffle(n-1-l,s) << (l-1);
 488                                 for (var uintL t = 0; t < tmax; t++) {
 489                                         var uintL i1 = (s << (l+1)) + t;
 490                                         var uintL i2 = i1 + tmax;
 491                                         var uintL i3 = i2 + tmax;
 492                                         var uintL i4 = i3 + tmax;
 493                                         // replace (x[i1],x[i2],x[i3],x[i4]) by
 494                                         // (x[i1]-x[i3] - x[i4]*gam,
 495                                         //  x[i3]*gam + x[i2]+x[i4]*(gam^2-1),
 496                                         //  x[i1]-x[i3] + x[i4]*gam,
 497                                         //  - x[i3]*gam + x[i2]+x[i4]*(gam^2-1))
 498                                         var fftr_real sum13;
 499                                         var fftr_real sum24;
 500                                         var fftr_real tmp3;
 501                                         var fftr_real tmp4;
 502                                         sum13 = x[i1] - x[i3];
 503                                         tmp3 = x[i3]*w[exp].gam.d;
 504                                         tmp4 = x[i4]*w[exp].gam.d;
 505                                         sum24 = x[i2]+x[i4]*w[exp].gam.e;
 506                                         x[i1] = sum13 - tmp4;
 507                                         x[i3] = sum13 + tmp4;
 508                                         x[i2] = sum24 + tmp3;
 509                                         x[i4] = sum24 - tmp3;
 510                                 }
 511                         }
 512                 }
 513                 /* l = 0 */ {
 514                         // replace (x[0],x[1]) by (x[0]+x[1], x[0]-x[1])
 515                         var fftr_real tmp;
 516                         tmp = x[1];
 517                         x[1] = x[0] - tmp;
 518                         x[0] = x[0] + tmp;
 519                 }
 520         }
 521         // Do an FFT of length N on y.
 522         if (!squaring) {
 523                 var uintL l;
 524                 for (l = n-1; l > 0; l--) {
 525                         /* s = 0 */ {
 526                                 var const uintL tmax = (uintL)1 << l;
 527                                 for (var uintL t = 0; t < tmax; t++) {
 528                                         var uintL i1 = t;
 529                                         var uintL i2 = i1 + tmax;
 530                                         // replace (y[i1],y[i2]) by
 531                                         // (y[i1] + y[i2], y[i1] - y[i2])
 532                                         var fftr_real tmp;
 533                                         tmp = y[i2];
 534                                         y[i2] = y[i1] - tmp;
 535                                         y[i1] = y[i1] + tmp;
 536                                 }
 537                         }
 538                         var const uintL smax = (uintL)1 << (n-1-l);
 539                         var const uintL tmax = (uintL)1 << (l-1);
 540                         for (var uintL s = 1; s < smax; s++) {
 541                                 var uintL exp = shuffle(n-1-l,s) << (l-1);
 542                                 for (var uintL t = 0; t < tmax; t++) {
 543                                         var uintL i1 = (s << (l+1)) + t;
 544                                         var uintL i2 = i1 + tmax;
 545                                         var uintL i3 = i2 + tmax;
 546                                         var uintL i4 = i3 + tmax;
 547                                         // replace (y[i1],y[i2],y[i3],y[i4]) by
 548                                         // (y[i1]-y[i3] - y[i4]*gam,
 549                                         //  y[i3]*gam + y[i2]+y[i4]*(gam^2-1),
 550                                         //  y[i1]-y[i3] + y[i4]*gam,
 551                                         //  - y[i3]*gam + y[i2]+y[i4]*(gam^2-1))
 552                                         var fftr_real sum13;
 553                                         var fftr_real sum24;
 554                                         var fftr_real tmp3;
 555                                         var fftr_real tmp4;
 556                                         sum13 = y[i1] - y[i3];
 557                                         tmp3 = y[i3]*w[exp].gam.d;
 558                                         tmp4 = y[i4]*w[exp].gam.d;
 559                                         sum24 = y[i2]+y[i4]*w[exp].gam.e;
 560                                         y[i1] = sum13 - tmp4;
 561                                         y[i3] = sum13 + tmp4;
 562                                         y[i2] = sum24 + tmp3;
 563                                         y[i4] = sum24 - tmp3;
 564                                 }
 565                         }
 566                 }
 567                 /* l = 0 */ {
 568                         // replace (y[0],y[1]) by (y[0]+y[1], y[0]-y[1])
 569                         var fftr_real tmp;
 570                         tmp = y[1];
 571                         y[1] = y[0] - tmp;
 572                         y[0] = y[0] + tmp;
 573                 }
 574         }
 575         // Multiply the transformed vectors into z.
 576         {
 577                 // Multiplication mod (z-1).
 578                 z[0] = x[0] * y[0];
 579                 // Multiplication mod (z+1).
 580                 z[1] = x[1] * y[1];
 581                 #if 0 // This needs w[1..2^(n-1)-1].
 582                 var const uintL smax = (uintL)1 << (n-1);
 583                 for (var uintL s = 1; s < smax; s++) {
 584                         // Multiplication mod (z^2 - 2 cos(2 pi j/2^n) z + 1)
 585                         // with j = shuffle(n-1,s).
 586                         var uintL exp = shuffle(n-1,s);
 587                         var fftr_real tmp0 = x[2*s] * y[2*s];
 588                         var fftr_real tmp1 = x[2*s] * y[2*s+1] + x[2*s+1] * y[2*s];
 589                         var fftr_real tmp2 = x[2*s+1] * y[2*s+1];
 590                         z[2*s] = tmp0 - tmp2;
 591                         z[2*s+1] = tmp1 + tmp2*w[exp].gam.d;
 592                 }
 593                 #else // This only needs w[1..2^(n-2)-1].
 594                 // Multiplication mod (z^2+1).
 595                 /* s = 1 */ {
 596                         var fftr_real tmp0 = x[2] * y[2];
 597                         var fftr_real tmp1 = x[2] * y[3] + x[3] * y[2];
 598                         var fftr_real tmp2 = x[3] * y[3];
 599                         z[2] = tmp0 - tmp2;
 600                         z[3] = tmp1;
 601                 }
 602                 var const uintL smax = (uintL)1 << (n-2);
 603                 for (var uintL s = 1; s < smax; s++) {
 604                         var uintL exp = shuffle(n-2,s);
 605                         // Multiplication mod (z^2 - 2 cos(2 pi j/2^n) z + 1)
 606                         // with j = shuffle(n-1,2*s) = shuffle(n-2,s).
 607                         var fftr_real gam = w[exp].gam.d;
 608                         {
 609                                 var fftr_real tmp0 = x[4*s] * y[4*s];
 610                                 var fftr_real tmp1 = x[4*s] * y[4*s+1] + x[4*s+1] * y[4*s];
 611                                 var fftr_real tmp2 = x[4*s+1] * y[4*s+1];
 612                                 z[4*s] = tmp0 - tmp2;
 613                                 z[4*s+1] = tmp1 + tmp2 * gam;
 614                         }
 615                         // Multiplication mod (z^2 - 2 cos(2 pi j/2^n) z + 1)
 616                         // with j = shuffle(n-1,2*s+1) = 2^(n-1) - shuffle(n-2,s).
 617                         {
 618                                 var fftr_real tmp0 = x[4*s+2] * y[4*s+2];
 619                                 var fftr_real tmp1 = x[4*s+2] * y[4*s+3] + x[4*s+3] * y[4*s+2];
 620                                 var fftr_real tmp2 = x[4*s+3] * y[4*s+3];
 621                                 z[4*s+2] = tmp0 - tmp2;
 622                                 z[4*s+3] = tmp1 - tmp2 * gam;
 623                         }
 624                 }
 625                 #endif
 626         }
 627         // Undo an FFT of length N on z.
 628         {
 629                 var uintL l;
 630                 /* l = 0 */ {
 631                         // replace (z[0],z[1]) by ((z[0]+z[1])/2, (z[0]-z[1])/2)
 632                         var fftr_real tmp;
 633                         tmp = z[1];
 634                         z[1] = (z[0] - tmp) * (fftr_real)0.5;
 635                         z[0] = (z[0] + tmp) * (fftr_real)0.5;
 636                 }
 637                 for (l = 1; l < n; l++) {
 638                         /* s = 0 */ {
 639                                 var const uintL tmax = (uintL)1 << l;
 640                                 for (var uintL t = 0; t < tmax; t++) {
 641                                         var uintL i1 = t;
 642                                         var uintL i2 = i1 + tmax;
 643                                         // replace (z[i1],z[i2]) by
 644                                         // ((z[i1]+z[i2])/2, (z[i1]-z[i2])/2)
 645                                         // Do the division by 2 later.
 646                                         var fftr_real tmp;
 647                                         tmp = z[i2];
 648                                         z[i2] = z[i1] - tmp;
 649                                         z[i1] = z[i1] + tmp;
 650                                 }
 651                         }
 652                         var const uintL smax = (uintL)1 << (n-1-l);
 653                         var const uintL tmax = (uintL)1 << (l-1);
 654                         for (var uintL s = 1; s < smax; s++) {
 655                                 var uintL exp = shuffle(n-1-l,s) << (l-1);
 656                                 for (var uintL t = 0; t < tmax; t++) {
 657                                         var uintL i1 = (s << (l+1)) + t;
 658                                         var uintL i2 = i1 + tmax;
 659                                         var uintL i3 = i2 + tmax;
 660                                         var uintL i4 = i3 + tmax;
 661                                         // replace (z[i1],z[i2],z[i3],z[i4]) by
 662                                         // ((z[i1]+z[i3]+(z[i2]-z[i4])/gam)/2,
 663                                         //  (z[i2]+z[i4]-(z[i3]-z[i1])/gam*(gam^2-1))/2,
 664                                         //  (z[i2]-z[i4])/(gam*2),
 665                                         //  (z[i3]-z[i1])/(gam*2))
 666                                         // Do the division by 2 later.
 667                                         var fftr_real sum13;
 668                                         var fftr_real sum24;
 669                                         var fftr_real tmp3;
 670                                         var fftr_real tmp4;
 671                                         sum13 = z[i1] + z[i3];
 672                                         sum24 = z[i2] + z[i4];
 673                                         tmp4 = (z[i3] - z[i1]) * w[exp].gam.f;
 674                                         tmp3 = (z[i2] - z[i4]) * w[exp].gam.f;
 675                                         z[i1] = sum13 + tmp3;
 676                                         z[i2] = sum24 - tmp4*w[exp].gam.e;
 677                                         z[i3] = tmp3;
 678                                         z[i4] = tmp4;
 679                                 }
 680                         }
 681                 }
 682                 // Do all divisions by 2 now.
 683                 {
 684                         var fftr_real f = (fftr_real)2.0 / (fftr_real)N; // 2^-(n-1)
 685                         for (i = 0; i < N; i++)
 686                                 z[i] = z[i]*f;
 687                 }
 688         }
 689 }
 690
 691 // For a given k >= 2, the maximum l is determined by
 692 // 2*l < m - 1/2 - 2*k - log_2(12*k-15) - (1 if e=1.0, 0 if e=0.5).
 693 // This is a decreasing function of k.
 694 #define max_l(k)  \
 695         (int)((fftr_real_mant_bits                      \
 696                - 2*(k)                                  \
 697                - ((k)<=2 ? 4 : (k)<=3 ? 5 : (k)<=5 ? 6 : (k)<=8 ? 7 : (k)<=16 ? 8 : (k)<=31 ? 9 : 10) \
 698                - (fftr_real_rounds == rounds_to_nearest ? 0 : 1)) \
 699               / 2)
 700 static int max_l_table[32+1] =
 701   {         0,         0,  max_l(2),  max_l(3),  max_l(4),  max_l(5),  max_l(6),
 702      max_l(7),  max_l(8),  max_l(9), max_l(10), max_l(11), max_l(12), max_l(13),
 703     max_l(14), max_l(15), max_l(16), max_l(17), max_l(18), max_l(19), max_l(20),
 704     max_l(21), max_l(22), max_l(23), max_l(24), max_l(25), max_l(26), max_l(27),
 705     max_l(28), max_l(29), max_l(30), max_l(31), max_l(32)
 706   };
 707
 708 // Split len uintD's below sourceptr into chunks of l bits, thus filling
 709 // N real numbers at x.
 710 static void fill_factor (uintL N, fftr_real* x, uintL l,
 711                          const uintD* sourceptr, uintL len)
 712 {
 713         var uintL i;
 714         if (max_l(2) > intDsize && l > intDsize) {
 715                 // l > intDsize
 716                 if (max_l(2) > 64 && l > 64) {
 717                         fprint(std::cerr, "FFT problem: l > 64 not supported by pow2_table\n");
 718                         cl_abort();
 719                 }
 720                 var fftr_real carry = 0;
 721                 var sintL carrybits = 0; // number of bits in carry (>=0, <l)
 722                 i = 0;
 723                 while (len > 0) {
 724                         var uintD digit = lsprefnext(sourceptr);
 725                         if (carrybits+intDsize >= l) {
 726                                 x[i] = carry + (fftr_real)(digit & bitm(l-carrybits)) * fftr_pow2_table[carrybits];
 727                                 i++;
 728                                 carry = (l-carrybits == intDsize ? (fftr_real)0 : (fftr_real)(digit >> (l-carrybits)));
 729                                 carrybits = carrybits+intDsize-l;
 730                         } else {
 731                                 carry = carry + (fftr_real)digit * fftr_pow2_table[carrybits];
 732                                 carrybits = carrybits+intDsize;
 733                         }
 734                         len--;
 735                 }
 736                 if (carrybits > 0) {
 737                         x[i] = carry;
 738                         i++;
 739                 }
 740                 if (i > N)
 741                         cl_abort();
 742         } else if (max_l(2) >= intDsize && l == intDsize) {
 743                 // l = intDsize
 744                 if (len > N)
 745                         cl_abort();
 746                 for (i = 0; i < len; i++) {
 747                         var uintD digit = lsprefnext(sourceptr);
 748                         x[i] = (fftr_real)digit;
 749                 }
 750         } else {
 751                 // l < intDsize
 752                 var const uintD l_mask = bit(l)-1;
 753                 var uintD carry = 0;
 754                 var sintL carrybits = 0; // number of bits in carry (>=0, <intDsize)
 755                 for (i = 0; i < N; i++) {
 756                         if (carrybits >= l) {
 757                                 x[i] = (fftr_real)(carry & l_mask);
 758                                 carry >>= l;
 759                                 carrybits -= l;
 760                         } else {
 761                                 if (len == 0)
 762                                         break;
 763                                 len--;
 764                                 var uintD digit = lsprefnext(sourceptr);
 765                                 x[i] = (fftr_real)((carry | (digit << carrybits)) & l_mask);
 766                                 carry = digit >> (l-carrybits);
 767                                 carrybits = intDsize - (l-carrybits);
 768                         }
 769                 }
 770                 while (carrybits > 0) {
 771                         if (!(i < N))
 772                                 cl_abort();
 773                         x[i] = (fftr_real)(carry & l_mask);
 774                         carry >>= l;
 775                         carrybits -= l;
 776                         i++;
 777                 }
 778                 if (len > 0)
 779                         cl_abort();
 780         }
 781         for ( ; i < N; i++)
 782                 x[i] = (fftr_real)0;
 783 }
 784
 785 // Given a not too large floating point number, round it to the nearest integer.
 786 static inline fftr_real fftr_fround (fftr_real x)
 787 {
 788         return
 789           #if (fftr_real_rounds == rounds_to_nearest)
 790             (x + (fftr_pow2(fftr_real_mant_bits-1)+fftr_pow2(fftr_real_mant_bits-2)))
 791             - (fftr_pow2(fftr_real_mant_bits-1)+fftr_pow2(fftr_real_mant_bits-2));
 792           #elif (fftr_real_rounds == rounds_to_infinity)
 793             (fftr_pow2(fftr_real_mant_bits-1)+fftr_pow2(fftr_real_mant_bits-2))
 794             - ((fftr_pow2(fftr_real_mant_bits-1)+fftr_pow2(fftr_real_mant_bits-2))
 795                - (x + (fftr_real)0.5));
 796           #else // rounds_to_zero, rounds_to_minus_infinity
 797             ((x + (fftr_real)0.5)
 798              + (fftr_pow2(fftr_real_mant_bits-1)+fftr_pow2(fftr_real_mant_bits-2)))
 799             - (fftr_pow2(fftr_real_mant_bits-1)+fftr_pow2(fftr_real_mant_bits-2));
 800           #endif
 801 }
 802
 803 // Given a not too large floating point number, round it down.
 804 static inline fftr_real fftr_ffloor (fftr_real x)
 805 {
 806         #if (fftr_real_rounds == rounds_to_nearest)
 807           var fftr_real y =
 808             (x + (fftr_pow2(fftr_real_mant_bits-1)+fftr_pow2(fftr_real_mant_bits-2)))
 809             - (fftr_pow2(fftr_real_mant_bits-1)+fftr_pow2(fftr_real_mant_bits-2));
 810           if (y <= x)
 811                 return y;
 812           else
 813                 return y - (fftr_real)1.0;
 814         #elif (fftr_real_rounds == rounds_to_infinity)
 815           return
 816             (fftr_pow2(fftr_real_mant_bits-1)+fftr_pow2(fftr_real_mant_bits-2))
 817             - ((fftr_pow2(fftr_real_mant_bits-1)+fftr_pow2(fftr_real_mant_bits-2))
 818                - x);
 819         #else // rounds_to_zero, rounds_to_minus_infinity
 820           return
 821             (x + (fftr_pow2(fftr_real_mant_bits-1)+fftr_pow2(fftr_real_mant_bits-2)))
 822             - (fftr_pow2(fftr_real_mant_bits-1)+fftr_pow2(fftr_real_mant_bits-2));
 823         #endif
 824 }
 825
 826 // Combine the N real numbers at z into uintD's below destptr.
 827 // The z[i] are known to be approximately integers >= 0, < N*2^(2*l).
 828 // Assumes room for floor(N*l/intDsize)+(1+ceiling((n+2*l)/intDsize)) uintD's
 829 // below destptr. Fills len digits and returns (destptr lspop len).
 830 static uintD* unfill_product (uintL n, uintL N, // N = 2^n
 831                               const fftr_real * z, uintL l,
 832                               uintD* destptr)
 833 {
 834         var uintL i;
 835         if (n + 2*l <= intDsize) {
 836                 // 2-digit carry is sufficient, l < intDsize
 837                 var uintD carry0 = 0;
 838                 var uintD carry1 = 0;
 839                 var uintL shift = 0; // shift next digit before adding it to the carry, >=0, <intDsize
 840                 for (i = 0; i < N; i++) {
 841                         // Fetch z[i] and round it to the nearest uintD.
 842                         var uintD digit = (uintD)(z[i] + (fftr_real)0.5);
 843                         if (shift > 0) {
 844                                 carry1 += digit >> (intDsize-shift);
 845                                 digit = digit << shift;
 846                         }
 847                         if ((carry0 += digit) < digit)
 848                                 carry1 += 1;
 849                         shift += l;
 850                         if (shift >= intDsize) {
 851                                 lsprefnext(destptr) = carry0;
 852                                 carry0 = carry1;
 853                                 carry1 = 0;
 854                                 shift -= intDsize;
 855                         }
 856                 }
 857                 lsprefnext(destptr) = carry0;
 858                 lsprefnext(destptr) = carry1;
 859         } else if (n + 2*l <= 2*intDsize) {
 860                 // 3-digit carry is sufficient, l < intDsize
 861                 #if HAVE_DD
 862                 var uintDD carry0 = 0;
 863                 var uintD carry1 = 0;
 864                 var uintL shift = 0; // shift next digit before adding it to the carry, >=0, <intDsize
 865                 for (i = 0; i < N; i++) {
 866                         // Fetch z[i] and round it to the nearest uintDD.
 867                         var uintDD digit = (uintDD)(z[i] + (fftr_real)0.5);
 868                         if (shift > 0) {
 869                                 carry1 += (uintD)(digit >> (2*intDsize-shift));
 870                                 digit = digit << shift;
 871                         }
 872                         if ((carry0 += digit) < digit)
 873                                 carry1 += 1;
 874                         shift += l;
 875                         if (shift >= intDsize) {
 876                                 lsprefnext(destptr) = lowD(carry0);
 877                                 carry0 = highlowDD(carry1,highD(carry0));
 878                                 carry1 = 0;
 879                                 shift -= intDsize;
 880                         }
 881                 }
 882                 lsprefnext(destptr) = lowD(carry0);
 883                 lsprefnext(destptr) = highD(carry0);
 884                 lsprefnext(destptr) = carry1;
 885                 #else
 886                 var uintD carry0 = 0;
 887                 var uintD carry1 = 0;
 888                 var uintD carry2 = 0;
 889                 var uintL shift = 0; // shift next digit before adding it to the carry, >=0, <intDsize
 890                 for (i = 0; i < N; i++) {
 891                         // Fetch z[i] and round it to the nearest uintDD.
 892                         var fftr_real zi = z[i] + (fftr_real)0.5;
 893                         var const fftr_real pow2_intDsize = fftr_pow2(intDsize);
 894                         var const fftr_real invpow2_intDsize = (fftr_real)1.0/pow2_intDsize;
 895                         var uintD digit1 = (uintD)(zi * invpow2_intDsize);
 896                         var uintD digit0 = (uintD)(zi - (fftr_real)digit1 * pow2_intDsize);
 897                         if (shift > 0) {
 898                                 carry2 += digit1 >> (intDsize-shift);
 899                                 digit1 = (digit1 << shift) | (digit0 >> (intDsize-shift));
 900                                 digit0 = digit0 << shift;
 901                         }
 902                         if ((carry0 += digit0) < digit0)
 903                                 if ((carry1 += 1) == 0)
 904                                         carry2 += 1;
 905                         if ((carry1 += digit1) < digit1)
 906                                 carry2 += 1;
 907                         shift += l;
 908                         if (shift >= intDsize) {
 909                                 lsprefnext(destptr) = carry0;
 910                                 carry0 = carry1;
 911                                 carry1 = carry2;
 912                                 carry2 = 0;
 913                                 shift -= intDsize;
 914                         }
 915                 }
 916                 lsprefnext(destptr) = carry0;
 917                 lsprefnext(destptr) = carry1;
 918                 lsprefnext(destptr) = carry2;
 919                 #endif
 920         } else {
 921                 // 1-digit+1-float carry is sufficient
 922                 var uintD carry0 = 0;
 923                 var fftr_real carry1 = 0;
 924                 var uintL shift = 0; // shift next digit before adding it to the carry, >=0, <intDsize
 925                 for (i = 0; i < N; i++) {
 926                         // Fetch z[i] and round it to the nearest integer.
 927                         var fftr_real digit = fftr_fround(z[i]);
 928                         #ifdef DEBUG_FFTR
 929                         if (!(digit >= (fftr_real)0
 930                               && z[i] > digit - (fftr_real)0.5
 931                               && z[i] < digit + (fftr_real)0.5))
 932                                 cl_abort();
 933                         #endif
 934                         if (shift > 0)
 935                                 digit = digit * fftr_pow2_table[shift];
 936                         var fftr_real digit1 = fftr_ffloor(digit*((fftr_real)1.0/fftr_pow2(intDsize)));
 937                         var uintD digit0 = (uintD)(digit - digit1*fftr_pow2(intDsize));
 938                         carry1 += digit1;
 939                         if ((carry0 += digit0) < digit0)
 940                                 carry1 += (fftr_real)1.0;
 941                         shift += l;
 942                         while (shift >= intDsize) {
 943                                 lsprefnext(destptr) = carry0;
 944                                 var fftr_real tmp = fftr_ffloor(carry1*((fftr_real)1.0/fftr_pow2(intDsize)));
 945                                 carry0 = (uintD)(carry1 - tmp*fftr_pow2(intDsize));
 946                                 carry1 = tmp;
 947                                 shift -= intDsize;
 948                         }
 949                 }
 950                 if (carry0 > 0 || carry1 > (fftr_real)0.0) {
 951                         lsprefnext(destptr) = carry0;
 952                         while (carry1 > (fftr_real)0.0) {
 953                                 var fftr_real tmp = fftr_ffloor(carry1*((fftr_real)1.0/fftr_pow2(intDsize)));
 954                                 lsprefnext(destptr) = (uintD)(carry1 - tmp*fftr_pow2(intDsize));
 955                                 carry1 = tmp;
 956                         }
 957                 }
 958         }
 959         return destptr;
 960 }
 961
 962 static inline void mulu_fftr_nocheck (const uintD* sourceptr1, uintC len1,
 963                                       const uintD* sourceptr2, uintC len2,
 964                                       uintD* destptr)
 965 // Es ist 2 <= len1 <= len2.
 966 {
 967         // We have to find parameters l and k such that
 968         // ceiling(len1*intDsize/l) + ceiling(len2*intDsize/l) - 1 <= 2^k,
 969         // and (12*k-15)*e*2^(2*l+2*k+1/2-m) < 1/2.
 970         // Try primarily to minimize k. Minimizing l buys you nothing.
 971         var uintL k;
 972         // Computing k: If len1 and len2 differ much, we'll split source2 -
 973         // hence for the moment just substitute len1 for len2.
 974         //
 975         // First approximation of k: A necessary condition for
 976         // 2*ceiling(len1*intDsize/l) - 1 <= 2^k
 977         // is  2*len1*intDsize/l_max - 1 <= 2^k.
 978         {
 979                 var const int l = max_l(2);
 980                 var uintL lhs = 2*ceiling((uintL)len1*intDsize,l) - 1; // >=1
 981                 if (lhs < 3)
 982                         k = 2;
 983                 else
 984                         integerlength32(lhs-1, k=); // k>=2
 985         }
 986         // Try whether this k is ok or whether we have to increase k.
 987         for ( ; ; k++) {
 988                 if (k >= sizeof(max_l_table)/sizeof(max_l_table[0])
 989                     || max_l_table[k] <= 0) {
 990                         fprint(std::cerr, "FFT problem: numbers too big, floating point precision not sufficient\n");
 991                         cl_abort();
 992                 }
 993                 if (2*ceiling((uintL)len1*intDsize,max_l_table[k])-1 <= ((uintL)1 << k))
 994                         break;
 995         }
 996         // We could try to reduce l, keeping the same k. But why should we?
 997         // Calculate the number of pieces in which source2 will have to be
 998         // split. Each of the pieces must satisfy
 999         // ceiling(len1*intDsize/l) + ceiling(len2*intDsize/l) - 1 <= 2^k,
1000         var uintL len2p;
1001         // Try once with k, once with k+1. Compare them.
1002         {
1003                 var uintL remaining_k = ((uintL)1 << k) + 1 - ceiling((uintL)len1*intDsize,max_l_table[k]);
1004                 var uintL max_piecelen_k = floor(remaining_k*max_l_table[k],intDsize);
1005                 var uintL numpieces_k = ceiling(len2,max_piecelen_k);
1006                 var uintL remaining_k1 = ((uintL)1 << (k+1)) + 1 - ceiling((uintL)len1*intDsize,max_l_table[k+1]);
1007                 var uintL max_piecelen_k1 = floor(remaining_k1*max_l_table[k+1],intDsize);
1008                 var uintL numpieces_k1 = ceiling(len2,max_piecelen_k1);
1009                 if (numpieces_k <= 2*numpieces_k1) {
1010                         // keep k
1011                         len2p = max_piecelen_k;
1012                 } else {
1013                         // choose k+1
1014                         k = k+1;
1015                         len2p = max_piecelen_k1;
1016                 }
1017         }
1018         var const uintL l = max_l_table[k];
1019         var const uintL n = k;
1020         var const uintL N = (uintL)1 << n;
1021         CL_ALLOCA_STACK;
1022         var fftr_real* const x = cl_alloc_array(fftr_real,N);
1023         var fftr_real* const y = cl_alloc_array(fftr_real,N);
1024         #ifdef DEBUG_FFTR
1025         var fftr_real* const z = cl_alloc_array(fftr_real,N);
1026         #else
1027         var fftr_real* const z = x; // put z in place of x - saves memory
1028         #endif
1029         var uintD* const tmpprod1 = cl_alloc_array(uintD,len1+1);
1030         var uintL tmpprod_len = floor(l<<n,intDsize)+6;
1031         var uintD* const tmpprod = cl_alloc_array(uintD,tmpprod_len);
1032         var uintL destlen = len1+len2;
1033         clear_loop_lsp(destptr,destlen);
1034         do {
1035                 if (len2p > len2)
1036                         len2p = len2;
1037                 if (len2p == 1) {
1038                         // cheap case
1039                         var uintD* tmpptr = arrayLSDptr(tmpprod1,len1+1);
1040                         mulu_loop_lsp(lspref(sourceptr2,0),sourceptr1,tmpptr,len1);
1041                         if (addto_loop_lsp(tmpptr,destptr,len1+1))
1042                                 if (inc_loop_lsp(destptr lspop (len1+1),destlen-(len1+1)))
1043                                         cl_abort();
1044                 } else {
1045                         var bool squaring = ((sourceptr1 == sourceptr2) && (len1 == len2p));
1046                         // Fill factor x.
1047                         fill_factor(N,x,l,sourceptr1,len1);
1048                         // Fill factor y.
1049                         if (!squaring)
1050                                 fill_factor(N,y,l,sourceptr2,len2p);
1051                         // Multiply.
1052                         if (!squaring)
1053                                 fftr_convolution(n,N, &x[0], &y[0], &z[0]);
1054                         else
1055                                 fftr_convolution(n,N, &x[0], &x[0], &z[0]);
1056                         #ifdef DEBUG_FFTR
1057                         // Check result.
1058                         {
1059                                 var fftr_real re_lo_limit = (fftr_real)(-0.5);
1060                                 var fftr_real re_hi_limit = (fftr_real)N * fftr_pow2_table[l] * fftr_pow2_table[l] + (fftr_real)0.5;
1061                                 for (var uintL i = 0; i < N; i++)
1062                                         if (!(z[i] > re_lo_limit
1063                                               && z[i] < re_hi_limit))
1064                                                 cl_abort();
1065                         }
1066                         #endif
1067                         var uintD* tmpLSDptr = arrayLSDptr(tmpprod,tmpprod_len);
1068                         var uintD* tmpMSDptr = unfill_product(n,N,z,l,tmpLSDptr);
1069                         var uintL tmplen =
1070                           #if CL_DS_BIG_ENDIAN_P
1071                             tmpLSDptr - tmpMSDptr;
1072                           #else
1073                             tmpMSDptr - tmpLSDptr;
1074                           #endif
1075                         if (tmplen > tmpprod_len)
1076                           cl_abort();
1077                         // Add result to destptr[-destlen..-1]:
1078                         if (tmplen > destlen) {
1079                                 if (test_loop_msp(tmpMSDptr,tmplen-destlen))
1080                                         cl_abort();
1081                                 tmplen = destlen;
1082                         }
1083                         if (addto_loop_lsp(tmpLSDptr,destptr,tmplen))
1084                                 if (inc_loop_lsp(destptr lspop tmplen,destlen-tmplen))
1085                                         cl_abort();
1086                 }
1087                 // Decrement len2.
1088                 destptr = destptr lspop len2p;
1089                 destlen -= len2p;
1090                 sourceptr2 = sourceptr2 lspop len2p;
1091                 len2 -= len2p;
1092         } while (len2 > 0);
1093 }
1094
1095 #ifndef _CHECKSUM
1096 #define _CHECKSUM
1097
1098 // Compute a checksum: number mod (2^intDsize-1).
1099 static uintD compute_checksum (const uintD* sourceptr, uintC len)
1100 {
1101         var uintD tmp = ~(uintD)0; // -1-(sum mod 2^intDsize-1), always >0
1102         do {
1103                 var uintD digit = lsprefnext(sourceptr);
1104                 if (digit < tmp)
1105                         tmp -= digit; // subtract digit
1106                 else
1107                         tmp -= digit+1; // subtract digit-(2^intDsize-1)
1108         } while (--len > 0);
1109         return ~tmp;
1110 }
1111
1112 // Multiply two checksums modulo (2^intDsize-1).
1113 static inline uintD multiply_checksum (uintD checksum1, uintD checksum2)
1114 {
1115         var uintD checksum;
1116         var uintD cksum_hi;
1117         #if HAVE_DD
1118         var uintDD cksum = muluD(checksum1,checksum2);
1119         cksum_hi = highD(cksum); checksum = lowD(cksum);
1120         #else
1121         muluD(checksum1,checksum2, cksum_hi =, checksum =);
1122         #endif
1123         if ((checksum += cksum_hi) + 1 <= cksum_hi)
1124                 checksum += 1;
1125         return checksum;
1126 }
1127
1128 #endif // _CHECKSUM
1129
1130 static void mulu_fftr (const uintD* sourceptr1, uintC len1,
1131                        const uintD* sourceptr2, uintC len2,
1132                        uintD* destptr)
1133 {
1134         // Compute checksums of the arguments and multiply them.
1135         var uintD checksum1 = compute_checksum(sourceptr1,len1);
1136         var uintD checksum2 = compute_checksum(sourceptr2,len2);
1137         var uintD checksum = multiply_checksum(checksum1,checksum2);
1138         mulu_fftr_nocheck(sourceptr1,len1,sourceptr2,len2,destptr);
1139         if (!(checksum == compute_checksum(destptr,len1+len2))) {
1140                 fprint(std::cerr, "FFT problem: checksum error\n");
1141                 cl_abort();
1142         }
1143 }