cannam@95: /*
cannam@95:  * Copyright (c) 2003, 2007-11 Matteo Frigo
cannam@95:  * Copyright (c) 2003, 2007-11 Massachusetts Institute of Technology
cannam@95:  *
cannam@95:  * This program is free software; you can redistribute it and/or modify
cannam@95:  * it under the terms of the GNU General Public License as published by
cannam@95:  * the Free Software Foundation; either version 2 of the License, or
cannam@95:  * (at your option) any later version.
cannam@95:  *
cannam@95:  * This program is distributed in the hope that it will be useful,
cannam@95:  * but WITHOUT ANY WARRANTY; without even the implied warranty of
cannam@95:  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
cannam@95:  * GNU General Public License for more details.
cannam@95:  *
cannam@95:  * You should have received a copy of the GNU General Public License
cannam@95:  * along with this program; if not, write to the Free Software
cannam@95:  * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA
cannam@95:  *
cannam@95:  */
cannam@95: 
cannam@95: 
cannam@95: /* Twiddle manipulation */
cannam@95: 
cannam@95: #include "ifftw.h"
cannam@95: #include <math.h>
cannam@95: 
cannam@95: #define HASHSZ 109
cannam@95: 
cannam@95: /* hash table of known twiddle factors */
cannam@95: static twid *twlist[HASHSZ];
cannam@95: 
cannam@95: static INT hash(INT n, INT r)
cannam@95: {
cannam@95:      INT h = n * 17 + r;
cannam@95: 
cannam@95:      if (h < 0) h = -h;
cannam@95: 
cannam@95:      return (h % HASHSZ);
cannam@95: }
cannam@95: 
cannam@95: static int equal_instr(const tw_instr *p, const tw_instr *q)
cannam@95: {
cannam@95:      if (p == q)
cannam@95:           return 1;
cannam@95: 
cannam@95:      for (;; ++p, ++q) {
cannam@95:           if (p->op != q->op)
cannam@95: 	       return 0;
cannam@95: 
cannam@95: 	  switch (p->op) {
cannam@95: 	      case TW_NEXT:
cannam@95: 		   return (p->v == q->v); /* p->i is ignored */
cannam@95: 
cannam@95: 	      case TW_FULL:
cannam@95: 	      case TW_HALF:
cannam@95: 		   if (p->v != q->v) return 0; /* p->i is ignored */
cannam@95: 		   break;
cannam@95: 
cannam@95: 	      default:
cannam@95: 		   if (p->v != q->v || p->i != q->i) return 0;
cannam@95: 		   break;
cannam@95: 	  }
cannam@95:      }
cannam@95:      A(0 /* can't happen */);
cannam@95: }
cannam@95: 
cannam@95: static int ok_twid(const twid *t, 
cannam@95: 		   enum wakefulness wakefulness,
cannam@95: 		   const tw_instr *q, INT n, INT r, INT m)
cannam@95: {
cannam@95:      return (wakefulness == t->wakefulness &&
cannam@95: 	     n == t->n &&
cannam@95: 	     r == t->r && 
cannam@95: 	     m <= t->m && 
cannam@95: 	     equal_instr(t->instr, q));
cannam@95: }
cannam@95: 
cannam@95: static twid *lookup(enum wakefulness wakefulness,
cannam@95: 		    const tw_instr *q, INT n, INT r, INT m)
cannam@95: {
cannam@95:      twid *p;
cannam@95: 
cannam@95:      for (p = twlist[hash(n,r)]; 
cannam@95: 	  p && !ok_twid(p, wakefulness, q, n, r, m); 
cannam@95: 	  p = p->cdr)
cannam@95:           ;
cannam@95:      return p;
cannam@95: }
cannam@95: 
cannam@95: static INT twlen0(INT r, const tw_instr *p, INT *vl)
cannam@95: {
cannam@95:      INT ntwiddle = 0;
cannam@95: 
cannam@95:      /* compute length of bytecode program */
cannam@95:      A(r > 0);
cannam@95:      for ( ; p->op != TW_NEXT; ++p) {
cannam@95: 	  switch (p->op) {
cannam@95: 	      case TW_FULL:
cannam@95: 		   ntwiddle += (r - 1) * 2;
cannam@95: 		   break;
cannam@95: 	      case TW_HALF:
cannam@95: 		   ntwiddle += (r - 1);
cannam@95: 		   break;
cannam@95: 	      case TW_CEXP:
cannam@95: 		   ntwiddle += 2;
cannam@95: 		   break;
cannam@95: 	      case TW_COS:
cannam@95: 	      case TW_SIN:
cannam@95: 		   ntwiddle += 1;
cannam@95: 		   break;
cannam@95: 	  }
cannam@95:      }
cannam@95: 
cannam@95:      *vl = (INT)p->v;
cannam@95:      return ntwiddle;
cannam@95: }
cannam@95: 
cannam@95: INT X(twiddle_length)(INT r, const tw_instr *p)
cannam@95: {
cannam@95:      INT vl;
cannam@95:      return twlen0(r, p, &vl);
cannam@95: }
cannam@95: 
cannam@95: static R *compute(enum wakefulness wakefulness,
cannam@95: 		  const tw_instr *instr, INT n, INT r, INT m)
cannam@95: {
cannam@95:      INT ntwiddle, j, vl;
cannam@95:      R *W, *W0;
cannam@95:      const tw_instr *p;
cannam@95:      triggen *t = X(mktriggen)(wakefulness, n);
cannam@95: 
cannam@95:      p = instr;
cannam@95:      ntwiddle = twlen0(r, p, &vl);
cannam@95: 
cannam@95:      A(m % vl == 0);
cannam@95: 
cannam@95:      W0 = W = (R *)MALLOC((ntwiddle * (m / vl)) * sizeof(R), TWIDDLES);
cannam@95: 
cannam@95:      for (j = 0; j < m; j += vl) {
cannam@95:           for (p = instr; p->op != TW_NEXT; ++p) {
cannam@95: 	       switch (p->op) {
cannam@95: 		   case TW_FULL: {
cannam@95: 			INT i;
cannam@95: 			for (i = 1; i < r; ++i) {
cannam@95: 			     A((j + (INT)p->v) * i < n);
cannam@95: 			     A((j + (INT)p->v) * i > -n);
cannam@95: 			     t->cexp(t, (j + (INT)p->v) * i, W);
cannam@95: 			     W += 2;
cannam@95: 			}
cannam@95: 			break;
cannam@95: 		   }
cannam@95: 
cannam@95: 		   case TW_HALF: {
cannam@95: 			INT i;
cannam@95: 			A((r % 2) == 1);
cannam@95: 			for (i = 1; i + i < r; ++i) {
cannam@95: 			     t->cexp(t, MULMOD(i, (j + (INT)p->v), n), W);
cannam@95: 			     W += 2;
cannam@95: 			}
cannam@95: 			break;
cannam@95: 		   }
cannam@95: 
cannam@95: 		   case TW_COS: {
cannam@95: 			R d[2];
cannam@95: 
cannam@95: 			A((j + (INT)p->v) * p->i < n);
cannam@95: 			A((j + (INT)p->v) * p->i > -n);
cannam@95: 			t->cexp(t, (j + (INT)p->v) * (INT)p->i, d);
cannam@95: 			*W++ = d[0];
cannam@95: 			break;
cannam@95: 		   }
cannam@95: 
cannam@95: 		   case TW_SIN: {
cannam@95: 			R d[2];
cannam@95: 
cannam@95: 			A((j + (INT)p->v) * p->i < n);
cannam@95: 			A((j + (INT)p->v) * p->i > -n);
cannam@95: 			t->cexp(t, (j + (INT)p->v) * (INT)p->i, d);
cannam@95: 			*W++ = d[1];
cannam@95: 			break;
cannam@95: 		   }
cannam@95: 
cannam@95: 		   case TW_CEXP:
cannam@95: 			A((j + (INT)p->v) * p->i < n);
cannam@95: 			A((j + (INT)p->v) * p->i > -n);
cannam@95: 			t->cexp(t, (j + (INT)p->v) * (INT)p->i, W);
cannam@95: 			W += 2;
cannam@95: 			break;
cannam@95: 	       }
cannam@95: 	  }
cannam@95:      }
cannam@95: 
cannam@95:      X(triggen_destroy)(t);
cannam@95:      return W0;
cannam@95: }
cannam@95: 
cannam@95: static void mktwiddle(enum wakefulness wakefulness,
cannam@95: 		      twid **pp, const tw_instr *instr, INT n, INT r, INT m)
cannam@95: {
cannam@95:      twid *p;
cannam@95:      INT h;
cannam@95: 
cannam@95:      if ((p = lookup(wakefulness, instr, n, r, m))) {
cannam@95:           ++p->refcnt;
cannam@95:      } else {
cannam@95: 	  p = (twid *) MALLOC(sizeof(twid), TWIDDLES);
cannam@95: 	  p->n = n;
cannam@95: 	  p->r = r;
cannam@95: 	  p->m = m;
cannam@95: 	  p->instr = instr;
cannam@95: 	  p->refcnt = 1;
cannam@95: 	  p->wakefulness = wakefulness;
cannam@95: 	  p->W = compute(wakefulness, instr, n, r, m);
cannam@95: 
cannam@95: 	  /* cons! onto twlist */
cannam@95: 	  h = hash(n, r);
cannam@95: 	  p->cdr = twlist[h];
cannam@95: 	  twlist[h] = p;
cannam@95:      }
cannam@95: 
cannam@95:      *pp = p;
cannam@95: }
cannam@95: 
cannam@95: static void twiddle_destroy(twid **pp)
cannam@95: {
cannam@95:      twid *p = *pp;
cannam@95:      twid **q;
cannam@95: 
cannam@95:      if ((--p->refcnt) == 0) {
cannam@95: 	  /* remove p from twiddle list */
cannam@95: 	  for (q = &twlist[hash(p->n, p->r)]; *q; q = &((*q)->cdr)) {
cannam@95: 	       if (*q == p) {
cannam@95: 		    *q = p->cdr;
cannam@95: 		    X(ifree)(p->W);
cannam@95: 		    X(ifree)(p);
cannam@95: 		    *pp = 0;
cannam@95: 		    return;
cannam@95: 	       }
cannam@95: 	  }
cannam@95: 	  A(0 /* can't happen */ );
cannam@95:      }
cannam@95: }
cannam@95: 
cannam@95: 
cannam@95: void X(twiddle_awake)(enum wakefulness wakefulness, twid **pp, 
cannam@95: 		      const tw_instr *instr, INT n, INT r, INT m)
cannam@95: {
cannam@95:      switch (wakefulness) {
cannam@95: 	 case SLEEPY: 
cannam@95: 	      twiddle_destroy(pp);
cannam@95: 	      break;
cannam@95: 	 default:
cannam@95: 	      mktwiddle(wakefulness, pp, instr, n, r, m);
cannam@95: 	      break;
cannam@95:      }
cannam@95: }