cannam@95: #include "config.h"
cannam@95: #include "bench.h"
cannam@95: #include <math.h>
cannam@95: 
cannam@95: #define DG unsigned short
cannam@95: #define ACC unsigned long
cannam@95: #define REAL bench_real
cannam@95: #define BITS_IN_REAL 53 /* mantissa */
cannam@95: 
cannam@95: #define SHFT 16
cannam@95: #define RADIX 65536L
cannam@95: #define IRADIX (1.0 / RADIX)
cannam@95: #define LO(x) ((x) & (RADIX - 1))
cannam@95: #define HI(x) ((x) >> SHFT)
cannam@95: #define HI_SIGNED(x) \
cannam@95:    ((((x) + (ACC)(RADIX >> 1) * RADIX) >> SHFT) - (RADIX >> 1))
cannam@95: #define ZEROEXP (-32768)
cannam@95: 
cannam@95: #define LEN 10
cannam@95: 
cannam@95: typedef struct {
cannam@95:      short sign;
cannam@95:      short expt;
cannam@95:      DG d[LEN]; 
cannam@95: } N[1];
cannam@95: 
cannam@95: #define EXA a->expt
cannam@95: #define EXB b->expt
cannam@95: #define EXC c->expt
cannam@95: 
cannam@95: #define AD a->d
cannam@95: #define BD b->d
cannam@95: 
cannam@95: #define SGNA a->sign
cannam@95: #define SGNB b->sign
cannam@95: 
cannam@95: static const N zero = {{ 1, ZEROEXP, {0} }};
cannam@95: 
cannam@95: static void cpy(const N a, N b)
cannam@95: {
cannam@95:      *b = *a;
cannam@95: }
cannam@95: 
cannam@95: static void fromreal(REAL x, N a)
cannam@95: {
cannam@95:      int i, e;
cannam@95: 
cannam@95:      cpy(zero, a);
cannam@95:      if (x == 0.0) return;
cannam@95:      
cannam@95:      if (x >= 0) { SGNA = 1; }
cannam@95:      else       { SGNA = -1; x = -x; }
cannam@95: 
cannam@95:      e = 0;
cannam@95:      while (x >= 1.0) { x *= IRADIX; ++e; }
cannam@95:      while (x < IRADIX) { x *= RADIX; --e; }
cannam@95:      EXA = e;
cannam@95:      
cannam@95:      for (i = LEN - 1; i >= 0 && x != 0.0; --i) {
cannam@95: 	  REAL y;
cannam@95: 
cannam@95: 	  x *= RADIX;
cannam@95: 	  y = (REAL) ((int) x);
cannam@95: 	  AD[i] = (DG)y;
cannam@95: 	  x -= y;
cannam@95:      }
cannam@95: }
cannam@95: 
cannam@95: static void fromshort(int x, N a)
cannam@95: {
cannam@95:      cpy(zero, a);
cannam@95: 
cannam@95:      if (x < 0) { x = -x; SGNA = -1; } 
cannam@95:      else { SGNA = 1; }
cannam@95:      EXA = 1;
cannam@95:      AD[LEN - 1] = x;
cannam@95: }
cannam@95: 
cannam@95: static void pack(DG *d, int e, int s, int l, N a)
cannam@95: {
cannam@95:      int i, j;
cannam@95: 
cannam@95:      for (i = l - 1; i >= 0; --i, --e) 
cannam@95: 	  if (d[i] != 0) 
cannam@95: 	       break;
cannam@95: 
cannam@95:      if (i < 0) {
cannam@95: 	  /* number is zero */
cannam@95: 	  cpy(zero, a);
cannam@95:      } else {
cannam@95: 	  EXA = e;
cannam@95: 	  SGNA = s;
cannam@95: 
cannam@95: 	  if (i >= LEN - 1) {
cannam@95: 	       for (j = LEN - 1; j >= 0; --i, --j)
cannam@95: 		    AD[j] = d[i];
cannam@95: 	  } else {
cannam@95: 	       for (j = LEN - 1; i >= 0; --i, --j)
cannam@95: 		    AD[j] = d[i];
cannam@95: 	       for ( ; j >= 0; --j)
cannam@95: 		    AD[j] = 0;
cannam@95: 	  }
cannam@95:      }
cannam@95: }
cannam@95: 
cannam@95: 
cannam@95: /* compare absolute values */
cannam@95: static int abscmp(const N a, const N b)
cannam@95: {
cannam@95:      int i;
cannam@95:      if (EXA > EXB) return 1;
cannam@95:      if (EXA < EXB) return -1;
cannam@95:      for (i = LEN - 1; i >= 0; --i) {
cannam@95: 	  if (AD[i] > BD[i])
cannam@95: 	       return 1;
cannam@95: 	  if (AD[i] < BD[i])
cannam@95: 	       return -1;
cannam@95:      }
cannam@95:      return 0;
cannam@95: }
cannam@95: 
cannam@95: static int eq(const N a, const N b)
cannam@95: {
cannam@95:      return (SGNA == SGNB) && (abscmp(a, b) == 0);
cannam@95: }
cannam@95: 
cannam@95: /* add magnitudes, for |a| >= |b| */
cannam@95: static void addmag0(int s, const N a, const N b, N c)
cannam@95: {
cannam@95:      int ia, ib;
cannam@95:      ACC r = 0;
cannam@95:      DG d[LEN + 1];
cannam@95: 
cannam@95:      for (ia = 0, ib = EXA - EXB; ib < LEN; ++ia, ++ib) {
cannam@95: 	  r += (ACC)AD[ia] + (ACC)BD[ib];
cannam@95: 	  d[ia] = LO(r);
cannam@95: 	  r = HI(r);
cannam@95:      }     
cannam@95:      for (; ia < LEN; ++ia) {
cannam@95: 	  r += (ACC)AD[ia];
cannam@95: 	  d[ia] = LO(r);
cannam@95: 	  r = HI(r);
cannam@95:      }
cannam@95:      d[ia] = LO(r);
cannam@95:      pack(d, EXA + 1, s * SGNA, LEN + 1, c);
cannam@95: }
cannam@95: 
cannam@95: static void addmag(int s, const N a, const N b, N c)
cannam@95: {
cannam@95:      if (abscmp(a, b) > 0) addmag0(1, a, b, c); else addmag0(s, b, a, c);
cannam@95: }
cannam@95: 
cannam@95: /* subtract magnitudes, for |a| >= |b| */
cannam@95: static void submag0(int s, const N a, const N b, N c)
cannam@95: {
cannam@95:      int ia, ib;
cannam@95:      ACC r = 0;
cannam@95:      DG d[LEN];
cannam@95: 
cannam@95:      for (ia = 0, ib = EXA - EXB; ib < LEN; ++ia, ++ib) {
cannam@95: 	  r += (ACC)AD[ia] - (ACC)BD[ib];
cannam@95: 	  d[ia] = LO(r);
cannam@95: 	  r = HI_SIGNED(r);
cannam@95:      }     
cannam@95:      for (; ia < LEN; ++ia) {
cannam@95: 	  r += (ACC)AD[ia];
cannam@95: 	  d[ia] = LO(r);
cannam@95: 	  r = HI_SIGNED(r);
cannam@95:      }
cannam@95: 
cannam@95:      pack(d, EXA, s * SGNA, LEN, c);
cannam@95: }
cannam@95: 
cannam@95: static void submag(int s, const N a, const N b, N c)
cannam@95: {
cannam@95:      if (abscmp(a, b) > 0) submag0(1, a, b, c); else submag0(s, b, a, c);
cannam@95: }
cannam@95: 
cannam@95: /* c = a + b */
cannam@95: static void add(const N a, const N b, N c)
cannam@95: {
cannam@95:      if (SGNA == SGNB) addmag(1, a, b, c); else submag(1, a, b, c);
cannam@95: }
cannam@95: 
cannam@95: static void sub(const N a, const N b, N c)
cannam@95: {
cannam@95:      if (SGNA == SGNB) submag(-1, a, b, c); else addmag(-1, a, b, c);
cannam@95: }
cannam@95: 
cannam@95: static void mul(const N a, const N b, N c)
cannam@95: {
cannam@95:      DG d[2 * LEN];
cannam@95:      int i, j, k;
cannam@95:      ACC r;
cannam@95: 
cannam@95:      for (i = 0; i < LEN; ++i)
cannam@95: 	  d[2 * i] = d[2 * i + 1] = 0;
cannam@95: 
cannam@95:      for (i = 0; i < LEN; ++i) {
cannam@95: 	  ACC ai = AD[i];
cannam@95: 	  if (ai) {
cannam@95: 	       r = 0;
cannam@95: 	       for (j = 0, k = i; j < LEN; ++j, ++k) {
cannam@95: 		    r += ai * (ACC)BD[j] + (ACC)d[k];
cannam@95: 		    d[k] = LO(r);
cannam@95: 		    r = HI(r);
cannam@95: 	       }
cannam@95: 	       d[k] = LO(r);
cannam@95: 	  }
cannam@95:      }
cannam@95: 
cannam@95:      pack(d, EXA + EXB, SGNA * SGNB, 2 * LEN, c);
cannam@95: }
cannam@95: 
cannam@95: static REAL toreal(const N a)
cannam@95: {
cannam@95:      REAL h, l, f;
cannam@95:      int i, bits;
cannam@95:      ACC r;
cannam@95:      DG sticky;
cannam@95: 
cannam@95:      if (EXA != ZEROEXP) {
cannam@95: 	  f = IRADIX;
cannam@95: 	  i = LEN;
cannam@95: 
cannam@95: 	  bits = 0;
cannam@95: 	  h = (r = AD[--i]) * f; f *= IRADIX;
cannam@95: 	  for (bits = 0; r > 0; ++bits)
cannam@95: 	       r >>= 1;
cannam@95: 
cannam@95: 	  /* first digit */
cannam@95: 	  while (bits + SHFT <= BITS_IN_REAL) {
cannam@95: 	       h += AD[--i] * f;  f *= IRADIX; bits += SHFT;
cannam@95: 	  }
cannam@95: 
cannam@95: 	  /* guard digit (leave one bit for sticky bit, hence `<' instead
cannam@95: 	     of `<=') */
cannam@95: 	  bits = 0; l = 0.0;
cannam@95: 	  while (bits + SHFT < BITS_IN_REAL) {
cannam@95: 	       l += AD[--i] * f;  f *= IRADIX; bits += SHFT;
cannam@95: 	  }
cannam@95: 	  
cannam@95: 	  /* sticky bit */
cannam@95: 	  sticky = 0;
cannam@95: 	  while (i > 0) 
cannam@95: 	       sticky |= AD[--i];
cannam@95: 
cannam@95: 	  if (sticky)
cannam@95: 	       l += (RADIX / 2) * f;
cannam@95: 
cannam@95: 	  h += l;
cannam@95: 
cannam@95: 	  for (i = 0; i < EXA; ++i) h *= (REAL)RADIX;
cannam@95: 	  for (i = 0; i > EXA; --i) h *= IRADIX;
cannam@95: 	  if (SGNA == -1) h = -h;
cannam@95: 	  return h;
cannam@95:      } else {
cannam@95: 	  return 0.0;
cannam@95:      }
cannam@95: }
cannam@95: 
cannam@95: static void neg(N a)
cannam@95: {
cannam@95:      SGNA = -SGNA;
cannam@95: }
cannam@95: 
cannam@95: static void inv(const N a, N x)
cannam@95: {
cannam@95:      N w, z, one, two;
cannam@95: 
cannam@95:      fromreal(1.0 / toreal(a), x); /* initial guess */
cannam@95:      fromshort(1, one);
cannam@95:      fromshort(2, two);
cannam@95: 
cannam@95:      for (;;) {
cannam@95: 	  /* Newton */
cannam@95: 	  mul(a, x, w);
cannam@95: 	  sub(two, w, z);
cannam@95: 	  if (eq(one, z)) break;
cannam@95: 	  mul(x, z, x);
cannam@95:      }
cannam@95: }
cannam@95: 
cannam@95: 
cannam@95: /* 2 pi */
cannam@95: static const N n2pi = {{
cannam@95:      1, 1,
cannam@95:      {18450, 59017, 1760, 5212, 9779, 4518, 2886, 54545, 18558, 6}
cannam@95: }};
cannam@95: 
cannam@95: /* 1 / 31! */
cannam@95: static const N i31fac = {{ 
cannam@95:      1, -7, 
cannam@95:      {28087, 45433, 51357, 24545, 14291, 3954, 57879, 8109, 38716, 41382}
cannam@95: }};
cannam@95: 
cannam@95: 
cannam@95: /* 1 / 32! */
cannam@95: static const N i32fac = {{
cannam@95:      1, -7,
cannam@95:      {52078, 60811, 3652, 39679, 37310, 47227, 28432, 57597, 13497, 1293}
cannam@95: }};
cannam@95: 
cannam@95: static void msin(const N a, N b)
cannam@95: {
cannam@95:      N a2, g, k;
cannam@95:      int i;
cannam@95: 
cannam@95:      cpy(i31fac, g);
cannam@95:      cpy(g, b);
cannam@95:      mul(a, a, a2);
cannam@95: 
cannam@95:      /* Taylor */
cannam@95:      for (i = 31; i > 1; i -= 2) {
cannam@95: 	  fromshort(i * (i - 1), k);
cannam@95: 	  mul(k, g, g);
cannam@95: 	  mul(a2, b, k);
cannam@95: 	  sub(g, k, b);
cannam@95:      }
cannam@95:      mul(a, b, b);
cannam@95: }
cannam@95: 
cannam@95: static void mcos(const N a, N b)
cannam@95: {
cannam@95:      N a2, g, k;
cannam@95:      int i;
cannam@95: 
cannam@95:      cpy(i32fac, g);
cannam@95:      cpy(g, b);
cannam@95:      mul(a, a, a2);
cannam@95: 
cannam@95:      /* Taylor */
cannam@95:      for (i = 32; i > 0; i -= 2) {
cannam@95: 	  fromshort(i * (i - 1), k);
cannam@95: 	  mul(k, g, g);
cannam@95: 	  mul(a2, b, k);
cannam@95: 	  sub(g, k, b);
cannam@95:      }
cannam@95: }
cannam@95: 
cannam@95: static void by2pi(REAL m, REAL n, N a)
cannam@95: {
cannam@95:      N b;
cannam@95: 
cannam@95:      fromreal(n, b);
cannam@95:      inv(b, a);
cannam@95:      fromreal(m, b);
cannam@95:      mul(a, b, a);
cannam@95:      mul(n2pi, a, a);
cannam@95: }
cannam@95: 
cannam@95: static void sin2pi(REAL m, REAL n, N a);
cannam@95: static void cos2pi(REAL m, REAL n, N a)
cannam@95: {
cannam@95:      N b;
cannam@95:      if (m < 0) cos2pi(-m, n, a);
cannam@95:      else if (m > n * 0.5) cos2pi(n - m, n, a);
cannam@95:      else if (m > n * 0.25) {sin2pi(m - n * 0.25, n, a); neg(a);}
cannam@95:      else if (m > n * 0.125) sin2pi(n * 0.25 - m, n, a);
cannam@95:      else { by2pi(m, n, b); mcos(b, a); }
cannam@95: }
cannam@95: 
cannam@95: static void sin2pi(REAL m, REAL n, N a)
cannam@95: {
cannam@95:      N b;
cannam@95:      if (m < 0)  {sin2pi(-m, n, a); neg(a);}
cannam@95:      else if (m > n * 0.5) {sin2pi(n - m, n, a); neg(a);}
cannam@95:      else if (m > n * 0.25) {cos2pi(m - n * 0.25, n, a);}
cannam@95:      else if (m > n * 0.125) {cos2pi(n * 0.25 - m, n, a);}
cannam@95:      else {by2pi(m, n, b); msin(b, a);}
cannam@95: }
cannam@95: 
cannam@95: /*----------------------------------------------------------------------*/
cannam@95: /* FFT stuff */
cannam@95: 
cannam@95: /* (r0 + i i0)(r1 + i i1) */
cannam@95: static void cmul(N r0, N i0, N r1, N i1, N r2, N i2)
cannam@95: {
cannam@95:      N s, t, q;
cannam@95:      mul(r0, r1, s);
cannam@95:      mul(i0, i1, t);
cannam@95:      sub(s, t, q);
cannam@95:      mul(r0, i1, s);
cannam@95:      mul(i0, r1, t);
cannam@95:      add(s, t, i2);
cannam@95:      cpy(q, r2);
cannam@95: }
cannam@95: 
cannam@95: /* (r0 - i i0)(r1 + i i1) */
cannam@95: static void cmulj(N r0, N i0, N r1, N i1, N r2, N i2)
cannam@95: {
cannam@95:      N s, t, q;
cannam@95:      mul(r0, r1, s);
cannam@95:      mul(i0, i1, t);
cannam@95:      add(s, t, q);
cannam@95:      mul(r0, i1, s);
cannam@95:      mul(i0, r1, t);
cannam@95:      sub(s, t, i2);
cannam@95:      cpy(q, r2);
cannam@95: }
cannam@95: 
cannam@95: static void mcexp(int m, int n, N r, N i)
cannam@95: {
cannam@95:      static int cached_n = -1;
cannam@95:      static N w[64][2];
cannam@95:      int k, j;
cannam@95:      if (n != cached_n) {
cannam@95: 	  for (j = 1, k = 0; j < n; j += j, ++k) {
cannam@95: 	       cos2pi(j, n, w[k][0]);
cannam@95: 	       sin2pi(j, n, w[k][1]);
cannam@95: 	  }
cannam@95: 	  cached_n = n;
cannam@95:      }
cannam@95: 
cannam@95:      fromshort(1, r);
cannam@95:      fromshort(0, i);
cannam@95:      if (m > 0) {
cannam@95: 	  for (k = 0; m; ++k, m >>= 1) 
cannam@95: 	       if (m & 1)
cannam@95: 		    cmul(w[k][0], w[k][1], r, i, r, i);
cannam@95:      } else {
cannam@95: 	  m = -m;
cannam@95: 	  for (k = 0; m; ++k, m >>= 1) 
cannam@95: 	       if (m & 1)
cannam@95: 		    cmulj(w[k][0], w[k][1], r, i, r, i);
cannam@95:      }
cannam@95: }
cannam@95: 
cannam@95: static void bitrev(int n, N *a)
cannam@95: {
cannam@95:      int i, j, m;
cannam@95:      for (i = j = 0; i < n - 1; ++i) {
cannam@95: 	  if (i < j) {
cannam@95: 	       N t;
cannam@95: 	       cpy(a[2*i], t); cpy(a[2*j], a[2*i]); cpy(t, a[2*j]);
cannam@95: 	       cpy(a[2*i+1], t); cpy(a[2*j+1], a[2*i+1]); cpy(t, a[2*j+1]);
cannam@95: 	  }
cannam@95: 
cannam@95: 	  /* bit reversed counter */
cannam@95: 	  m = n; do { m >>= 1; j ^= m; } while (!(j & m));
cannam@95:      }
cannam@95: }
cannam@95: 
cannam@95: static void fft0(int n, N *a, int sign)
cannam@95: {
cannam@95:      int i, j, k;
cannam@95: 
cannam@95:      bitrev(n, a);
cannam@95:      for (i = 1; i < n; i = 2 * i) {
cannam@95: 	  for (j = 0; j < i; ++j) {
cannam@95: 	       N wr, wi;
cannam@95: 	       mcexp(sign * (int)j, 2 * i, wr, wi);
cannam@95: 	       for (k = j; k < n; k += 2 * i) {
cannam@95: 		    N *a0 = a + 2 * k;
cannam@95: 		    N *a1 = a0 + 2 * i;
cannam@95: 		    N r0, i0, r1, i1, t0, t1, xr, xi;
cannam@95: 		    cpy(a0[0], r0); cpy(a0[1], i0);
cannam@95: 		    cpy(a1[0], r1); cpy(a1[1], i1);
cannam@95: 		    mul(r1, wr, t0); mul(i1, wi, t1); sub(t0, t1, xr);
cannam@95: 		    mul(r1, wi, t0); mul(i1, wr, t1); add(t0, t1, xi);
cannam@95: 		    add(r0, xr, a0[0]);  add(i0, xi, a0[1]);
cannam@95: 		    sub(r0, xr, a1[0]);  sub(i0, xi, a1[1]);
cannam@95: 	       }
cannam@95: 	  }
cannam@95:      }
cannam@95: }
cannam@95: 
cannam@95: /* a[2*k]+i*a[2*k+1] = exp(2*pi*i*k^2/(2*n)) */
cannam@95: static void bluestein_sequence(int n, N *a)
cannam@95: {
cannam@95:      int k, ksq, n2 = 2 * n;
cannam@95: 
cannam@95:      ksq = 1; /* (-1)^2 */
cannam@95:      for (k = 0; k < n; ++k) {
cannam@95: 	  /* careful with overflow */
cannam@95: 	  ksq = ksq + 2*k - 1; while (ksq > n2) ksq -= n2;
cannam@95: 	  mcexp(ksq, n2, a[2*k], a[2*k+1]);
cannam@95:      }
cannam@95: }
cannam@95: 
cannam@95: static int pow2_atleast(int x)
cannam@95: {
cannam@95:      int h;
cannam@95:      for (h = 1; h < x; h = 2 * h)
cannam@95: 	  ;
cannam@95:      return h;
cannam@95: }
cannam@95: 
cannam@95: static N *cached_bluestein_w = 0;
cannam@95: static N *cached_bluestein_y = 0;
cannam@95: static int cached_bluestein_n = -1;
cannam@95: 
cannam@95: static void bluestein(int n, N *a)
cannam@95: {
cannam@95:      int nb = pow2_atleast(2 * n);
cannam@95:      N *b = (N *)bench_malloc(2 * nb * sizeof(N));
cannam@95:      N *w = cached_bluestein_w;
cannam@95:      N *y = cached_bluestein_y;
cannam@95:      N nbinv;
cannam@95:      int i;
cannam@95: 
cannam@95:      fromreal(1.0 / nb, nbinv); /* exact because nb = 2^k */
cannam@95: 
cannam@95:      if (cached_bluestein_n != n) {
cannam@95: 	  if (w) bench_free(w);
cannam@95: 	  if (y) bench_free(y);
cannam@95: 	  w = (N *)bench_malloc(2 * n * sizeof(N));
cannam@95: 	  y = (N *)bench_malloc(2 * nb * sizeof(N));
cannam@95: 	  cached_bluestein_n = n;
cannam@95: 	  cached_bluestein_w = w;
cannam@95: 	  cached_bluestein_y = y;
cannam@95: 
cannam@95: 	  bluestein_sequence(n, w);
cannam@95: 	  for (i = 0; i < 2*nb; ++i)  cpy(zero, y[i]);
cannam@95: 
cannam@95: 	  for (i = 0; i < n; ++i) {
cannam@95: 	       cpy(w[2*i], y[2*i]);
cannam@95: 	       cpy(w[2*i+1], y[2*i+1]);
cannam@95: 	  }
cannam@95: 	  for (i = 1; i < n; ++i) {
cannam@95: 	       cpy(w[2*i], y[2*(nb-i)]);
cannam@95: 	       cpy(w[2*i+1], y[2*(nb-i)+1]);
cannam@95: 	  }
cannam@95: 
cannam@95: 	  fft0(nb, y, -1);
cannam@95:      }
cannam@95: 
cannam@95:      for (i = 0; i < 2*nb; ++i)  cpy(zero, b[i]);
cannam@95:      
cannam@95:      for (i = 0; i < n; ++i) 
cannam@95: 	  cmulj(w[2*i], w[2*i+1], a[2*i], a[2*i+1], b[2*i], b[2*i+1]);
cannam@95: 
cannam@95:      /* scaled convolution b * y */
cannam@95:      fft0(nb, b, -1);
cannam@95: 
cannam@95:      for (i = 0; i < nb; ++i) 
cannam@95: 	  cmul(b[2*i], b[2*i+1], y[2*i], y[2*i+1], b[2*i], b[2*i+1]);
cannam@95:      fft0(nb, b, 1);
cannam@95: 
cannam@95:      for (i = 0; i < n; ++i) {
cannam@95: 	  cmulj(w[2*i], w[2*i+1], b[2*i], b[2*i+1], a[2*i], a[2*i+1]);
cannam@95: 	  mul(nbinv, a[2*i], a[2*i]);
cannam@95: 	  mul(nbinv, a[2*i+1], a[2*i+1]);
cannam@95:      }
cannam@95: 
cannam@95:      bench_free(b);
cannam@95: }
cannam@95: 
cannam@95: static void swapri(int n, N *a)
cannam@95: {
cannam@95:      int i;
cannam@95:      for (i = 0; i < n; ++i) {
cannam@95: 	  N t;
cannam@95: 	  cpy(a[2 * i], t);
cannam@95: 	  cpy(a[2 * i + 1], a[2 * i]);
cannam@95: 	  cpy(t, a[2 * i + 1]);
cannam@95:      }
cannam@95: }
cannam@95: 
cannam@95: static void fft1(int n, N *a, int sign)
cannam@95: {
cannam@95:      if (power_of_two(n)) {
cannam@95: 	  fft0(n, a, sign);
cannam@95:      } else {
cannam@95: 	  if (sign == 1) swapri(n, a);
cannam@95: 	  bluestein(n, a);
cannam@95: 	  if (sign == 1) swapri(n, a);
cannam@95:      }
cannam@95: }
cannam@95: 
cannam@95: static void fromrealv(int n, bench_complex *a, N *b)
cannam@95: {
cannam@95:      int i;
cannam@95: 
cannam@95:      for (i = 0; i < n; ++i) {
cannam@95: 	  fromreal(c_re(a[i]), b[2 * i]);
cannam@95: 	  fromreal(c_im(a[i]), b[2 * i + 1]);
cannam@95:      }
cannam@95: }
cannam@95: 
cannam@95: static void compare(int n, N *a, N *b, double *err)
cannam@95: {
cannam@95:      int i;
cannam@95:      double e1, e2, einf;
cannam@95:      double n1, n2, ninf;
cannam@95: 
cannam@95:      e1 = e2 = einf = 0.0;
cannam@95:      n1 = n2 = ninf = 0.0;
cannam@95: 
cannam@95: #    define DO(x1, x2, xinf, var) { 			\
cannam@95:      double d = var;					\
cannam@95:      if (d < 0) d = -d;					\
cannam@95:      x1 += d; x2 += d * d; if (d > xinf) xinf = d;	\
cannam@95: }
cannam@95: 	  
cannam@95:      for (i = 0; i < 2 * n; ++i) {
cannam@95: 	  N dd;
cannam@95: 	  sub(a[i], b[i], dd);
cannam@95: 	  DO(n1, n2, ninf, toreal(a[i]));
cannam@95: 	  DO(e1, e2, einf, toreal(dd));
cannam@95:      }
cannam@95: 
cannam@95: #    undef DO
cannam@95:      err[0] = e1 / n1;
cannam@95:      err[1] = sqrt(e2 / n2);
cannam@95:      err[2] = einf / ninf;
cannam@95: }
cannam@95: 
cannam@95: void fftaccuracy(int n, bench_complex *a, bench_complex *ffta,
cannam@95: 		 int sign, double err[6])
cannam@95: {
cannam@95:      N *b = (N *)bench_malloc(2 * n * sizeof(N));
cannam@95:      N *fftb = (N *)bench_malloc(2 * n * sizeof(N));
cannam@95:      N mn, ninv;
cannam@95:      int i;
cannam@95: 
cannam@95:      fromreal(n, mn); inv(mn, ninv);
cannam@95: 
cannam@95:      /* forward error */
cannam@95:      fromrealv(n, a, b); fromrealv(n, ffta, fftb);
cannam@95:      fft1(n, b, sign);
cannam@95:      compare(n, b, fftb, err);
cannam@95: 
cannam@95:      /* backward error */
cannam@95:      fromrealv(n, a, b); fromrealv(n, ffta, fftb);
cannam@95:      for (i = 0; i < 2 * n; ++i) mul(fftb[i], ninv, fftb[i]);
cannam@95:      fft1(n, fftb, -sign);
cannam@95:      compare(n, b, fftb, err + 3);
cannam@95: 
cannam@95:      bench_free(fftb);
cannam@95:      bench_free(b);
cannam@95: }
cannam@95: 
cannam@95: void fftaccuracy_done(void)
cannam@95: {
cannam@95:      if (cached_bluestein_w) bench_free(cached_bluestein_w);
cannam@95:      if (cached_bluestein_y) bench_free(cached_bluestein_y);
cannam@95:      cached_bluestein_w = 0;
cannam@95:      cached_bluestein_y = 0;
cannam@95:      cached_bluestein_n = -1;
cannam@95: }