Annotation of OpenXM_contrib2/asir2000/fft/dft.c, Revision 1.1.1.1
1.1 noro 1: /* $OpenXM: OpenXM/src/asir99/fft/dft.c,v 1.1.1.1 1999/11/10 08:12:27 noro Exp $ */
2: /*
3: * store arith modulus
4: * -------------------------------------------------------
5: * IEEE: float (23+1) <-> double (52+1) 24 bits
6: * int (32) <-> double (52+1) 26 bits
7: * 370: float (24) <-> double (56) 24/25 bits
8: * int (32) <-> double (56) 28 bits
9: */
10:
11:
12: #include "dft.h"
13:
14:
15: /*
16: #define C_DFT_FORE C_DFTfore
17: #define C_DFT_BACK C_DFTback
18: #define C_PREP_ALPHA C_prep_alpha
19: */
20:
21:
22: /*****************************************************/
23:
24: void C_DFT_FORE( in, nin, i1, K, powa,
25: #ifdef POWA_STRIDE
26: a1,
27: #endif
28: out, o1, P, Pinv, wk )
29: ModNum in[], powa[], out[], wk[];
30: ModNum P;
31: int nin, i1, K, o1;
32: #ifdef POWA_STRIDE
33: int a1;
34: #endif
35: double Pinv;
36: /*
37: * Let N = 2^K, and \alpha be a primitive N-th root of 1.
38: * out[0:N-1][o1] = DFT of in[0:nin-1][i1] using \alpha modulo P.
39: * The powers \alpha (to the degrees upto N/2=2^{K-1} at least) are given in powa[*],
40: * (or powa[a1**]).
41: * nin <= N assumed.
42: * wk[] must have 2^K entries at least.
43: */
44: {
45: int istep, ostep, i, j, k, idist, odist, n, Ndiv2n, K_k,
46: #ifndef OPT
47: K_k_1,
48: #endif
49: #ifdef POWA_STRIDE
50: aj,
51: #endif
52: inwk;
53: int idistwk, idistout, iostep, oostep;
54: ModNum *qi, *qo, *qis, *qos, a, t0, tn;
55:
56: /*
57: * for k = K-1, ..., 0 do the following:
58: *
59: * Out[2^{K-k}*i + j] <= t0 + \alpha^{2^k * j} * tn;
60: * Out[2^{K-k}*i + j + 2^{K-k-1}] <= t0 - \alpha^{2^k * j} * tn;
61: *
62: * where t0 = In[2^{K-k-1}*i + j] and
63: * tn = In[2^{K-k-1}*i + j + 2^{K-1}],
64: * for 0 <= i < 2^k and 0 <= j < 2^{K-k-1}.
65: *
66: */
67: /* In the following, n = 2^k and Ndiv2n = 2^{K-k-1}.
68: */
69: n = 1 << (k = K - 1);
70: idistwk = n, idistout = o1 << k;
71: /* when k = K-1, Out[2*i] <= In[i] + In[i+2^{K-1}],
72: * and Out[2*i+1] <= In[i] - In[i+2^{K-1}], for 0 <= i < 2^{K-1}.
73: */
74: if ( k&1 ) qo = wk, ostep = 2, inwk = 1, odist = 1;
75: else qo = out, ostep = o1 << 1, inwk = 0, odist = o1;
76: qi = in;
77: /**/
78: if ( nin <= n ) {
79: for ( i = nin; i > 0; i--, qi += i1, qo += ostep ) qo[0] = qo[odist] = qi[0];
80: for ( i = n - nin; i > 0; i--, qo += ostep ) qo[0] = qo[odist] = 0;
81: } else {
82: idist = i1 << k; /* distance = 2^{K-1} */
83: for ( j = i = nin - n; i > 0; i--, qi += i1, qo += ostep ) {
84: t0 = qi[0], tn = qi[idist];
85: if ( tn == (ModNum)0 ) qo[0] = qo[odist] = t0;
86: else {
87: qo[0] = AplusBmodP( t0, tn, P );
88: qo[odist] = A_BmodP( t0, tn, P );
89: }
90: }
91: for ( i = n - j; i > 0; i--, qi += i1, qo += ostep )
92: qo[0] = qo[odist] = qi[0];
93: }
94: /******/
95: #ifndef J_INNER
96: for ( K_k = 1, Ndiv2n = 1; --k > 0; inwk = 1 - inwk ) {
97: #ifndef OPT
98: K_k = (K_k_1 = K_k) + 1;
99: #endif
100: n >>= 1; /* == 2^k */
101: Ndiv2n <<= 1; /* == 2^{K-k-1} */
102: if ( inwk )
103: qi = wk, qo = out, iostep = 1, oostep = o1, idist = idistwk;
104: else
105: qi = out, qo = wk, iostep = o1, oostep = 1, idist = idistout;
106: qis = qi, qos = qo;
107: /**/
108: istep = ostep; /* = iostep << K_k_1; */
109: #ifndef OPT
110: ostep = oostep << K_k;
111: odist = oostep << K_k_1;
112: #else
113: odist = oostep << K_k;
114: ++K_k;
115: ostep = oostep << K_k;
116: #endif /* OPT */
117: /* j = 0 */
118: for ( i = n; i-- > 0; qi += istep, qo += ostep ) {
119: t0 = qi[0], tn = qi[idist];
120: if ( tn == (ModNum)0 ) qo[0] = qo[odist] = t0;
121: else {
122: qo[0] = AplusBmodP( t0, tn, P );
123: qo[odist] = A_BmodP( t0, tn, P );
124: }
125: }
126: /**/
127: #ifdef POWA_STRIDE
128: for ( aj = a1, j = 1; j < Ndiv2n; aj += a1, j++ ) {
129: a = powa[aj << k];
130: #else
131: for ( j = 1; j < Ndiv2n; j++ ) {
132: a = powa[j << k];
133: #endif
134: qi = (qis += iostep), qo = (qos += oostep);
135: for ( i = 0; i < n; i++, qi += istep, qo += ostep ) {
136: t0 = qi[0], tn = qi[idist];
137: if ( tn == (ModNum)0 ) qo[0] = qo[odist] = t0;
138: else {
139: AxBmodP( tn, ModNum, tn, a, P, Pinv );
140: qo[0] = AplusBmodP( t0, tn, P );
141: qo[odist] = A_BmodP( t0, tn, P );
142: }
143: }
144: }
145: }
146: /*
147: * When k = 0, for i = 0 (0 <= i < 2^0=1) and 0 <= j < 2^{K-1},
148: * Out[j] <= t0 + \alpha^{j}*tn and Out[j+2^{K-1}] <= t0 - \alpha^{j}*tn,
149: * where t0 = In[j] and tn = In[j + 2^{K-1}].
150: */
151: qi = wk, qo = out, idist = idistwk, odist = idistout;
152: /* == 2^{K-1}, == o1 << (K-1) */
153: t0 = qi[0], tn = qi[idist];
154: if ( tn == (ModNum)0 ) qo[0] = qo[odist] = t0;
155: else {
156: qo[0] = AplusBmodP( t0, tn, P );
157: qo[odist] = A_BmodP( t0, tn, P );
158: }
159: #ifdef POWA_STRIDE
160: for ( k = o1, aj = a1, j = 1; j < idist; aj += a1, j++ ) {
161: #else
162: for ( k = o1, j = 1; j < idist; j++ ) {
163: #endif
164: qi++, qo += k;
165: t0 = qi[0], tn = qi[idist];
166: if ( tn == (ModNum)0 ) qo[0] = qo[odist] = t0;
167: else {
168: #ifdef POWA_STRIDE
169: a = powa[aj];
170: #else
171: a = powa[j];
172: #endif
173: AxBmodP( tn, ModNum, tn, a, P, Pinv );
174: qo[0] = AplusBmodP( t0, tn, P );
175: qo[odist] = A_BmodP( t0, tn, P );
176: }
177: }
178: #else
179: for ( K_k = 1, Ndiv2n = 1; --k >= 0; inwk = 1 - inwk ) {
180: #ifndef OPT
181: K_k = (K_k_1 = K_k) + 1;
182: #endif
183: n >>= 1; /* == 2^k */
184: Ndiv2n <<= 1; /* == 2^{K-k-1} */
185: if ( inwk ) qi = wk, qo = out, istep = 1, ostep = o1, idist = idistwk;
186: else qi = out, qo = wk, istep = o1, ostep = 1, idist = idistout;
187: qis = qi, qos = qo;
188: /**/
189: iostep = oostep; /* = istep << K_k_1; */
190: #ifndef OPT
191: oostep = ostep << K_k;
192: odist = ostep << K_k_1;
193: #else
194: odist = ostep << K_k;
195: ++K_k;
196: oostep = ostep << K_k;
197: #endif /* OPT */
198: for ( i = n; i-- > 0; qis += iostep, qos += oostep ) {
199: /* j = 0 */
200: t0 = qis[0], tn = qis[idist];
201: if ( tn == (ModNum)0 ) qos[0] = qos[odist] = t0;
202: else {
203: qos[0] = AplusBmodP( t0, tn, P );
204: qos[odist] = A_BmodP( t0, tn, P );
205: }
206: #ifdef POWA_STRIDE
207: for ( aj = a1, j = 1, qi = qis, qo = qos; j < Ndiv2n; aj += a1, j++ ) {
208: #else
209: for ( j = 1, qi = qis, qo = qos; j < Ndiv2n; j++ ) {
210: #endif
211: qi += istep, qo += ostep;
212: t0 = qi[0], tn = qi[idist];
213: if ( tn == (ModNum)0 ) qo[0] = qo[odist] = t0;
214: else {
215: #ifdef POWA_STRIDE
216: a = powa[aj << k];
217: #else
218: a = powa[j << k];
219: #endif
220: AxBmodP( tn, ModNum, tn, a, P, Pinv );
221: qo[0] = AplusBmodP( t0, tn, P );
222: qo[odist] = A_BmodP( t0, tn, P );
223: }
224: }
225: }
226: }
227: #endif
228: }
229:
230: /*********************************************************************
231: *
232: * NOTE:
233: * (1) Let \alpha be a primitive N-th root of unity modulo P,
234: * where N is even,
235: * and t be an integer s.t. 0 <= t < N/2.
236: * Then, 1/\alpha^t is given by (- \alpha^s) with s = N/2 - t.
237: *
238: * (2) Let P be a prime s.t. P = 2^n*Q + 1, where Q is an odd integer.
239: * Then, for 0 < s <= n, 2^{-s} \equiv - 2^{n-s}*Q \bmod P,
240: * because 1 \equiv - 2^n*Q \bmod P.
241: *
242: **********************************************************************/
243:
244: void C_DFT_BACK( in, N, i1, log2N, powa,
245: #ifdef POWA_STRIDE
246: a1,
247: #endif
248: out, o1, osi, nout, Ninv, P, Pinv, wk )
249: ModNum in[], powa[], out[], wk[];
250: int N, log2N, osi, nout, i1, o1;
251: #ifdef POWA_STRIDE
252: int a1;
253: #endif
254: ModNum P, Ninv;
255: double Pinv;
256: /*
257: * Let K denote log2N. Let N = 2^K, and \alpha be a primitive N-th root of 1.
258: * This routine computes the inverse discrete-Fourier transform of in[0:N-1][i1]
259: * using \alpha^{-1} modulo P, and store the `nout' transformed data from the
260: * `osi'-th in out[0:nout-1][o1]. 0 <= osi < N-1 and 0 < nout <= N.
261: * The powers of \alpha (to the degrees upto N/2=2^{K-1} at least) are given in powa[*] (or powa[a1* *]).
262: * The area wk[] is used to contain the intermediate transformed data, and thus,
263: * must have 2*N entries at least. Notice that `out' cannot be used because
264: * its space amount (`nout') may not be sufficient for this purpose.
265: */
266: {
267: int i, j, k, n, is, os, istep, ostep, idist, inwk, halfN;
268: #ifdef POWA_STRIDE
269: int aj, halfNa;
270: #endif
271: ModNum *qi, *qo, *qis, *qos, a, t0, tn, tt;
272:
273: /*
274: * for k = K-1, ..., 0 do the following:
275: *
276: * Out[2^{K-k}*i + j] <= t0 + \alpha^{- 2^k * j} * tn;
277: * = t0 - \alpha^{N/2 - 2^k * j}*tn;
278: * Out[2^{K-k}*i + j + 2^{K-k-1}] <= t0 - \alpha^{- 2^k * j} * tn;
279: * = t0 + \alpha^{N/2 - 2^k * j}*tn;
280: *
281: * where t0 = In[2^{K-k-1}*i + j] and
282: * tn = In[2^{K-k-1}*i + j + 2^{K-1}],
283: * for 0 <= i < 2^k and 0 <= j < 2^{K-k-1}.
284: *
285: */
286: if ( log2N == 1 ) {
287: /* K = 1, k = 0, 0 <= i < 1, 0 <= j < 1.
288: * Out[0] <= t0 + tn, Out[1] <= t0 - tn,
289: * where t0 = In[0] and Tn = In[1].
290: */
291: t0 = in[0], tn = in[i1];
292: if ( osi == 0 ) {
293: /* out[0] = (t0 + tn)*Ninv; */
294: tt = tn == (ModNum)0 ? t0 : AplusBmodP( t0, tn, P );
295: if ( tt != (ModNum)0 ) {
296: AxBmodP( tt, ModNum, tt, Ninv, P, Pinv );
297: }
298: out[0] = tt;
299: i = 1;
300: } else i = 0;
301: /**/
302: if ( osi + nout >= 2 ) {
303: /* out[osi == 0 ? 1 : 0] = (t0 - tn)*Ninv; */
304: tt = tn == (ModNum)0 ? t0 : A_BmodP( t0, tn, P );
305: if ( tt != (ModNum)0 ) {
306: AxBmodP( tt, ModNum, tt, Ninv, P, Pinv );
307: }
308: out[i] = tt;
309: }
310: return;
311: }
312: /****/
313: halfN = n = 1 << (k = log2N - 1); /* == 2^{K-1} */
314: #ifdef POWA_STRIDE
315: halfNa = a1 << k;
316: #endif
317: /*
318: * when k = K-1,
319: * Out[2*i] <= t0 + tn, and Out[2*i+1] <= t0 - tn,
320: * where t0 = In[i] and tn = In[i+2^{K-1}],
321: * for 0 <= i < 2^{K-1}.
322: */
323: qi = in, istep = i1, idist = i1 << k, qo = wk, inwk = 0;
324: for ( i = n; i-- > 0; qi += istep, qo += 2 ) {
325: t0 = qi[0], tn = qi[idist];
326: if ( tn == (ModNum)0 ) qo[0] = qo[1] = t0;
327: else {
328: qo[0] = AplusBmodP( t0, tn, P );
329: qo[1] = A_BmodP( t0, tn, P );
330: }
331: }
332: #ifdef EBUG
333: fprintf( stderr, "::: DFT^{-1} ::: after the first step\n" );
334: for ( qo = wk, i = 0; i < N; i++ ) {
335: if ( (i%5) == 0 ) fprintf( stderr, "\t" );
336: fprintf( stderr, "%10d,", (int)wk[i] );
337: if ( (i%5) == 4 ) fprintf( stderr, "\n" );
338: }
339: if ( (i%5) != 0 ) fprintf( stderr, "\n" );
340: #endif
341: /**/
342: idist = halfN;
343: for ( ostep = 2; --k > 0; inwk = 1 - inwk ) {
344: n >>= 1;
345: istep = ostep; /* == 2^{K-k-1} */
346: ostep <<= 1; /* == 2^{K-k} */
347: if ( inwk ) qi = &wk[N], qo = wk;
348: else qi = wk, qo = &wk[N];
349: qis = qi, qos = qo;
350: /*
351: * for j = 0,
352: * Out[2^{K-k}*i] <= t0 + tn, and Out[2^{K-k}*i + 2^{K-k-1}] <= t0 - tn,
353: * where t0 = In[2^{K-k-1}*i] and tn = In[2^{K-k-1}*i + 2^{K-1}].
354: */
355: for ( i = n, is = os = 0; i-- > 0; qi += istep, qo += ostep ) {
356: t0 = qi[0], tn = qi[idist];
357: if ( tn == (ModNum)0 ) qo[0] = qo[istep] = t0;
358: else {
359: qo[0] = AplusBmodP( t0, tn, P );
360: qo[istep] = A_BmodP( t0, tn, P );
361: }
362: }
363: #ifdef POWA_STRIDE
364: for ( aj = a1, j = 1; j < istep; aj += a1, j++ ) {
365: /*** a = P - powa[halfNa - (aj << k)]; ***/
366: a = powa[halfNa - (aj << k)];
367: #else
368: for ( j = 1; j < istep; j++ ) {
369: /*** a = P - powa[halfN - (j << k)]; ***/
370: a = powa[halfN - (j << k)];
371: #endif
372: qi = ++qis, qo = ++qos;
373: for ( i = n; i-- > 0; qi += istep, qo += ostep ) {
374: t0 = qi[0], tn = qi[idist];
375: if ( tn == (ModNum)0 ) qo[0] = qo[istep] = t0;
376: else {
377: AxBmodP( tn, ModNum, tn, a, P, Pinv );
378: /*** qo[0] = AplusBmodP( t0, tn, P ); ***/
379: /*** qo[istep] = A_BmodP( t0, tn, P ); ***/
380: qo[0] = A_BmodP( t0, tn, P );
381: qo[istep] = AplusBmodP( t0, tn, P );
382: }
383: }
384: }
385: }
386: #if 0
387: final:
388: #endif
389: /*
390: * The final step of k = 0. i can take only the value 0.
391: *
392: * Out[j] <= t0 + \alpha^{-j} * tn;
393: * Out[j + 2^{K-1}] <= t0 - \alpha^{-j} * tn;
394: *
395: * where t0 = In[j] and tn = In[j + 2^{K-1}],
396: * for 0 <= j < 2^{K-1}.
397: *
398: */
399: qo = out, qi = &wk[inwk ? N : 0];
400: if ( (n = osi + nout) > N ) nout = (n = N) - osi; /* true nout */
401: if ( (k = nout - halfN) <= 0 ) { /* no overlap, i.e., no +- */
402: if ( osi <= 0 ) {
403: t0 = qi[0], tn = qi[idist];
404: if ( tn != (ModNum)0 ) t0 = AplusBmodP( t0, tn, P );
405: if ( t0 == (ModNum)0 ) qo[0] = t0;
406: else {
407: AxBmodP( qo[0], ModNum, t0, Ninv, P, Pinv );
408: }
409: #ifdef POWA_STRIDE
410: aj = a1;
411: #endif
412: j = 1, k = n, i = 0, qo += o1;
413: } else {
414: j = osi;
415: #ifdef POWA_STRIDE
416: aj = osi*a1;
417: #endif
418: if ( n <= halfN ) i = 0, k = n; /* lower only */
419: else if ( j == halfN ) goto L_halfN; /* start with j = N/2 */
420: else if ( j > halfN ) goto L_upper; /* upper only */
421: else i = n - halfN + 1, k = halfN; /* we have lower, and i upper */
422: }
423: } else {
424: os = o1*halfN;
425: #ifdef EBUG
426: fprintf( stderr, "::::: DFT^{-1} ::: os=%d*%d=%d\n", o1, halfN, os );
427: #endif
428: if ( osi <= 0 ) {
429: t0 = qi[0], tn = qi[idist];
430: if ( tn == (ModNum)0 ) {
431: if ( t0 == (ModNum)0 ) tt = t0;
432: else { AxBmodP( tt, ModNum, t0, Ninv, P, Pinv ); }
433: qo[0] = qo[os] = tt;
434: } else {
435: tt = AplusBmodP( t0, tn, P );
436: if ( tt == (ModNum)0 ) qo[0] = tt;
437: else { AxBmodP( qo[0], ModNum, tt, Ninv, P, Pinv ); }
438: tt = A_BmodP( t0, tn, P );
439: if ( tt == (ModNum)0 ) qo[os] = tt;
440: else { AxBmodP( qo[os], ModNum, tt, Ninv, P, Pinv ); }
441: }
442: #ifdef EBUG
443: fprintf( stderr, "::::: DFT^{-1} ::: [0]=%d, [%d]=%d\n", (int)qo[0], os, (int)qo[os] );
444: #endif
445: #ifdef POWA_STRIDE
446: j = 1, aj = a1, n = halfN, i = 0, qo += o1, k--;
447: } else j = osi, aj = osi*a1, i = (n -= k) - halfN;
448: /**/
449: for ( ; k-- > 0; aj += a1, j++, qo += o1 ) {
450: #else
451: j = 1, n = halfN, i = 0, qo += o1, k--;
452: } else j = osi, i = (n -= k) - halfN;
453: /**/
454: for ( ; k-- > 0; j++, qo += o1 ) {
455: #endif
456: t0 = qi[j], tn = qi[j+idist];
457: if ( tn == (ModNum)0 ) {
458: if ( t0 != (ModNum)0 ) {
459: AxBmodP( t0, ModNum, t0, Ninv, P, Pinv );
460: }
461: qo[0] = qo[os] = t0;
462: } else {
463: #ifdef POWA_STRIDE
464: /*** a = P - powa[halfNa - aj]; ***/
465: a = powa[halfNa - aj];
466: #else
467: /*** a = P - powa[halfN - j]; ***/
468: a = powa[halfN - j];
469: #endif
470: AxBmodP( tn, ModNum, tn, a, P, Pinv );
471: /*** tt = AplusBmodP( t0, tn, P ); ***/
472: tt = A_BmodP( t0, tn, P );
473: if ( tt == (ModNum)0 ) qo[0] = tt;
474: else { AxBmodP( qo[0], ModNum, tt, Ninv, P, Pinv ); }
475: /*** tt = A_BmodP( t0, tn, P ); ***/
476: tt = AplusBmodP( t0, tn, P );
477: if ( tt == (ModNum)0 ) qo[os] = tt;
478: else { AxBmodP( qo[os], ModNum, tt, Ninv, P, Pinv ); }
479: }
480: }
481: k = halfN;
482: }
483: /*
484: * At this point, k is an upper limit < N/2 for j, i is the number of j's
485: * > N/2, and n is the real upper limit of j.
486: */
487: #ifdef POWA_STRIDE
488: for ( ; j < k; aj += a1, j++, qo += o1 ) {
489: #else
490: for ( ; j < k; j++, qo += o1 ) {
491: #endif
492: t0 = qi[j], tn = qi[j+idist];
493: if ( tn == (ModNum)0 ) {
494: if ( t0 == (ModNum)0 ) qo[0] = t0;
495: else { AxBmodP( qo[0], ModNum, t0, Ninv, P, Pinv ); }
496: } else {
497: #ifdef POWA_STRIDE
498: a = P - powa[halfNa - aj];
499: #else
500: a = P - powa[halfN - j];
501: #endif
502: AxBplusCmodP( tt, ModNum, tn, a, t0, P, Pinv );
503: if ( tt == (ModNum)0 ) qo[0] = tt;
504: else { AxBmodP( qo[0], ModNum, tt, Ninv, P, Pinv ); }
505: }
506: }
507: if ( i <= 0 ) return;
508: /**/
509: L_halfN:
510: t0 = qi[0], tn = qi[idist];
511: tt = tn == (ModNum)0 ? t0 : A_BmodP( t0, tn, P );
512: if ( tt == (ModNum)0 ) qo[0] = tt;
513: else { AxBmodP( qo[0], ModNum, tt, Ninv, P, Pinv ); }
514: qo += o1, j++;
515: /**/
516: L_upper:
517: #ifdef POWA_STRIDE
518: for ( n -= halfN, aj = a1, j = 1; j < n; aj += a1, j++, qo += o1 ) {
519: #else
520: for ( n -= halfN, j = 1; j < n; j++, qo += o1 ) {
521: #endif
522: t0 = qi[j], tn = qi[j+idist];
523: if ( tn == (ModNum)0 ) {
524: if ( t0 == (ModNum)0 ) qo[0] = t0;
525: else { AxBmodP( qo[0], ModNum, t0, Ninv, P, Pinv ); }
526: } else {
527: #ifdef POWA_STRIDE
528: a = powa[halfNa - aj];
529: #else
530: a = powa[halfN - j];
531: #endif
532: AxBplusCmodP( tt, ModNum, tn, a, t0, P, Pinv );
533: if ( tt == (ModNum)0 ) qo[0] = tt;
534: else { AxBmodP( qo[0], ModNum, tt, Ninv, P, Pinv ); }
535: }
536: }
537: }
538:
539: #if USE_FLOAT
540: void C_PREP_ALPHA( r, log2ord, log2k, n, tbl, P, Pinv )
541: ModNum r, tbl[], P;
542: int log2ord, log2k, n;
543: double Pinv;
544: /*
545: * Let K and k denote log2ord and log2k, respectively.
546: * Let r be a primitive (2^K)-th root of unity in Z/(P), where P is a prime.
547: * Compute a = r^{2^{K-k}}, a primitive (2^k)-th root of unity, and
548: * its powers to the degrees upto n in tbl[*].
549: */
550: {
551: int i;
552: ModNum *q;
553: double t, w, dp = (double) P;
554:
555: for ( w = (double) r, i = log2ord - log2k; i > 0; i-- ) {
556: w *= w;
557: if ( w < dp ) continue;
558: w -= (dp*(double)((int) (w*Pinv)));
559: }
560: tbl[0] = 1;
561: tbl[1] = (ModNum)(t = w);
562: for ( q = &tbl[2], i = n - 1; i > 0; i-- ) {
563: t *= w;
564: if ( t >= dp ) t -= (dp*(double)((int) (t*Pinv)));
565: *q++ = t;
566: }
567: }
568: #else
569:
570: void C_PREP_ALPHA( r, log2ord, log2k, n, tbl, P, Pinv )
571: ModNum r, tbl[], P;
572: int log2ord, log2k, n;
573: double Pinv;
574: /*
575: * Let K and k denote log2ord and log2k, respectively.
576: * Let r be a primitive (2^K)-th root of unity in Z/(P), where P is a prime.
577: * Compute a = r^{2^{K-k}}, a primitive (2^k)-th root of unity, and
578: * its powers to the degrees upto n in tbl[*].
579: */
580: {
581: int i;
582: ModNum *q;
583: ModNum t, w, s;
584:
585: for ( w = r, i = log2ord - log2k; i > 0; i-- ) {
586: AxBmodP( t, ModNum, w, w, P, Pinv );
587: w = t;
588: }
589: tbl[0] = 1;
590: tbl[1] = (ModNum)(t = w);
591: for ( q = &tbl[2], i = n - 1; i > 0; i-- ) {
592: AxBmodP( s, ModNum, t, w, P, Pinv );
593: t = s;
594: *q++ = t;
595: }
596: }
597: #endif
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to dft.c CVS log ![[TXT]](/icons/cvsweb/text.gif){kind=icon}
![[DIR]](/icons/cvsweb/dir.gif){kind=icon}
Up to [local] / OpenXM_contrib2 / asir2000 / fft