view Lib/fftw-3.2.1/kernel/trig.c @ 1:e86e9c111b29

Updates stuff that potentially fixes the memory leak and also makes it work on Windows and Linux (Need to test). Still have to fix fftw include for linux in Jucer.

/*
 * Copyright (c) 2003, 2007-8 Matteo Frigo
 * Copyright (c) 2003, 2007-8 Massachusetts Institute of Technology
 *
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 2 of the License, or
 * (at your option) any later version.
 *
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
 *
 */


/* trigonometric functions */
#include "ifftw.h"
#include <math.h>

#ifdef TRIGREAL_IS_LONG_DOUBLE
#  define COS cosl
#  define SIN sinl
#  define KTRIG(x) (x##L)
#  ifndef HAVE_DECL_SINL
     extern long double sinl(long double x);
#  endif
#  ifndef HAVE_DECL_COSL
     extern long double cosl(long double x);
#  endif
#else
#  define COS cos
#  define SIN sin
#  define KTRIG(x) (x)
#endif

static const trigreal K2PI =
    KTRIG(6.2831853071795864769252867665590057683943388);
#define by2pi(m, n) ((K2PI * (m)) / (n))

/*
 * Improve accuracy by reducing x to range [0..1/8]
 * before multiplication by 2 * PI.
 */

static void real_cexp(INT m, INT n, trigreal *out)
{
     trigreal theta, c, s, t;
     unsigned octant = 0;
     INT quarter_n = n;

     n += n; n += n;
     m += m; m += m;

     if (m < 0) m += n;
     if (m > n - m) { m = n - m; octant |= 4; }
     if (m - quarter_n > 0) { m = m - quarter_n; octant |= 2; }
     if (m > quarter_n - m) { m = quarter_n - m; octant |= 1; }

     theta = by2pi(m, n);
     c = COS(theta); s = SIN(theta);

     if (octant & 1) { t = c; c = s; s = t; }
     if (octant & 2) { t = c; c = -s; s = t; }
     if (octant & 4) { s = -s; }

     out[0] = c; 
     out[1] = s; 
}

static INT choose_twshft(INT n)
{
     INT log2r = 0;
     while (n > 0) {
	  ++log2r;
	  n /= 4;
     }
     return log2r;
}

static void cexpl_sqrtn_table(triggen *p, INT m, trigreal *res)
{
     m += p->n * (m < 0);

     {
	  INT m0 = m & p->twmsk;
	  INT m1 = m >> p->twshft;
	  trigreal wr0 = p->W0[2 * m0];
	  trigreal wi0 = p->W0[2 * m0 + 1];
	  trigreal wr1 = p->W1[2 * m1];
	  trigreal wi1 = p->W1[2 * m1 + 1];

	  res[0] = wr1 * wr0 - wi1 * wi0;
	  res[1] = wi1 * wr0 + wr1 * wi0;
     }
}

/* multiply (xr, xi) by exp(FFT_SIGN * 2*pi*i*m/n) */
static void rotate_sqrtn_table(triggen *p, INT m, R xr, R xi, R *res)
{
     m += p->n * (m < 0);

     {
	  INT m0 = m & p->twmsk;
	  INT m1 = m >> p->twshft;
	  trigreal wr0 = p->W0[2 * m0];
	  trigreal wi0 = p->W0[2 * m0 + 1];
	  trigreal wr1 = p->W1[2 * m1];
	  trigreal wi1 = p->W1[2 * m1 + 1];
	  trigreal wr = wr1 * wr0 - wi1 * wi0;
	  trigreal wi = wi1 * wr0 + wr1 * wi0;

#if FFT_SIGN == -1
	  res[0] = xr * wr + xi * wi;
	  res[1] = xi * wr - xr * wi;
#else
	  res[0] = xr * wr - xi * wi;
	  res[1] = xi * wr + xr * wi;
#endif
     }
}

static void cexpl_sincos(triggen *p, INT m, trigreal *res)
{
     real_cexp(m, p->n, res);
}

static void cexp_zero(triggen *p, INT m, R *res)
{
     UNUSED(p); UNUSED(m);
     res[0] = 0;
     res[1] = 0;
}

static void cexpl_zero(triggen *p, INT m, trigreal *res)
{
     UNUSED(p); UNUSED(m);
     res[0] = 0;
     res[1] = 0;
}

static void cexp_generic(triggen *p, INT m, R *res)
{
     trigreal resl[2];
     p->cexpl(p, m, resl);
     res[0] = (R)resl[0];
     res[1] = (R)resl[1];
}

static void rotate_generic(triggen *p, INT m, R xr, R xi, R *res)
{
     trigreal w[2];
     p->cexpl(p, m, w);
     res[0] = xr * w[0] - xi * (FFT_SIGN * w[1]);
     res[1] = xi * w[0] + xr * (FFT_SIGN * w[1]);
}

triggen *X(mktriggen)(enum wakefulness wakefulness, INT n)
{
     INT i, n0, n1;
     triggen *p = (triggen *)MALLOC(sizeof(*p), TWIDDLES);

     p->n = n;
     p->W0 = p->W1 = 0;
     p->cexp = 0;
     p->rotate = 0;

     switch (wakefulness) {
	 case SLEEPY:
	      A(0 /* can't happen */);
	      break;

	 case AWAKE_SQRTN_TABLE: {
	      INT twshft = choose_twshft(n);

	      p->twshft = twshft;
	      p->twradix = ((INT)1) << twshft;
	      p->twmsk = p->twradix - 1;

	      n0 = p->twradix;
	      n1 = (n + n0 - 1) / n0;

	      p->W0 = (trigreal *)MALLOC(n0 * 2 * sizeof(trigreal), TWIDDLES);
	      p->W1 = (trigreal *)MALLOC(n1 * 2 * sizeof(trigreal), TWIDDLES);

	      for (i = 0; i < n0; ++i) 
		   real_cexp(i, n, p->W0 + 2 * i);

	      for (i = 0; i < n1; ++i) 
		   real_cexp(i * p->twradix, n, p->W1 + 2 * i);

	      p->cexpl = cexpl_sqrtn_table;
	      p->rotate = rotate_sqrtn_table;
	      break;
	 }

	 case AWAKE_SINCOS: 
	      p->cexpl = cexpl_sincos;
	      break;

	 case AWAKE_ZERO: 
	      p->cexp = cexp_zero;
	      p->cexpl = cexpl_zero;
	      break;
     }

     if (!p->cexp) {
	  if (sizeof(trigreal) == sizeof(R))
	       p->cexp = (void (*)(triggen *, INT, R *))p->cexpl;
	  else
	       p->cexp = cexp_generic;
     }
     if (!p->rotate)     
	       p->rotate = rotate_generic;
     return p;
}

void X(triggen_destroy)(triggen *p)
{
     X(ifree0)(p->W0);
     X(ifree0)(p->W1);
     X(ifree)(p);
}