Mercurial > hg > sv-dependency-builds

annotate src/fftw-3.3.5/libbench2/mp.c @ 84:08ae793730bd

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Add null config files
author Chris Cannam
date Mon, 02 Mar 2020 14:03:47 +0000
parents 2cd0e3b3e1fd
children
rev   line source
Chris@42 1 #include "config.h"
Chris@42 2 #include "bench.h"
Chris@42 3 #include <math.h>
Chris@42 4
Chris@42 5 #define DG unsigned short
Chris@42 6 #define ACC unsigned long
Chris@42 7 #define REAL bench_real
Chris@42 8 #define BITS_IN_REAL 53 /* mantissa */
Chris@42 9
Chris@42 10 #define SHFT 16
Chris@42 11 #define RADIX 65536L
Chris@42 12 #define IRADIX (1.0 / RADIX)
Chris@42 13 #define LO(x) ((x) & (RADIX - 1))
Chris@42 14 #define HI(x) ((x) >> SHFT)
Chris@42 15 #define HI_SIGNED(x) \
Chris@42 16 ((((x) + (ACC)(RADIX >> 1) * RADIX) >> SHFT) - (RADIX >> 1))
Chris@42 17 #define ZEROEXP (-32768)
Chris@42 18
Chris@42 19 #define LEN 10
Chris@42 20
Chris@42 21 typedef struct {
Chris@42 22 short sign;
Chris@42 23 short expt;
Chris@42 24 DG d[LEN];
Chris@42 25 } N[1];
Chris@42 26
Chris@42 27 #define EXA a->expt
Chris@42 28 #define EXB b->expt
Chris@42 29 #define EXC c->expt
Chris@42 30
Chris@42 31 #define AD a->d
Chris@42 32 #define BD b->d
Chris@42 33
Chris@42 34 #define SGNA a->sign
Chris@42 35 #define SGNB b->sign
Chris@42 36
Chris@42 37 static const N zero = {{ 1, ZEROEXP, {0} }};
Chris@42 38
Chris@42 39 static void cpy(const N a, N b)
Chris@42 40 {
Chris@42 41 *b = *a;
Chris@42 42 }
Chris@42 43
Chris@42 44 static void fromreal(REAL x, N a)
Chris@42 45 {
Chris@42 46 int i, e;
Chris@42 47
Chris@42 48 cpy(zero, a);
Chris@42 49 if (x == 0.0) return;
Chris@42 50
Chris@42 51 if (x >= 0) { SGNA = 1; }
Chris@42 52 else { SGNA = -1; x = -x; }
Chris@42 53
Chris@42 54 e = 0;
Chris@42 55 while (x >= 1.0) { x *= IRADIX; ++e; }
Chris@42 56 while (x < IRADIX) { x *= RADIX; --e; }
Chris@42 57 EXA = e;
Chris@42 58
Chris@42 59 for (i = LEN - 1; i >= 0 && x != 0.0; --i) {
Chris@42 60 REAL y;
Chris@42 61
Chris@42 62 x *= RADIX;
Chris@42 63 y = (REAL) ((int) x);
Chris@42 64 AD[i] = (DG)y;
Chris@42 65 x -= y;
Chris@42 66 }
Chris@42 67 }
Chris@42 68
Chris@42 69 static void fromshort(int x, N a)
Chris@42 70 {
Chris@42 71 cpy(zero, a);
Chris@42 72
Chris@42 73 if (x < 0) { x = -x; SGNA = -1; }
Chris@42 74 else { SGNA = 1; }
Chris@42 75 EXA = 1;
Chris@42 76 AD[LEN - 1] = x;
Chris@42 77 }
Chris@42 78
Chris@42 79 static void pack(DG *d, int e, int s, int l, N a)
Chris@42 80 {
Chris@42 81 int i, j;
Chris@42 82
Chris@42 83 for (i = l - 1; i >= 0; --i, --e)
Chris@42 84 if (d[i] != 0)
Chris@42 85 break;
Chris@42 86
Chris@42 87 if (i < 0) {
Chris@42 88 /* number is zero */
Chris@42 89 cpy(zero, a);
Chris@42 90 } else {
Chris@42 91 EXA = e;
Chris@42 92 SGNA = s;
Chris@42 93
Chris@42 94 if (i >= LEN - 1) {
Chris@42 95 for (j = LEN - 1; j >= 0; --i, --j)
Chris@42 96 AD[j] = d[i];
Chris@42 97 } else {
Chris@42 98 for (j = LEN - 1; i >= 0; --i, --j)
Chris@42 99 AD[j] = d[i];
Chris@42 100 for ( ; j >= 0; --j)
Chris@42 101 AD[j] = 0;
Chris@42 102 }
Chris@42 103 }
Chris@42 104 }
Chris@42 105
Chris@42 106
Chris@42 107 /* compare absolute values */
Chris@42 108 static int abscmp(const N a, const N b)
Chris@42 109 {
Chris@42 110 int i;
Chris@42 111 if (EXA > EXB) return 1;
Chris@42 112 if (EXA < EXB) return -1;
Chris@42 113 for (i = LEN - 1; i >= 0; --i) {
Chris@42 114 if (AD[i] > BD[i])
Chris@42 115 return 1;
Chris@42 116 if (AD[i] < BD[i])
Chris@42 117 return -1;
Chris@42 118 }
Chris@42 119 return 0;
Chris@42 120 }
Chris@42 121
Chris@42 122 static int eq(const N a, const N b)
Chris@42 123 {
Chris@42 124 return (SGNA == SGNB) && (abscmp(a, b) == 0);
Chris@42 125 }
Chris@42 126
Chris@42 127 /* add magnitudes, for |a| >= |b| */
Chris@42 128 static void addmag0(int s, const N a, const N b, N c)
Chris@42 129 {
Chris@42 130 int ia, ib;
Chris@42 131 ACC r = 0;
Chris@42 132 DG d[LEN + 1];
Chris@42 133
Chris@42 134 for (ia = 0, ib = EXA - EXB; ib < LEN; ++ia, ++ib) {
Chris@42 135 r += (ACC)AD[ia] + (ACC)BD[ib];
Chris@42 136 d[ia] = LO(r);
Chris@42 137 r = HI(r);
Chris@42 138 }
Chris@42 139 for (; ia < LEN; ++ia) {
Chris@42 140 r += (ACC)AD[ia];
Chris@42 141 d[ia] = LO(r);
Chris@42 142 r = HI(r);
Chris@42 143 }
Chris@42 144 d[ia] = LO(r);
Chris@42 145 pack(d, EXA + 1, s * SGNA, LEN + 1, c);
Chris@42 146 }
Chris@42 147
Chris@42 148 static void addmag(int s, const N a, const N b, N c)
Chris@42 149 {
Chris@42 150 if (abscmp(a, b) > 0) addmag0(1, a, b, c); else addmag0(s, b, a, c);
Chris@42 151 }
Chris@42 152
Chris@42 153 /* subtract magnitudes, for |a| >= |b| */
Chris@42 154 static void submag0(int s, const N a, const N b, N c)
Chris@42 155 {
Chris@42 156 int ia, ib;
Chris@42 157 ACC r = 0;
Chris@42 158 DG d[LEN];
Chris@42 159
Chris@42 160 for (ia = 0, ib = EXA - EXB; ib < LEN; ++ia, ++ib) {
Chris@42 161 r += (ACC)AD[ia] - (ACC)BD[ib];
Chris@42 162 d[ia] = LO(r);
Chris@42 163 r = HI_SIGNED(r);
Chris@42 164 }
Chris@42 165 for (; ia < LEN; ++ia) {
Chris@42 166 r += (ACC)AD[ia];
Chris@42 167 d[ia] = LO(r);
Chris@42 168 r = HI_SIGNED(r);
Chris@42 169 }
Chris@42 170
Chris@42 171 pack(d, EXA, s * SGNA, LEN, c);
Chris@42 172 }
Chris@42 173
Chris@42 174 static void submag(int s, const N a, const N b, N c)
Chris@42 175 {
Chris@42 176 if (abscmp(a, b) > 0) submag0(1, a, b, c); else submag0(s, b, a, c);
Chris@42 177 }
Chris@42 178
Chris@42 179 /* c = a + b */
Chris@42 180 static void add(const N a, const N b, N c)
Chris@42 181 {
Chris@42 182 if (SGNA == SGNB) addmag(1, a, b, c); else submag(1, a, b, c);
Chris@42 183 }
Chris@42 184
Chris@42 185 static void sub(const N a, const N b, N c)
Chris@42 186 {
Chris@42 187 if (SGNA == SGNB) submag(-1, a, b, c); else addmag(-1, a, b, c);
Chris@42 188 }
Chris@42 189
Chris@42 190 static void mul(const N a, const N b, N c)
Chris@42 191 {
Chris@42 192 DG d[2 * LEN];
Chris@42 193 int i, j, k;
Chris@42 194 ACC r;
Chris@42 195
Chris@42 196 for (i = 0; i < LEN; ++i)
Chris@42 197 d[2 * i] = d[2 * i + 1] = 0;
Chris@42 198
Chris@42 199 for (i = 0; i < LEN; ++i) {
Chris@42 200 ACC ai = AD[i];
Chris@42 201 if (ai) {
Chris@42 202 r = 0;
Chris@42 203 for (j = 0, k = i; j < LEN; ++j, ++k) {
Chris@42 204 r += ai * (ACC)BD[j] + (ACC)d[k];
Chris@42 205 d[k] = LO(r);
Chris@42 206 r = HI(r);
Chris@42 207 }
Chris@42 208 d[k] = LO(r);
Chris@42 209 }
Chris@42 210 }
Chris@42 211
Chris@42 212 pack(d, EXA + EXB, SGNA * SGNB, 2 * LEN, c);
Chris@42 213 }
Chris@42 214
Chris@42 215 static REAL toreal(const N a)
Chris@42 216 {
Chris@42 217 REAL h, l, f;
Chris@42 218 int i, bits;
Chris@42 219 ACC r;
Chris@42 220 DG sticky;
Chris@42 221
Chris@42 222 if (EXA != ZEROEXP) {
Chris@42 223 f = IRADIX;
Chris@42 224 i = LEN;
Chris@42 225
Chris@42 226 bits = 0;
Chris@42 227 h = (r = AD[--i]) * f; f *= IRADIX;
Chris@42 228 for (bits = 0; r > 0; ++bits)
Chris@42 229 r >>= 1;
Chris@42 230
Chris@42 231 /* first digit */
Chris@42 232 while (bits + SHFT <= BITS_IN_REAL) {
Chris@42 233 h += AD[--i] * f; f *= IRADIX; bits += SHFT;
Chris@42 234 }
Chris@42 235
Chris@42 236 /* guard digit (leave one bit for sticky bit, hence `<' instead
Chris@42 237 of `<=') */
Chris@42 238 bits = 0; l = 0.0;
Chris@42 239 while (bits + SHFT < BITS_IN_REAL) {
Chris@42 240 l += AD[--i] * f; f *= IRADIX; bits += SHFT;
Chris@42 241 }
Chris@42 242
Chris@42 243 /* sticky bit */
Chris@42 244 sticky = 0;
Chris@42 245 while (i > 0)
Chris@42 246 sticky |= AD[--i];
Chris@42 247
Chris@42 248 if (sticky)
Chris@42 249 l += (RADIX / 2) * f;
Chris@42 250
Chris@42 251 h += l;
Chris@42 252
Chris@42 253 for (i = 0; i < EXA; ++i) h *= (REAL)RADIX;
Chris@42 254 for (i = 0; i > EXA; --i) h *= IRADIX;
Chris@42 255 if (SGNA == -1) h = -h;
Chris@42 256 return h;
Chris@42 257 } else {
Chris@42 258 return 0.0;
Chris@42 259 }
Chris@42 260 }
Chris@42 261
Chris@42 262 static void neg(N a)
Chris@42 263 {
Chris@42 264 SGNA = -SGNA;
Chris@42 265 }
Chris@42 266
Chris@42 267 static void inv(const N a, N x)
Chris@42 268 {
Chris@42 269 N w, z, one, two;
Chris@42 270
Chris@42 271 fromreal(1.0 / toreal(a), x); /* initial guess */
Chris@42 272 fromshort(1, one);
Chris@42 273 fromshort(2, two);
Chris@42 274
Chris@42 275 for (;;) {
Chris@42 276 /* Newton */
Chris@42 277 mul(a, x, w);
Chris@42 278 sub(two, w, z);
Chris@42 279 if (eq(one, z)) break;
Chris@42 280 mul(x, z, x);
Chris@42 281 }
Chris@42 282 }
Chris@42 283
Chris@42 284
Chris@42 285 /* 2 pi */
Chris@42 286 static const N n2pi = {{
Chris@42 287 1, 1,
Chris@42 288 {18450, 59017, 1760, 5212, 9779, 4518, 2886, 54545, 18558, 6}
Chris@42 289 }};
Chris@42 290
Chris@42 291 /* 1 / 31! */
Chris@42 292 static const N i31fac = {{
Chris@42 293 1, -7,
Chris@42 294 {28087, 45433, 51357, 24545, 14291, 3954, 57879, 8109, 38716, 41382}
Chris@42 295 }};
Chris@42 296
Chris@42 297
Chris@42 298 /* 1 / 32! */
Chris@42 299 static const N i32fac = {{
Chris@42 300 1, -7,
Chris@42 301 {52078, 60811, 3652, 39679, 37310, 47227, 28432, 57597, 13497, 1293}
Chris@42 302 }};
Chris@42 303
Chris@42 304 static void msin(const N a, N b)
Chris@42 305 {
Chris@42 306 N a2, g, k;
Chris@42 307 int i;
Chris@42 308
Chris@42 309 cpy(i31fac, g);
Chris@42 310 cpy(g, b);
Chris@42 311 mul(a, a, a2);
Chris@42 312
Chris@42 313 /* Taylor */
Chris@42 314 for (i = 31; i > 1; i -= 2) {
Chris@42 315 fromshort(i * (i - 1), k);
Chris@42 316 mul(k, g, g);
Chris@42 317 mul(a2, b, k);
Chris@42 318 sub(g, k, b);
Chris@42 319 }
Chris@42 320 mul(a, b, b);
Chris@42 321 }
Chris@42 322
Chris@42 323 static void mcos(const N a, N b)
Chris@42 324 {
Chris@42 325 N a2, g, k;
Chris@42 326 int i;
Chris@42 327
Chris@42 328 cpy(i32fac, g);
Chris@42 329 cpy(g, b);
Chris@42 330 mul(a, a, a2);
Chris@42 331
Chris@42 332 /* Taylor */
Chris@42 333 for (i = 32; i > 0; i -= 2) {
Chris@42 334 fromshort(i * (i - 1), k);
Chris@42 335 mul(k, g, g);
Chris@42 336 mul(a2, b, k);
Chris@42 337 sub(g, k, b);
Chris@42 338 }
Chris@42 339 }
Chris@42 340
Chris@42 341 static void by2pi(REAL m, REAL n, N a)
Chris@42 342 {
Chris@42 343 N b;
Chris@42 344
Chris@42 345 fromreal(n, b);
Chris@42 346 inv(b, a);
Chris@42 347 fromreal(m, b);
Chris@42 348 mul(a, b, a);
Chris@42 349 mul(n2pi, a, a);
Chris@42 350 }
Chris@42 351
Chris@42 352 static void sin2pi(REAL m, REAL n, N a);
Chris@42 353 static void cos2pi(REAL m, REAL n, N a)
Chris@42 354 {
Chris@42 355 N b;
Chris@42 356 if (m < 0) cos2pi(-m, n, a);
Chris@42 357 else if (m > n * 0.5) cos2pi(n - m, n, a);
Chris@42 358 else if (m > n * 0.25) {sin2pi(m - n * 0.25, n, a); neg(a);}
Chris@42 359 else if (m > n * 0.125) sin2pi(n * 0.25 - m, n, a);
Chris@42 360 else { by2pi(m, n, b); mcos(b, a); }
Chris@42 361 }
Chris@42 362
Chris@42 363 static void sin2pi(REAL m, REAL n, N a)
Chris@42 364 {
Chris@42 365 N b;
Chris@42 366 if (m < 0) {sin2pi(-m, n, a); neg(a);}
Chris@42 367 else if (m > n * 0.5) {sin2pi(n - m, n, a); neg(a);}
Chris@42 368 else if (m > n * 0.25) {cos2pi(m - n * 0.25, n, a);}
Chris@42 369 else if (m > n * 0.125) {cos2pi(n * 0.25 - m, n, a);}
Chris@42 370 else {by2pi(m, n, b); msin(b, a);}
Chris@42 371 }
Chris@42 372
Chris@42 373 /*----------------------------------------------------------------------*/
Chris@42 374 /* FFT stuff */
Chris@42 375
Chris@42 376 /* (r0 + i i0)(r1 + i i1) */
Chris@42 377 static void cmul(N r0, N i0, N r1, N i1, N r2, N i2)
Chris@42 378 {
Chris@42 379 N s, t, q;
Chris@42 380 mul(r0, r1, s);
Chris@42 381 mul(i0, i1, t);
Chris@42 382 sub(s, t, q);
Chris@42 383 mul(r0, i1, s);
Chris@42 384 mul(i0, r1, t);
Chris@42 385 add(s, t, i2);
Chris@42 386 cpy(q, r2);
Chris@42 387 }
Chris@42 388
Chris@42 389 /* (r0 - i i0)(r1 + i i1) */
Chris@42 390 static void cmulj(N r0, N i0, N r1, N i1, N r2, N i2)
Chris@42 391 {
Chris@42 392 N s, t, q;
Chris@42 393 mul(r0, r1, s);
Chris@42 394 mul(i0, i1, t);
Chris@42 395 add(s, t, q);
Chris@42 396 mul(r0, i1, s);
Chris@42 397 mul(i0, r1, t);
Chris@42 398 sub(s, t, i2);
Chris@42 399 cpy(q, r2);
Chris@42 400 }
Chris@42 401
Chris@42 402 static void mcexp(int m, int n, N r, N i)
Chris@42 403 {
Chris@42 404 static int cached_n = -1;
Chris@42 405 static N w[64][2];
Chris@42 406 int k, j;
Chris@42 407 if (n != cached_n) {
Chris@42 408 for (j = 1, k = 0; j < n; j += j, ++k) {
Chris@42 409 cos2pi(j, n, w[k][0]);
Chris@42 410 sin2pi(j, n, w[k][1]);
Chris@42 411 }
Chris@42 412 cached_n = n;
Chris@42 413 }
Chris@42 414
Chris@42 415 fromshort(1, r);
Chris@42 416 fromshort(0, i);
Chris@42 417 if (m > 0) {
Chris@42 418 for (k = 0; m; ++k, m >>= 1)
Chris@42 419 if (m & 1)
Chris@42 420 cmul(w[k][0], w[k][1], r, i, r, i);
Chris@42 421 } else {
Chris@42 422 m = -m;
Chris@42 423 for (k = 0; m; ++k, m >>= 1)
Chris@42 424 if (m & 1)
Chris@42 425 cmulj(w[k][0], w[k][1], r, i, r, i);
Chris@42 426 }
Chris@42 427 }
Chris@42 428
Chris@42 429 static void bitrev(int n, N *a)
Chris@42 430 {
Chris@42 431 int i, j, m;
Chris@42 432 for (i = j = 0; i < n - 1; ++i) {
Chris@42 433 if (i < j) {
Chris@42 434 N t;
Chris@42 435 cpy(a[2*i], t); cpy(a[2*j], a[2*i]); cpy(t, a[2*j]);
Chris@42 436 cpy(a[2*i+1], t); cpy(a[2*j+1], a[2*i+1]); cpy(t, a[2*j+1]);
Chris@42 437 }
Chris@42 438
Chris@42 439 /* bit reversed counter */
Chris@42 440 m = n; do { m >>= 1; j ^= m; } while (!(j & m));
Chris@42 441 }
Chris@42 442 }
Chris@42 443
Chris@42 444 static void fft0(int n, N *a, int sign)
Chris@42 445 {
Chris@42 446 int i, j, k;
Chris@42 447
Chris@42 448 bitrev(n, a);
Chris@42 449 for (i = 1; i < n; i = 2 * i) {
Chris@42 450 for (j = 0; j < i; ++j) {
Chris@42 451 N wr, wi;
Chris@42 452 mcexp(sign * (int)j, 2 * i, wr, wi);
Chris@42 453 for (k = j; k < n; k += 2 * i) {
Chris@42 454 N *a0 = a + 2 * k;
Chris@42 455 N *a1 = a0 + 2 * i;
Chris@42 456 N r0, i0, r1, i1, t0, t1, xr, xi;
Chris@42 457 cpy(a0[0], r0); cpy(a0[1], i0);
Chris@42 458 cpy(a1[0], r1); cpy(a1[1], i1);
Chris@42 459 mul(r1, wr, t0); mul(i1, wi, t1); sub(t0, t1, xr);
Chris@42 460 mul(r1, wi, t0); mul(i1, wr, t1); add(t0, t1, xi);
Chris@42 461 add(r0, xr, a0[0]); add(i0, xi, a0[1]);
Chris@42 462 sub(r0, xr, a1[0]); sub(i0, xi, a1[1]);
Chris@42 463 }
Chris@42 464 }
Chris@42 465 }
Chris@42 466 }
Chris@42 467
Chris@42 468 /* a[2*k]+i*a[2*k+1] = exp(2*pi*i*k^2/(2*n)) */
Chris@42 469 static void bluestein_sequence(int n, N *a)
Chris@42 470 {
Chris@42 471 int k, ksq, n2 = 2 * n;
Chris@42 472
Chris@42 473 ksq = 1; /* (-1)^2 */
Chris@42 474 for (k = 0; k < n; ++k) {
Chris@42 475 /* careful with overflow */
Chris@42 476 ksq = ksq + 2*k - 1; while (ksq > n2) ksq -= n2;
Chris@42 477 mcexp(ksq, n2, a[2*k], a[2*k+1]);
Chris@42 478 }
Chris@42 479 }
Chris@42 480
Chris@42 481 static int pow2_atleast(int x)
Chris@42 482 {
Chris@42 483 int h;
Chris@42 484 for (h = 1; h < x; h = 2 * h)
Chris@42 485 ;
Chris@42 486 return h;
Chris@42 487 }
Chris@42 488
Chris@42 489 static N *cached_bluestein_w = 0;
Chris@42 490 static N *cached_bluestein_y = 0;
Chris@42 491 static int cached_bluestein_n = -1;
Chris@42 492
Chris@42 493 static void bluestein(int n, N *a)
Chris@42 494 {
Chris@42 495 int nb = pow2_atleast(2 * n);
Chris@42 496 N *b = (N *)bench_malloc(2 * nb * sizeof(N));
Chris@42 497 N *w = cached_bluestein_w;
Chris@42 498 N *y = cached_bluestein_y;
Chris@42 499 N nbinv;
Chris@42 500 int i;
Chris@42 501
Chris@42 502 fromreal(1.0 / nb, nbinv); /* exact because nb = 2^k */
Chris@42 503
Chris@42 504 if (cached_bluestein_n != n) {
Chris@42 505 if (w) bench_free(w);
Chris@42 506 if (y) bench_free(y);
Chris@42 507 w = (N *)bench_malloc(2 * n * sizeof(N));
Chris@42 508 y = (N *)bench_malloc(2 * nb * sizeof(N));
Chris@42 509 cached_bluestein_n = n;
Chris@42 510 cached_bluestein_w = w;
Chris@42 511 cached_bluestein_y = y;
Chris@42 512
Chris@42 513 bluestein_sequence(n, w);
Chris@42 514 for (i = 0; i < 2*nb; ++i) cpy(zero, y[i]);
Chris@42 515
Chris@42 516 for (i = 0; i < n; ++i) {
Chris@42 517 cpy(w[2*i], y[2*i]);
Chris@42 518 cpy(w[2*i+1], y[2*i+1]);
Chris@42 519 }
Chris@42 520 for (i = 1; i < n; ++i) {
Chris@42 521 cpy(w[2*i], y[2*(nb-i)]);
Chris@42 522 cpy(w[2*i+1], y[2*(nb-i)+1]);
Chris@42 523 }
Chris@42 524
Chris@42 525 fft0(nb, y, -1);
Chris@42 526 }
Chris@42 527
Chris@42 528 for (i = 0; i < 2*nb; ++i) cpy(zero, b[i]);
Chris@42 529
Chris@42 530 for (i = 0; i < n; ++i)
Chris@42 531 cmulj(w[2*i], w[2*i+1], a[2*i], a[2*i+1], b[2*i], b[2*i+1]);
Chris@42 532
Chris@42 533 /* scaled convolution b * y */
Chris@42 534 fft0(nb, b, -1);
Chris@42 535
Chris@42 536 for (i = 0; i < nb; ++i)
Chris@42 537 cmul(b[2*i], b[2*i+1], y[2*i], y[2*i+1], b[2*i], b[2*i+1]);
Chris@42 538 fft0(nb, b, 1);
Chris@42 539
Chris@42 540 for (i = 0; i < n; ++i) {
Chris@42 541 cmulj(w[2*i], w[2*i+1], b[2*i], b[2*i+1], a[2*i], a[2*i+1]);
Chris@42 542 mul(nbinv, a[2*i], a[2*i]);
Chris@42 543 mul(nbinv, a[2*i+1], a[2*i+1]);
Chris@42 544 }
Chris@42 545
Chris@42 546 bench_free(b);
Chris@42 547 }
Chris@42 548
Chris@42 549 static void swapri(int n, N *a)
Chris@42 550 {
Chris@42 551 int i;
Chris@42 552 for (i = 0; i < n; ++i) {
Chris@42 553 N t;
Chris@42 554 cpy(a[2 * i], t);
Chris@42 555 cpy(a[2 * i + 1], a[2 * i]);
Chris@42 556 cpy(t, a[2 * i + 1]);
Chris@42 557 }
Chris@42 558 }
Chris@42 559
Chris@42 560 static void fft1(int n, N *a, int sign)
Chris@42 561 {
Chris@42 562 if (power_of_two(n)) {
Chris@42 563 fft0(n, a, sign);
Chris@42 564 } else {
Chris@42 565 if (sign == 1) swapri(n, a);
Chris@42 566 bluestein(n, a);
Chris@42 567 if (sign == 1) swapri(n, a);
Chris@42 568 }
Chris@42 569 }
Chris@42 570
Chris@42 571 static void fromrealv(int n, bench_complex *a, N *b)
Chris@42 572 {
Chris@42 573 int i;
Chris@42 574
Chris@42 575 for (i = 0; i < n; ++i) {
Chris@42 576 fromreal(c_re(a[i]), b[2 * i]);
Chris@42 577 fromreal(c_im(a[i]), b[2 * i + 1]);
Chris@42 578 }
Chris@42 579 }
Chris@42 580
Chris@42 581 static void compare(int n, N *a, N *b, double *err)
Chris@42 582 {
Chris@42 583 int i;
Chris@42 584 double e1, e2, einf;
Chris@42 585 double n1, n2, ninf;
Chris@42 586
Chris@42 587 e1 = e2 = einf = 0.0;
Chris@42 588 n1 = n2 = ninf = 0.0;
Chris@42 589
Chris@42 590 # define DO(x1, x2, xinf, var) { \
Chris@42 591 double d = var; \
Chris@42 592 if (d < 0) d = -d; \
Chris@42 593 x1 += d; x2 += d * d; if (d > xinf) xinf = d; \
Chris@42 594 }
Chris@42 595
Chris@42 596 for (i = 0; i < 2 * n; ++i) {
Chris@42 597 N dd;
Chris@42 598 sub(a[i], b[i], dd);
Chris@42 599 DO(n1, n2, ninf, toreal(a[i]));
Chris@42 600 DO(e1, e2, einf, toreal(dd));
Chris@42 601 }
Chris@42 602
Chris@42 603 # undef DO
Chris@42 604 err[0] = e1 / n1;
Chris@42 605 err[1] = sqrt(e2 / n2);
Chris@42 606 err[2] = einf / ninf;
Chris@42 607 }
Chris@42 608
Chris@42 609 void fftaccuracy(int n, bench_complex *a, bench_complex *ffta,
Chris@42 610 int sign, double err[6])
Chris@42 611 {
Chris@42 612 N *b = (N *)bench_malloc(2 * n * sizeof(N));
Chris@42 613 N *fftb = (N *)bench_malloc(2 * n * sizeof(N));
Chris@42 614 N mn, ninv;
Chris@42 615 int i;
Chris@42 616
Chris@42 617 fromreal(n, mn); inv(mn, ninv);
Chris@42 618
Chris@42 619 /* forward error */
Chris@42 620 fromrealv(n, a, b); fromrealv(n, ffta, fftb);
Chris@42 621 fft1(n, b, sign);
Chris@42 622 compare(n, b, fftb, err);
Chris@42 623
Chris@42 624 /* backward error */
Chris@42 625 fromrealv(n, a, b); fromrealv(n, ffta, fftb);
Chris@42 626 for (i = 0; i < 2 * n; ++i) mul(fftb[i], ninv, fftb[i]);
Chris@42 627 fft1(n, fftb, -sign);
Chris@42 628 compare(n, b, fftb, err + 3);
Chris@42 629
Chris@42 630 bench_free(fftb);
Chris@42 631 bench_free(b);
Chris@42 632 }
Chris@42 633
Chris@42 634 void fftaccuracy_done(void)
Chris@42 635 {
Chris@42 636 if (cached_bluestein_w) bench_free(cached_bluestein_w);
Chris@42 637 if (cached_bluestein_y) bench_free(cached_bluestein_y);
Chris@42 638 cached_bluestein_w = 0;
Chris@42 639 cached_bluestein_y = 0;
Chris@42 640 cached_bluestein_n = -1;
Chris@42 641 }