diff src/fftw-3.3.8/kernel/twiddle.c @ 167:bd3cc4d1df30

Add FFTW 3.3.8 source, and a Linux build

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/src/fftw-3.3.8/kernel/twiddle.c	Tue Nov 19 14:52:55 2019 +0000
@@ -0,0 +1,256 @@
+/*
+ * Copyright (c) 2003, 2007-14 Matteo Frigo
+ * Copyright (c) 2003, 2007-14 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., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA
+ *
+ */
+
+
+/* Twiddle manipulation */
+
+#include "kernel/ifftw.h"
+#include <math.h>
+
+#define HASHSZ 109
+
+/* hash table of known twiddle factors */
+static twid *twlist[HASHSZ];
+
+static INT hash(INT n, INT r)
+{
+     INT h = n * 17 + r;
+
+     if (h < 0) h = -h;
+
+     return (h % HASHSZ);
+}
+
+static int equal_instr(const tw_instr *p, const tw_instr *q)
+{
+     if (p == q)
+          return 1;
+
+     for (;; ++p, ++q) {
+          if (p->op != q->op)
+	       return 0;
+
+	  switch (p->op) {
+	      case TW_NEXT:
+		   return (p->v == q->v); /* p->i is ignored */
+
+	      case TW_FULL:
+	      case TW_HALF:
+		   if (p->v != q->v) return 0; /* p->i is ignored */
+		   break;
+
+	      default:
+		   if (p->v != q->v || p->i != q->i) return 0;
+		   break;
+	  }
+     }
+     A(0 /* can't happen */);
+}
+
+static int ok_twid(const twid *t, 
+		   enum wakefulness wakefulness,
+		   const tw_instr *q, INT n, INT r, INT m)
+{
+     return (wakefulness == t->wakefulness &&
+	     n == t->n &&
+	     r == t->r && 
+	     m <= t->m && 
+	     equal_instr(t->instr, q));
+}
+
+static twid *lookup(enum wakefulness wakefulness,
+		    const tw_instr *q, INT n, INT r, INT m)
+{
+     twid *p;
+
+     for (p = twlist[hash(n,r)]; 
+	  p && !ok_twid(p, wakefulness, q, n, r, m); 
+	  p = p->cdr)
+          ;
+     return p;
+}
+
+static INT twlen0(INT r, const tw_instr *p, INT *vl)
+{
+     INT ntwiddle = 0;
+
+     /* compute length of bytecode program */
+     A(r > 0);
+     for ( ; p->op != TW_NEXT; ++p) {
+	  switch (p->op) {
+	      case TW_FULL:
+		   ntwiddle += (r - 1) * 2;
+		   break;
+	      case TW_HALF:
+		   ntwiddle += (r - 1);
+		   break;
+	      case TW_CEXP:
+		   ntwiddle += 2;
+		   break;
+	      case TW_COS:
+	      case TW_SIN:
+		   ntwiddle += 1;
+		   break;
+	  }
+     }
+
+     *vl = (INT)p->v;
+     return ntwiddle;
+}
+
+INT X(twiddle_length)(INT r, const tw_instr *p)
+{
+     INT vl;
+     return twlen0(r, p, &vl);
+}
+
+static R *compute(enum wakefulness wakefulness,
+		  const tw_instr *instr, INT n, INT r, INT m)
+{
+     INT ntwiddle, j, vl;
+     R *W, *W0;
+     const tw_instr *p;
+     triggen *t = X(mktriggen)(wakefulness, n);
+
+     p = instr;
+     ntwiddle = twlen0(r, p, &vl);
+
+     A(m % vl == 0);
+
+     W0 = W = (R *)MALLOC((ntwiddle * (m / vl)) * sizeof(R), TWIDDLES);
+
+     for (j = 0; j < m; j += vl) {
+          for (p = instr; p->op != TW_NEXT; ++p) {
+	       switch (p->op) {
+		   case TW_FULL: {
+			INT i;
+			for (i = 1; i < r; ++i) {
+			     A((j + (INT)p->v) * i < n);
+			     A((j + (INT)p->v) * i > -n);
+			     t->cexp(t, (j + (INT)p->v) * i, W);
+			     W += 2;
+			}
+			break;
+		   }
+
+		   case TW_HALF: {
+			INT i;
+			A((r % 2) == 1);
+			for (i = 1; i + i < r; ++i) {
+			     t->cexp(t, MULMOD(i, (j + (INT)p->v), n), W);
+			     W += 2;
+			}
+			break;
+		   }
+
+		   case TW_COS: {
+			R d[2];
+
+			A((j + (INT)p->v) * p->i < n);
+			A((j + (INT)p->v) * p->i > -n);
+			t->cexp(t, (j + (INT)p->v) * (INT)p->i, d);
+			*W++ = d[0];
+			break;
+		   }
+
+		   case TW_SIN: {
+			R d[2];
+
+			A((j + (INT)p->v) * p->i < n);
+			A((j + (INT)p->v) * p->i > -n);
+			t->cexp(t, (j + (INT)p->v) * (INT)p->i, d);
+			*W++ = d[1];
+			break;
+		   }
+
+		   case TW_CEXP:
+			A((j + (INT)p->v) * p->i < n);
+			A((j + (INT)p->v) * p->i > -n);
+			t->cexp(t, (j + (INT)p->v) * (INT)p->i, W);
+			W += 2;
+			break;
+	       }
+	  }
+     }
+
+     X(triggen_destroy)(t);
+     return W0;
+}
+
+static void mktwiddle(enum wakefulness wakefulness,
+		      twid **pp, const tw_instr *instr, INT n, INT r, INT m)
+{
+     twid *p;
+     INT h;
+
+     if ((p = lookup(wakefulness, instr, n, r, m))) {
+          ++p->refcnt;
+     } else {
+	  p = (twid *) MALLOC(sizeof(twid), TWIDDLES);
+	  p->n = n;
+	  p->r = r;
+	  p->m = m;
+	  p->instr = instr;
+	  p->refcnt = 1;
+	  p->wakefulness = wakefulness;
+	  p->W = compute(wakefulness, instr, n, r, m);
+
+	  /* cons! onto twlist */
+	  h = hash(n, r);
+	  p->cdr = twlist[h];
+	  twlist[h] = p;
+     }
+
+     *pp = p;
+}
+
+static void twiddle_destroy(twid **pp)
+{
+     twid *p = *pp;
+     twid **q;
+
+     if ((--p->refcnt) == 0) {
+	  /* remove p from twiddle list */
+	  for (q = &twlist[hash(p->n, p->r)]; *q; q = &((*q)->cdr)) {
+	       if (*q == p) {
+		    *q = p->cdr;
+		    X(ifree)(p->W);
+		    X(ifree)(p);
+		    *pp = 0;
+		    return;
+	       }
+	  }
+	  A(0 /* can't happen */ );
+     }
+}
+
+
+void X(twiddle_awake)(enum wakefulness wakefulness, twid **pp, 
+		      const tw_instr *instr, INT n, INT r, INT m)
+{
+     switch (wakefulness) {
+	 case SLEEPY: 
+	      twiddle_destroy(pp);
+	      break;
+	 default:
+	      mktwiddle(wakefulness, pp, instr, n, r, m);
+	      break;
+     }
+}