src/base/digitseq/cl_DS_mul_fftcs.h

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