← Back to branch summary

~diresu/blender/blender-command-port

« back to all changes in this revision

Viewing changes to extern/fftw/kernel/twiddle.c

  • Committer: theeth
  • Date: 2008-10-14 16:52:04 UTC
  • Revision ID: vcs-imports@canonical.com-20081014165204-r32w2gm6s0osvdhn
copy back trunk

Show diffs side-by-side

added added

removed removed

Lines of Context:
extern/fftw/kernel/twiddle.c
 
1
/*
 
2
 * Copyright (c) 2003, 2006 Matteo Frigo
 
3
 * Copyright (c) 2003, 2006 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., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
 
18
 *
 
19
 */
 
20
 
 
21
/* $Id: twiddle.c,v 1.33 2006-01-05 22:39:02 athena Exp $ */
 
22
 
 
23
/* Twiddle manipulation */
 
24
 
 
25
#include "ifftw.h"
 
26
#include <math.h>
 
27
 
 
28
#define HASHSZ 109
 
29
 
 
30
/* hash table of known twiddle factors */
 
31
static twid *twlist[HASHSZ];
 
32
 
 
33
static INT hash(INT n, INT r)
 
34
{
 
35
     INT h = n * 17 + r;
 
36
 
 
37
     if (h < 0) h = -h;
 
38
 
 
39
     return (h % HASHSZ);
 
40
}
 
41
 
 
42
static int equal_instr(const tw_instr *p, const tw_instr *q)
 
43
{
 
44
     if (p == q)
 
45
          return 1;
 
46
 
 
47
     for (;; ++p, ++q) {
 
48
          if (p->op != q->op)
 
49
               return 0;
 
50
 
 
51
          switch (p->op) {
 
52
              case TW_NEXT:
 
53
                   return 1;
 
54
 
 
55
              case TW_FULL:
 
56
              case TW_HALF:
 
57
                   if (p->v != q->v) return 0; /* p->i is ignored */
 
58
                   break;
 
59
 
 
60
              default:
 
61
                   if (p->v != q->v || p->i != q->i) return 0;
 
62
                   break;
 
63
          }
 
64
     }
 
65
     A(0 /* can't happen */);
 
66
}
 
67
 
 
68
static int ok_twid(const twid *t, 
 
69
                   enum wakefulness wakefulness,
 
70
                   const tw_instr *q, INT n, INT r, INT m)
 
71
{
 
72
     return (wakefulness == t->wakefulness &&
 
73
             n == t->n &&
 
74
             r == t->r && 
 
75
             m <= t->m && 
 
76
             equal_instr(t->instr, q));
 
77
}
 
78
 
 
79
static twid *lookup(enum wakefulness wakefulness,
 
80
                    const tw_instr *q, INT n, INT r, INT m)
 
81
{
 
82
     twid *p;
 
83
 
 
84
     for (p = twlist[hash(n,r)]; 
 
85
          p && !ok_twid(p, wakefulness, q, n, r, m); 
 
86
          p = p->cdr)
 
87
          ;
 
88
     return p;
 
89
}
 
90
 
 
91
static INT twlen0(INT r, const tw_instr *p, INT *vl)
 
92
{
 
93
     INT ntwiddle = 0;
 
94
 
 
95
     /* compute length of bytecode program */
 
96
     A(r > 0);
 
97
     for ( ; p->op != TW_NEXT; ++p) {
 
98
          switch (p->op) {
 
99
              case TW_FULL:
 
100
                   ntwiddle += (r - 1) * 2;
 
101
                   break;
 
102
              case TW_HALF:
 
103
                   ntwiddle += (r - 1);
 
104
                   break;
 
105
              case TW_CEXP:
 
106
                   ntwiddle += 2;
 
107
                   break;
 
108
              case TW_COS:
 
109
              case TW_SIN:
 
110
                   ntwiddle += 1;
 
111
                   break;
 
112
          }
 
113
     }
 
114
 
 
115
     *vl = (INT)p->v;
 
116
     return ntwiddle;
 
117
}
 
118
 
 
119
INT X(twiddle_length)(INT r, const tw_instr *p)
 
120
{
 
121
     INT vl;
 
122
     return twlen0(r, p, &vl);
 
123
}
 
124
 
 
125
static R *compute(enum wakefulness wakefulness,
 
126
                  const tw_instr *instr, INT n, INT r, INT m)
 
