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