1 /*
2 * Copyright (c) 2003, 2007-14 Matteo Frigo
3 * Copyright (c) 2003, 2007-14 Massachusetts Institute of Technology
4 *
5 * This program is free software; you can redistribute it and/or modify
6 * it under the terms of the GNU General Public License as published by
7 * the Free Software Foundation; either version 2 of the License, or
8 * (at your option) any later version.
9 *
10 * This program is distributed in the hope that it will be useful,
11 * but WITHOUT ANY WARRANTY; without even the implied warranty of
12 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13 * GNU General Public License for more details.
14 *
15 * You should have received a copy of the GNU General Public License
16 * along with this program; if not, write to the Free Software
17 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
18 *
19 */
20
21
22 /* Twiddle manipulation */
23
24 #include "kernel/ifftw.h"
25 #include <math.h>
26
27 #define HASHSZ 109
28
29 /* hash table of known twiddle factors */
30 static twid *twlist[HASHSZ];
31
hash(INT n,INT r)32 static INT hash(INT n, INT r)
33 {
34 INT h = n * 17 + r;
35
36 if (h < 0) h = -h;
37
38 return (h % HASHSZ);
39 }
40
equal_instr(const tw_instr * p,const tw_instr * q)41 static int equal_instr(const tw_instr *p, const tw_instr *q)
42 {
43 if (p == q)
44 return 1;
45
46 for (;; ++p, ++q) {
47 if (p->op != q->op)
48 return 0;
49
50 switch (p->op) {
51 case TW_NEXT:
52 return (p->v == q->v); /* p->i is ignored */
53
54 case TW_FULL:
55 case TW_HALF:
56 if (p->v != q->v) return 0; /* p->i is ignored */
57 break;
58
59 default:
60 if (p->v != q->v || p->i != q->i) return 0;
61 break;
62 }
63 }
64 A(0 /* can't happen */);
65 }
66
ok_twid(const twid * t,enum wakefulness wakefulness,const tw_instr * q,INT n,INT r,INT m)67 static int ok_twid(const twid *t,
68 enum wakefulness wakefulness,
69 const tw_instr *q, INT n, INT r, INT m)
70 {
71 return (wakefulness == t->wakefulness &&
72 n == t->n &&
73 r == t->r &&
74 m <= t->m &&
75 equal_instr(t->instr, q));
76 }
77
lookup(enum wakefulness wakefulness,const tw_instr * q,INT n,INT r,INT m)78 static twid *lookup(enum wakefulness wakefulness,
79 const tw_instr *q, INT n, INT r, INT m)
80 {
81 twid *p;
82
83 for (p = twlist[hash(n,r)];
84 p && !ok_twid(p, wakefulness, q, n, r, m);
85 p = p->cdr)
86 ;
87 return p;
88 }
89
twlen0(INT r,const tw_instr * p,INT * vl)90 static INT twlen0(INT r, const tw_instr *p, INT *vl)
91 {
92 INT ntwiddle = 0;
93
94 /* compute length of bytecode program */
95 A(r > 0);
96 for ( ; p->op != TW_NEXT; ++p) {
97 switch (p->op) {
98 case TW_FULL:
99 ntwiddle += (r - 1) * 2;
100 break;
101 case TW_HALF:
102 ntwiddle += (r - 1);
103 break;
104 case TW_CEXP:
105 ntwiddle += 2;
106 break;
107 case TW_COS:
108 case TW_SIN:
109 ntwiddle += 1;
110 break;
111 }
112 }
113
114 *vl = (INT)p->v;
115 return ntwiddle;
116 }
117
X(twiddle_length)118 INT X(twiddle_length)(INT r, const tw_instr *p)
119 {
120 INT vl;
121 return twlen0(r, p, &vl);
122 }
123
compute(enum wakefulness wakefulness,const tw_instr * instr,INT n,INT r,INT m)124 static R *compute(enum wakefulness wakefulness,
125 const tw_instr *instr, INT n, INT r, INT m)
126 {
127 INT ntwiddle, j, vl;
128 R *W, *W0;
129 const tw_instr *p;
130 triggen *t = X(mktriggen)(wakefulness, n);
131
132 p = instr;
133 ntwiddle = twlen0(r, p, &vl);
134
135 A(m % vl == 0);
136
137 W0 = W = (R *)MALLOC((ntwiddle * (m / vl)) * sizeof(R), TWIDDLES);
138
139 for (j = 0; j < m; j += vl) {
140 for (p = instr; p->op != TW_NEXT; ++p) {
141 switch (p->op) {
142 case TW_FULL: {
143 INT i;
144 for (i = 1; i < r; ++i) {
145 A((j + (INT)p->v) * i < n);
146 A((j + (INT)p->v) * i > -n);
147 t->cexp(t, (j + (INT)p->v) * i, W);
148 W += 2;
149 }
150 break;
151 }
152
153 case TW_HALF: {
154 INT i;
155 A((r % 2) == 1);
156 for (i = 1; i + i < r; ++i) {
157 t->cexp(t, MULMOD(i, (j + (INT)p->v), n), W);
158 W += 2;
159 }
160 break;
161 }
162
163 case TW_COS: {
164 R d[2];
165
166 A((j + (INT)p->v) * p->i < n);
167 A((j + (INT)p->v) * p->i > -n);
168 t->cexp(t, (j + (INT)p->v) * (INT)p->i, d);
169 *W++ = d[0];
170 break;
171 }
172
173 case TW_SIN: {
174 R d[2];
175
176 A((j + (INT)p->v) * p->i < n);
177 A((j + (INT)p->v) * p->i > -n);
178 t->cexp(t, (j + (INT)p->v) * (INT)p->i, d);
179 *W++ = d[1];
180 break;
181 }
182
183 case TW_CEXP:
184 A((j + (INT)p->v) * p->i < n);
185 A((j + (INT)p->v) * p->i > -n);
186 t->cexp(t, (j + (INT)p->v) * (INT)p->i, W);
187 W += 2;
188 break;
189 }
190 }
191 }
192
193 X(triggen_destroy)(t);
194 return W0;
195 }
196
mktwiddle(enum wakefulness wakefulness,twid ** pp,const tw_instr * instr,INT n,INT r,INT m)197 static void mktwiddle(enum wakefulness wakefulness,
198 twid **pp, const tw_instr *instr, INT n, INT r, INT m)
199 {
200 twid *p;
201 INT h;
202
203 if ((p = lookup(wakefulness, instr, n, r, m))) {
204 ++p->refcnt;
205 } else {
206 p = (twid *) MALLOC(sizeof(twid), TWIDDLES);
207 p->n = n;
208 p->r = r;
209 p->m = m;
210 p->instr = instr;
211 p->refcnt = 1;
212 p->wakefulness = wakefulness;
213 p->W = compute(wakefulness, instr, n, r, m);
214
215 /* cons! onto twlist */
216 h = hash(n, r);
217 p->cdr = twlist[h];
218 twlist[h] = p;
219 }
220
221 *pp = p;
222 }
223
twiddle_destroy(twid ** pp)224 static void twiddle_destroy(twid **pp)
225 {
226 twid *p = *pp;
227 twid **q;
228
229 if ((--p->refcnt) == 0) {
230 /* remove p from twiddle list */
231 for (q = &twlist[hash(p->n, p->r)]; *q; q = &((*q)->cdr)) {
232 if (*q == p) {
233 *q = p->cdr;
234 X(ifree)(p->W);
235 X(ifree)(p);
236 *pp = 0;
237 return;
238 }
239 }
240 A(0 /* can't happen */ );
241 }
242 }
243
244
X(twiddle_awake)245 void X(twiddle_awake)(enum wakefulness wakefulness, twid **pp,
246 const tw_instr *instr, INT n, INT r, INT m)
247 {
248 switch (wakefulness) {
249 case SLEEPY:
250 twiddle_destroy(pp);
251 break;
252 default:
253 mktwiddle(wakefulness, pp, instr, n, r, m);
254 break;
255 }
256 }
257