127
{
 
128
     INT ntwiddle, j, vl;
 
129
     R *W, *W0;
 
130
     const tw_instr *p;
 
131
     triggen *t = X(mktriggen)(wakefulness, n);
 
132
 
 
133
     p = instr;
 
134
     ntwiddle = twlen0(r, p, &vl);
 
135
 
 
136
     A(m % vl == 0);
 
137
 
 
138
     W0 = W = (R *)MALLOC((ntwiddle * (m / vl)) * sizeof(R), TWIDDLES);
 
139
 
 
140
     for (j = 0; j < m; j += vl) {
 
141
          for (p = instr; p->op != TW_NEXT; ++p) {
 
142
               switch (p->op) {
 
143
                   case TW_FULL: {
 
144
                        INT i;
 
145
                        A(m * r <= n);
 
146
                        for (i = 1; i < r; ++i) {
 
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(m * r <= n);
 
167
                        t->cexp(t, (j + (INT)p->v) * (INT)p->i, d);
 
168
                        *W++ = d[0];
 
169
                        break;
 
170
 
 
171
                   }
 
172
 
 
173
                   case TW_SIN: {
 
174
                        R d[2];
 
175
 
 
176
                        A(m * r <= n);
 
177
                        t->cexp(t, (j + (INT)p->v) * (INT)p->i, d);
 
178
                        *W++ = d[1];
 
179
                        break;
 
180
                   }
 
181
 
 
182
                   case TW_CEXP:
 
183
                        A(m * r <= n);
 
184
                        t->cexp(t, (j + (INT)p->v) * (INT)p->i, W);
 
185
                        W += 2;
 
186
                        break;
 
187
               }
 
188
          }
 
189
     }
 
190
 
 
191
     X(triggen_destroy)(t);
 
192
     return W0;
 
193
}
 
194
 
 
195
static void mktwiddle(enum wakefulness wakefulness,
 
196
                      twid **pp, const tw_instr *instr, INT n, INT r, INT m)
 
197
{
 
198
     twid *p;
 
199
     INT h;
 
200
 
 
201
     if ((p = lookup(wakefulness, instr, n, r, m))) {
 
202
          ++p->refcnt;
 
203
     } else {
 
204
          p = (twid *) MALLOC(sizeof(twid), TWIDDLES);
 
205
          p->n = n;
 
206
          p->r = r;
 
207
          p->m = m;
 
208
          p->instr = instr;
 
209
          p->refcnt = 1;
 
210
          p->wakefulness = wakefulness;
 
211
          p->W = compute(wakefulness, instr, n, r, m);
 
212
 
 
213
          /* cons! onto twlist */
 
214
          h = hash(n, r);
 
215
          p->cdr = twlist[h];
 
216
          twlist[h] = p;
 
217
     }
 
218
 
 
219
     *pp = p;
 
220
}
 
221
 
 
222
static void twiddle_destroy(twid **pp)
 
223
{
 
224
     twid *p = *pp;
 
225
     twid **q;
 
226
 
 
227
     if ((--p->refcnt) == 0) {
 
228
          /* remove p from twiddle list */
 
229
          for (q = &twlist[hash(p->n, p->r)]; *q; q = &((*q)->cdr)) {
 
230
               if (*q == p) {
 
231
                    *q = p->cdr;
 
232
                    X(ifree)(p->W);
 
233
                    X(ifree)(p);
 
234
                    *pp = 0;
 
235
                    return;
 
236
               }
 
237
          }
 
238
          A(0 /* can't happen */ );
 
239
     }
 
240
}
 
241
 
 
242
 
 
243
void X(twiddle_awake)(enum wakefulness wakefulness, twid **pp, 
 
244
                      const tw_instr *instr, INT n, INT r, INT m)
 
245
{
 
246
     switch (wakefulness) {
 
247
         case SLEEPY: 
 
248
              twiddle_destroy(pp);
 
249
              break;
 
250
         default:
 
251
              mktwiddle(wakefulness, pp, instr, n, r, m);
 
252
              break;
 
253
     }
 
254
}
 
255
 
 
256
/* return a pointer to twiddles (0 if none) starting at mstart out of m */
 
257
const R *X(twiddle_shift)(const twid *p, INT mstart)
 
258
{
 
259
     if (p) {
 
260
          INT ntwiddle, vl;
 
261
 
 
262
          ntwiddle = twlen0(p->r, p->instr, &vl);
 
263
          A((mstart % vl) == 0);
 
264
          return (p->W + (mstart / vl) * ntwiddle);
 
265
     } else {
 
266
          return 0;
 
267
     }
 
268
}