Chris@10: /*
Chris@10:  * Copyright (c) 2003, 2007-11 Matteo Frigo
Chris@10:  * Copyright (c) 2003, 2007-11 Massachusetts Institute of Technology
Chris@10:  *
Chris@10:  * This program is free software; you can redistribute it and/or modify
Chris@10:  * it under the terms of the GNU General Public License as published by
Chris@10:  * the Free Software Foundation; either version 2 of the License, or
Chris@10:  * (at your option) any later version.
Chris@10:  *
Chris@10:  * This program is distributed in the hope that it will be useful,
Chris@10:  * but WITHOUT ANY WARRANTY; without even the implied warranty of
Chris@10:  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
Chris@10:  * GNU General Public License for more details.
Chris@10:  *
Chris@10:  * You should have received a copy of the GNU General Public License
Chris@10:  * along with this program; if not, write to the Free Software
Chris@10:  * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA
Chris@10:  *
Chris@10:  */
Chris@10: 
Chris@10: 
Chris@10: /* trigonometric functions */
Chris@10: #include "ifftw.h"
Chris@10: #include <math.h>
Chris@10: 
Chris@10: #if defined(TRIGREAL_IS_LONG_DOUBLE)
Chris@10: #  define COS cosl
Chris@10: #  define SIN sinl
Chris@10: #  define KTRIG(x) (x##L)
Chris@10: #  ifndef HAVE_DECL_SINL
Chris@10:      extern long double sinl(long double x);
Chris@10: #  endif
Chris@10: #  ifndef HAVE_DECL_COSL
Chris@10:      extern long double cosl(long double x);
Chris@10: #  endif
Chris@10: #elif defined(TRIGREAL_IS_QUAD)
Chris@10: #  define COS cosq
Chris@10: #  define SIN sinq
Chris@10: #  define KTRIG(x) (x##Q)
Chris@10:    extern __float128 sinq(__float128 x);
Chris@10:    extern __float128 cosq(__float128 x);
Chris@10: #else
Chris@10: #  define COS cos
Chris@10: #  define SIN sin
Chris@10: #  define KTRIG(x) (x)
Chris@10: #endif
Chris@10: 
Chris@10: static const trigreal K2PI =
Chris@10:     KTRIG(6.2831853071795864769252867665590057683943388);
Chris@10: #define by2pi(m, n) ((K2PI * (m)) / (n))
Chris@10: 
Chris@10: /*
Chris@10:  * Improve accuracy by reducing x to range [0..1/8]
Chris@10:  * before multiplication by 2 * PI.
Chris@10:  */
Chris@10: 
Chris@10: static void real_cexp(INT m, INT n, trigreal *out)
Chris@10: {
Chris@10:      trigreal theta, c, s, t;
Chris@10:      unsigned octant = 0;
Chris@10:      INT quarter_n = n;
Chris@10: 
Chris@10:      n += n; n += n;
Chris@10:      m += m; m += m;
Chris@10: 
Chris@10:      if (m < 0) m += n;
Chris@10:      if (m > n - m) { m = n - m; octant |= 4; }
Chris@10:      if (m - quarter_n > 0) { m = m - quarter_n; octant |= 2; }
Chris@10:      if (m > quarter_n - m) { m = quarter_n - m; octant |= 1; }
Chris@10: 
Chris@10:      theta = by2pi(m, n);
Chris@10:      c = COS(theta); s = SIN(theta);
Chris@10: 
Chris@10:      if (octant & 1) { t = c; c = s; s = t; }
Chris@10:      if (octant & 2) { t = c; c = -s; s = t; }
Chris@10:      if (octant & 4) { s = -s; }
Chris@10: 
Chris@10:      out[0] = c; 
Chris@10:      out[1] = s; 
Chris@10: }
Chris@10: 
Chris@10: static INT choose_twshft(INT n)
Chris@10: {
Chris@10:      INT log2r = 0;
Chris@10:      while (n > 0) {
Chris@10: 	  ++log2r;
Chris@10: 	  n /= 4;
Chris@10:      }
Chris@10:      return log2r;
Chris@10: }
Chris@10: 
Chris@10: static void cexpl_sqrtn_table(triggen *p, INT m, trigreal *res)
Chris@10: {
Chris@10:      m += p->n * (m < 0);
Chris@10: 
Chris@10:      {
Chris@10: 	  INT m0 = m & p->twmsk;
Chris@10: 	  INT m1 = m >> p->twshft;
Chris@10: 	  trigreal wr0 = p->W0[2 * m0];
Chris@10: 	  trigreal wi0 = p->W0[2 * m0 + 1];
Chris@10: 	  trigreal wr1 = p->W1[2 * m1];
Chris@10: 	  trigreal wi1 = p->W1[2 * m1 + 1];
Chris@10: 
Chris@10: 	  res[0] = wr1 * wr0 - wi1 * wi0;
Chris@10: 	  res[1] = wi1 * wr0 + wr1 * wi0;
Chris@10:      }
Chris@10: }
Chris@10: 
Chris@10: /* multiply (xr, xi) by exp(FFT_SIGN * 2*pi*i*m/n) */
Chris@10: static void rotate_sqrtn_table(triggen *p, INT m, R xr, R xi, R *res)
Chris@10: {
Chris@10:      m += p->n * (m < 0);
Chris@10: 
Chris@10:      {
Chris@10: 	  INT m0 = m & p->twmsk;
Chris@10: 	  INT m1 = m >> p->twshft;
Chris@10: 	  trigreal wr0 = p->W0[2 * m0];
Chris@10: 	  trigreal wi0 = p->W0[2 * m0 + 1];
Chris@10: 	  trigreal wr1 = p->W1[2 * m1];
Chris@10: 	  trigreal wi1 = p->W1[2 * m1 + 1];
Chris@10: 	  trigreal wr = wr1 * wr0 - wi1 * wi0;
Chris@10: 	  trigreal wi = wi1 * wr0 + wr1 * wi0;
Chris@10: 
Chris@10: #if FFT_SIGN == -1
Chris@10: 	  res[0] = xr * wr + xi * wi;
Chris@10: 	  res[1] = xi * wr - xr * wi;
Chris@10: #else
Chris@10: 	  res[0] = xr * wr - xi * wi;
Chris@10: 	  res[1] = xi * wr + xr * wi;
Chris@10: #endif
Chris@10:      }
Chris@10: }
Chris@10: 
Chris@10: static void cexpl_sincos(triggen *p, INT m, trigreal *res)
Chris@10: {
Chris@10:      real_cexp(m, p->n, res);
Chris@10: }
Chris@10: 
Chris@10: static void cexp_zero(triggen *p, INT m, R *res)
Chris@10: {
Chris@10:      UNUSED(p); UNUSED(m);
Chris@10:      res[0] = 0;
Chris@10:      res[1] = 0;
Chris@10: }
Chris@10: 
Chris@10: static void cexpl_zero(triggen *p, INT m, trigreal *res)
Chris@10: {
Chris@10:      UNUSED(p); UNUSED(m);
Chris@10:      res[0] = 0;
Chris@10:      res[1] = 0;
Chris@10: }
Chris@10: 
Chris@10: static void cexp_generic(triggen *p, INT m, R *res)
Chris@10: {
Chris@10:      trigreal resl[2];
Chris@10:      p->cexpl(p, m, resl);
Chris@10:      res[0] = (R)resl[0];
Chris@10:      res[1] = (R)resl[1];
Chris@10: }
Chris@10: 
Chris@10: static void rotate_generic(triggen *p, INT m, R xr, R xi, R *res)
Chris@10: {
Chris@10:      trigreal w[2];
Chris@10:      p->cexpl(p, m, w);
Chris@10:      res[0] = xr * w[0] - xi * (FFT_SIGN * w[1]);
Chris@10:      res[1] = xi * w[0] + xr * (FFT_SIGN * w[1]);
Chris@10: }
Chris@10: 
Chris@10: triggen *X(mktriggen)(enum wakefulness wakefulness, INT n)
Chris@10: {
Chris@10:      INT i, n0, n1;
Chris@10:      triggen *p = (triggen *)MALLOC(sizeof(*p), TWIDDLES);
Chris@10: 
Chris@10:      p->n = n;
Chris@10:      p->W0 = p->W1 = 0;
Chris@10:      p->cexp = 0;
Chris@10:      p->rotate = 0;
Chris@10: 
Chris@10:      switch (wakefulness) {
Chris@10: 	 case SLEEPY:
Chris@10: 	      A(0 /* can't happen */);
Chris@10: 	      break;
Chris@10: 
Chris@10: 	 case AWAKE_SQRTN_TABLE: {
Chris@10: 	      INT twshft = choose_twshft(n);
Chris@10: 
Chris@10: 	      p->twshft = twshft;
Chris@10: 	      p->twradix = ((INT)1) << twshft;
Chris@10: 	      p->twmsk = p->twradix - 1;
Chris@10: 
Chris@10: 	      n0 = p->twradix;
Chris@10: 	      n1 = (n + n0 - 1) / n0;
Chris@10: 
Chris@10: 	      p->W0 = (trigreal *)MALLOC(n0 * 2 * sizeof(trigreal), TWIDDLES);
Chris@10: 	      p->W1 = (trigreal *)MALLOC(n1 * 2 * sizeof(trigreal), TWIDDLES);
Chris@10: 
Chris@10: 	      for (i = 0; i < n0; ++i) 
Chris@10: 		   real_cexp(i, n, p->W0 + 2 * i);
Chris@10: 
Chris@10: 	      for (i = 0; i < n1; ++i) 
Chris@10: 		   real_cexp(i * p->twradix, n, p->W1 + 2 * i);
Chris@10: 
Chris@10: 	      p->cexpl = cexpl_sqrtn_table;
Chris@10: 	      p->rotate = rotate_sqrtn_table;
Chris@10: 	      break;
Chris@10: 	 }
Chris@10: 
Chris@10: 	 case AWAKE_SINCOS: 
Chris@10: 	      p->cexpl = cexpl_sincos;
Chris@10: 	      break;
Chris@10: 
Chris@10: 	 case AWAKE_ZERO: 
Chris@10: 	      p->cexp = cexp_zero;
Chris@10: 	      p->cexpl = cexpl_zero;
Chris@10: 	      break;
Chris@10:      }
Chris@10: 
Chris@10:      if (!p->cexp) {
Chris@10: 	  if (sizeof(trigreal) == sizeof(R))
Chris@10: 	       p->cexp = (void (*)(triggen *, INT, R *))p->cexpl;
Chris@10: 	  else
Chris@10: 	       p->cexp = cexp_generic;
Chris@10:      }
Chris@10:      if (!p->rotate)     
Chris@10: 	       p->rotate = rotate_generic;
Chris@10:      return p;
Chris@10: }
Chris@10: 
Chris@10: void X(triggen_destroy)(triggen *p)
Chris@10: {
Chris@10:      X(ifree0)(p->W0);
Chris@10:      X(ifree0)(p->W1);
Chris@10:      X(ifree)(p);
Chris@10: }