Annotation of OpenXM_contrib2/asir2018/fft/dft.c, Revision 1.2
1.1 noro 1: /*
2: * Copyright (c) 1994-2000 FUJITSU LABORATORIES LIMITED
3: * All rights reserved.
4: *
5: * FUJITSU LABORATORIES LIMITED ("FLL") hereby grants you a limited,
6: * non-exclusive and royalty-free license to use, copy, modify and
7: * redistribute, solely for non-commercial and non-profit purposes, the
8: * computer program, "Risa/Asir" ("SOFTWARE"), subject to the terms and
9: * conditions of this Agreement. For the avoidance of doubt, you acquire
10: * only a limited right to use the SOFTWARE hereunder, and FLL or any
11: * third party developer retains all rights, including but not limited to
12: * copyrights, in and to the SOFTWARE.
13: *
14: * (1) FLL does not grant you a license in any way for commercial
15: * purposes. You may use the SOFTWARE only for non-commercial and
16: * non-profit purposes only, such as academic, research and internal
17: * business use.
18: * (2) The SOFTWARE is protected by the Copyright Law of Japan and
19: * international copyright treaties. If you make copies of the SOFTWARE,
20: * with or without modification, as permitted hereunder, you shall affix
21: * to all such copies of the SOFTWARE the above copyright notice.
22: * (3) An explicit reference to this SOFTWARE and its copyright owner
23: * shall be made on your publication or presentation in any form of the
24: * results obtained by use of the SOFTWARE.
25: * (4) In the event that you modify the SOFTWARE, you shall notify FLL by
26: * e-mail at risa-admin@sec.flab.fujitsu.co.jp of the detailed specification
27: * for such modification or the source code of the modified part of the
28: * SOFTWARE.
29: *
30: * THE SOFTWARE IS PROVIDED AS IS WITHOUT ANY WARRANTY OF ANY KIND. FLL
31: * MAKES ABSOLUTELY NO WARRANTIES, EXPRESSED, IMPLIED OR STATUTORY, AND
32: * EXPRESSLY DISCLAIMS ANY IMPLIED WARRANTY OF MERCHANTABILITY, FITNESS
33: * FOR A PARTICULAR PURPOSE OR NONINFRINGEMENT OF THIRD PARTIES'
34: * RIGHTS. NO FLL DEALER, AGENT, EMPLOYEES IS AUTHORIZED TO MAKE ANY
35: * MODIFICATIONS, EXTENSIONS, OR ADDITIONS TO THIS WARRANTY.
36: * UNDER NO CIRCUMSTANCES AND UNDER NO LEGAL THEORY, TORT, CONTRACT,
37: * OR OTHERWISE, SHALL FLL BE LIABLE TO YOU OR ANY OTHER PERSON FOR ANY
38: * DIRECT, INDIRECT, SPECIAL, INCIDENTAL, PUNITIVE OR CONSEQUENTIAL
39: * DAMAGES OF ANY CHARACTER, INCLUDING, WITHOUT LIMITATION, DAMAGES
40: * ARISING OUT OF OR RELATING TO THE SOFTWARE OR THIS AGREEMENT, DAMAGES
41: * FOR LOSS OF GOODWILL, WORK STOPPAGE, OR LOSS OF DATA, OR FOR ANY
42: * DAMAGES, EVEN IF FLL SHALL HAVE BEEN INFORMED OF THE POSSIBILITY OF
43: * SUCH DAMAGES, OR FOR ANY CLAIM BY ANY OTHER PARTY. EVEN IF A PART
44: * OF THE SOFTWARE HAS BEEN DEVELOPED BY A THIRD PARTY, THE THIRD PARTY
45: * DEVELOPER SHALL HAVE NO LIABILITY IN CONNECTION WITH THE USE,
46: * PERFORMANCE OR NON-PERFORMANCE OF THE SOFTWARE.
47: *
1.2 ! noro 48: * $OpenXM: OpenXM_contrib2/asir2018/fft/dft.c,v 1.1 2018/09/19 05:45:07 noro Exp $
1.1 noro 49: */
50: /*
51: * store arith modulus
52: * -------------------------------------------------------
53: * IEEE: float (23+1) <-> double (52+1) 24 bits
54: * int (32) <-> double (52+1) 26 bits
55: * 370: float (24) <-> double (56) 24/25 bits
56: * int (32) <-> double (56) 28 bits
57: */
58:
59:
60: #include "dft.h"
61:
62: /*
63: #define C_DFT_FORE C_DFTfore
64: #define C_DFT_BACK C_DFTback
65: #define C_PREP_ALPHA C_prep_alpha
66: */
67:
68:
69: /*****************************************************/
70:
71: void C_DFT_FORE( in, nin, i1, K, powa,
72: #ifdef POWA_STRIDE
73: a1,
74: #endif
75: out, o1, P, Pinv, wk )
76: ModNum in[], powa[], out[], wk[];
77: ModNum P;
78: int nin, i1, K, o1;
79: #ifdef POWA_STRIDE
80: int a1;
81: #endif
82: double Pinv;
83: /*
84: * Let N = 2^K, and \alpha be a primitive N-th root of 1.
85: * out[0:N-1][o1] = DFT of in[0:nin-1][i1] using \alpha modulo P.
86: * The powers \alpha (to the degrees upto N/2=2^{K-1} at least) are given in powa[*],
87: * (or powa[a1**]).
88: * nin <= N assumed.
89: * wk[] must have 2^K entries at least.
90: */
91: {
92: int istep, ostep, i, j, k, idist, odist, n, Ndiv2n, K_k,
93: #ifndef OPT
94: K_k_1,
95: #endif
96: #ifdef POWA_STRIDE
97: aj,
98: #endif
99: inwk;
100: int idistwk, idistout, iostep, oostep;
101: ModNum *qi, *qo, *qis, *qos, a, t0, tn;
102:
103: /*
104: * for k = K-1, ..., 0 do the following:
105: *
106: * Out[2^{K-k}*i + j] <= t0 + \alpha^{2^k * j} * tn;
107: * Out[2^{K-k}*i + j + 2^{K-k-1}] <= t0 - \alpha^{2^k * j} * tn;
108: *
109: * where t0 = In[2^{K-k-1}*i + j] and
110: * tn = In[2^{K-k-1}*i + j + 2^{K-1}],
111: * for 0 <= i < 2^k and 0 <= j < 2^{K-k-1}.
112: *
113: */
114: /* In the following, n = 2^k and Ndiv2n = 2^{K-k-1}.
115: */
116: n = 1 << (k = K - 1);
117: idistwk = n, idistout = o1 << k;
118: /* when k = K-1, Out[2*i] <= In[i] + In[i+2^{K-1}],
119: * and Out[2*i+1] <= In[i] - In[i+2^{K-1}], for 0 <= i < 2^{K-1}.
120: */
121: if ( k&1 ) qo = wk, ostep = 2, inwk = 1, odist = 1;
122: else qo = out, ostep = o1 << 1, inwk = 0, odist = o1;
123: qi = in;
124: /**/
125: if ( nin <= n ) {
126: for ( i = nin; i > 0; i--, qi += i1, qo += ostep ) qo[0] = qo[odist] = qi[0];
127: for ( i = n - nin; i > 0; i--, qo += ostep ) qo[0] = qo[odist] = 0;
128: } else {
129: idist = i1 << k; /* distance = 2^{K-1} */
130: for ( j = i = nin - n; i > 0; i--, qi += i1, qo += ostep ) {
131: t0 = qi[0], tn = qi[idist];
132: if ( tn == (ModNum)0 ) qo[0] = qo[odist] = t0;
133: else {
134: qo[0] = AplusBmodP( t0, tn, P );
135: qo[odist] = A_BmodP( t0, tn, P );
136: }
137: }
138: for ( i = n - j; i > 0; i--, qi += i1, qo += ostep )
139: qo[0] = qo[odist] = qi[0];
140: }
141: /******/
142: #ifndef J_INNER
143: for ( K_k = 1, Ndiv2n = 1; --k > 0; inwk = 1 - inwk ) {
144: #ifndef OPT
145: K_k = (K_k_1 = K_k) + 1;
146: #endif
147: n >>= 1; /* == 2^k */
148: Ndiv2n <<= 1; /* == 2^{K-k-1} */
149: if ( inwk )
150: qi = wk, qo = out, iostep = 1, oostep = o1, idist = idistwk;
151: else
152: qi = out, qo = wk, iostep = o1, oostep = 1, idist = idistout;
153: qis = qi, qos = qo;
154: /**/
155: istep = ostep; /* = iostep << K_k_1; */
156: #ifndef OPT
157: ostep = oostep << K_k;
158: odist = oostep << K_k_1;
159: #else
160: odist = oostep << K_k;
161: ++K_k;
162: ostep = oostep << K_k;
163: #endif /* OPT */
164: /* j = 0 */
165: for ( i = n; i-- > 0; qi += istep, qo += ostep ) {
166: t0 = qi[0], tn = qi[idist];
167: if ( tn == (ModNum)0 ) qo[0] = qo[odist] = t0;
168: else {
169: qo[0] = AplusBmodP( t0, tn, P );
170: qo[odist] = A_BmodP( t0, tn, P );
171: }
172: }
173: /**/
174: #ifdef POWA_STRIDE
175: for ( aj = a1, j = 1; j < Ndiv2n; aj += a1, j++ ) {
176: a = powa[aj << k];
177: #else
178: for ( j = 1; j < Ndiv2n; j++ ) {
179: a = powa[j << k];
180: #endif
181: qi = (qis += iostep), qo = (qos += oostep);
182: for ( i = 0; i < n; i++, qi += istep, qo += ostep ) {
183: t0 = qi[0], tn = qi[idist];
184: if ( tn == (ModNum)0 ) qo[0] = qo[odist] = t0;
185: else {
186: AxBmodP( tn, ModNum, tn, a, P, Pinv );
187: qo[0] = AplusBmodP( t0, tn, P );
188: qo[odist] = A_BmodP( t0, tn, P );
189: }
190: }
191: }
192: }
193: /*
194: * When k = 0, for i = 0 (0 <= i < 2^0=1) and 0 <= j < 2^{K-1},
195: * Out[j] <= t0 + \alpha^{j}*tn and Out[j+2^{K-1}] <= t0 - \alpha^{j}*tn,
196: * where t0 = In[j] and tn = In[j + 2^{K-1}].
197: */
198: qi = wk, qo = out, idist = idistwk, odist = idistout;
199: /* == 2^{K-1}, == o1 << (K-1) */
200: t0 = qi[0], tn = qi[idist];
201: if ( tn == (ModNum)0 ) qo[0] = qo[odist] = t0;
202: else {
203: qo[0] = AplusBmodP( t0, tn, P );
204: qo[odist] = A_BmodP( t0, tn, P );
205: }
206: #ifdef POWA_STRIDE
207: for ( k = o1, aj = a1, j = 1; j < idist; aj += a1, j++ ) {
208: #else
209: for ( k = o1, j = 1; j < idist; j++ ) {
210: #endif
211: qi++, qo += k;
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];
217: #else
218: a = powa[j];
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: #else
226: for ( K_k = 1, Ndiv2n = 1; --k >= 0; inwk = 1 - inwk ) {
227: #ifndef OPT
228: K_k = (K_k_1 = K_k) + 1;
229: #endif
230: n >>= 1; /* == 2^k */
231: Ndiv2n <<= 1; /* == 2^{K-k-1} */
232: if ( inwk ) qi = wk, qo = out, istep = 1, ostep = o1, idist = idistwk;
233: else qi = out, qo = wk, istep = o1, ostep = 1, idist = idistout;
234: qis = qi, qos = qo;
235: /**/
236: iostep = oostep; /* = istep << K_k_1; */
237: #ifndef OPT
238: oostep = ostep << K_k;
239: odist = ostep << K_k_1;
240: #else
241: odist = ostep << K_k;
242: ++K_k;
243: oostep = ostep << K_k;
244: #endif /* OPT */
245: for ( i = n; i-- > 0; qis += iostep, qos += oostep ) {
246: /* j = 0 */
247: t0 = qis[0], tn = qis[idist];
248: if ( tn == (ModNum)0 ) qos[0] = qos[odist] = t0;
249: else {
250: qos[0] = AplusBmodP( t0, tn, P );
251: qos[odist] = A_BmodP( t0, tn, P );
252: }
253: #ifdef POWA_STRIDE
254: for ( aj = a1, j = 1, qi = qis, qo = qos; j < Ndiv2n; aj += a1, j++ ) {
255: #else
256: for ( j = 1, qi = qis, qo = qos; j < Ndiv2n; j++ ) {
257: #endif
258: qi += istep, qo += ostep;
259: t0 = qi[0], tn = qi[idist];
260: if ( tn == (ModNum)0 ) qo[0] = qo[odist] = t0;
261: else {
262: #ifdef POWA_STRIDE
263: a = powa[aj << k];
264: #else
265: a = powa[j << k];
266: #endif
267: AxBmodP( tn, ModNum, tn, a, P, Pinv );
268: qo[0] = AplusBmodP( t0, tn, P );
269: qo[odist] = A_BmodP( t0, tn, P );
270: }
271: }
272: }
273: }
274: #endif
275: }
276:
277: /*********************************************************************
278: *
279: * NOTE:
280: * (1) Let \alpha be a primitive N-th root of unity modulo P,
281: * where N is even,
282: * and t be an integer s.t. 0 <= t < N/2.
283: * Then, 1/\alpha^t is given by (- \alpha^s) with s = N/2 - t.
284: *
285: * (2) Let P be a prime s.t. P = 2^n*Q + 1, where Q is an odd integer.
286: * Then, for 0 < s <= n, 2^{-s} \equiv - 2^{n-s}*Q \bmod P,
287: * because 1 \equiv - 2^n*Q \bmod P.
288: *
289: **********************************************************************/
290:
291: void C_DFT_BACK( in, N, i1, log2N, powa,
292: #ifdef POWA_STRIDE
293: a1,
294: #endif
295: out, o1, osi, nout, Ninv, P, Pinv, wk )
296: ModNum in[], powa[], out[], wk[];
297: int N, log2N, osi, nout, i1, o1;
298: #ifdef POWA_STRIDE
299: int a1;
300: #endif
301: ModNum P, Ninv;
302: double Pinv;
303: /*
304: * Let K denote log2N. Let N = 2^K, and \alpha be a primitive N-th root of 1.
305: * This routine computes the inverse discrete-Fourier transform of in[0:N-1][i1]
306: * using \alpha^{-1} modulo P, and store the `nout' transformed data from the
307: * `osi'-th in out[0:nout-1][o1]. 0 <= osi < N-1 and 0 < nout <= N.
308: * The powers of \alpha (to the degrees upto N/2=2^{K-1} at least) are given in powa[*] (or powa[a1* *]).
309: * The area wk[] is used to contain the intermediate transformed data, and thus,
310: * must have 2*N entries at least. Notice that `out' cannot be used because
311: * its space amount (`nout') may not be sufficient for this purpose.
312: */
313: {
314: int i, j, k, n, is, os, istep, ostep, idist, inwk, halfN;
315: #ifdef POWA_STRIDE
316: int aj, halfNa;
317: #endif
318: ModNum *qi, *qo, *qis, *qos, a, t0, tn, tt;
319:
320: /*
321: * for k = K-1, ..., 0 do the following:
322: *
323: * Out[2^{K-k}*i + j] <= t0 + \alpha^{- 2^k * j} * tn;
324: * = t0 - \alpha^{N/2 - 2^k * j}*tn;
325: * Out[2^{K-k}*i + j + 2^{K-k-1}] <= t0 - \alpha^{- 2^k * j} * tn;
326: * = t0 + \alpha^{N/2 - 2^k * j}*tn;
327: *
328: * where t0 = In[2^{K-k-1}*i + j] and
329: * tn = In[2^{K-k-1}*i + j + 2^{K-1}],
330: * for 0 <= i < 2^k and 0 <= j < 2^{K-k-1}.
331: *
332: */
333: if ( log2N == 1 ) {
334: /* K = 1, k = 0, 0 <= i < 1, 0 <= j < 1.
335: * Out[0] <= t0 + tn, Out[1] <= t0 - tn,
336: * where t0 = In[0] and Tn = In[1].
337: */
338: t0 = in[0], tn = in[i1];
339: if ( osi == 0 ) {
340: /* out[0] = (t0 + tn)*Ninv; */
341: tt = tn == (ModNum)0 ? t0 : AplusBmodP( t0, tn, P );
342: if ( tt != (ModNum)0 ) {
343: AxBmodP( tt, ModNum, tt, Ninv, P, Pinv );
344: }
345: out[0] = tt;
346: i = 1;
347: } else i = 0;
348: /**/
349: if ( osi + nout >= 2 ) {
350: /* out[osi == 0 ? 1 : 0] = (t0 - tn)*Ninv; */
351: tt = tn == (ModNum)0 ? t0 : A_BmodP( t0, tn, P );
352: if ( tt != (ModNum)0 ) {
353: AxBmodP( tt, ModNum, tt, Ninv, P, Pinv );
354: }
355: out[i] = tt;
356: }
357: return;
358: }
359: /****/
360: halfN = n = 1 << (k = log2N - 1); /* == 2^{K-1} */
361: #ifdef POWA_STRIDE
362: halfNa = a1 << k;
363: #endif
364: /*
365: * when k = K-1,
366: * Out[2*i] <= t0 + tn, and Out[2*i+1] <= t0 - tn,
367: * where t0 = In[i] and tn = In[i+2^{K-1}],
368: * for 0 <= i < 2^{K-1}.
369: */
370: qi = in, istep = i1, idist = i1 << k, qo = wk, inwk = 0;
371: for ( i = n; i-- > 0; qi += istep, qo += 2 ) {
372: t0 = qi[0], tn = qi[idist];
373: if ( tn == (ModNum)0 ) qo[0] = qo[1] = t0;
374: else {
375: qo[0] = AplusBmodP( t0, tn, P );
376: qo[1] = A_BmodP( t0, tn, P );
377: }
378: }
379: #ifdef EBUG
380: fprintf( stderr, "::: DFT^{-1} ::: after the first step\n" );
381: for ( qo = wk, i = 0; i < N; i++ ) {
382: if ( (i%5) == 0 ) fprintf( stderr, "\t" );
383: fprintf( stderr, "%10d,", (int)wk[i] );
384: if ( (i%5) == 4 ) fprintf( stderr, "\n" );
385: }
386: if ( (i%5) != 0 ) fprintf( stderr, "\n" );
387: #endif
388: /**/
389: idist = halfN;
390: for ( ostep = 2; --k > 0; inwk = 1 - inwk ) {
391: n >>= 1;
392: istep = ostep; /* == 2^{K-k-1} */
393: ostep <<= 1; /* == 2^{K-k} */
394: if ( inwk ) qi = &wk[N], qo = wk;
395: else qi = wk, qo = &wk[N];
396: qis = qi, qos = qo;
397: /*
398: * for j = 0,
399: * Out[2^{K-k}*i] <= t0 + tn, and Out[2^{K-k}*i + 2^{K-k-1}] <= t0 - tn,
400: * where t0 = In[2^{K-k-1}*i] and tn = In[2^{K-k-1}*i + 2^{K-1}].
401: */
402: for ( i = n, is = os = 0; i-- > 0; qi += istep, qo += ostep ) {
403: t0 = qi[0], tn = qi[idist];
404: if ( tn == (ModNum)0 ) qo[0] = qo[istep] = t0;
405: else {
406: qo[0] = AplusBmodP( t0, tn, P );
407: qo[istep] = A_BmodP( t0, tn, P );
408: }
409: }
410: #ifdef POWA_STRIDE
411: for ( aj = a1, j = 1; j < istep; aj += a1, j++ ) {
412: /*** a = P - powa[halfNa - (aj << k)]; ***/
413: a = powa[halfNa - (aj << k)];
414: #else
415: for ( j = 1; j < istep; j++ ) {
416: /*** a = P - powa[halfN - (j << k)]; ***/
417: a = powa[halfN - (j << k)];
418: #endif
419: qi = ++qis, qo = ++qos;
420: for ( i = n; i-- > 0; qi += istep, qo += ostep ) {
421: t0 = qi[0], tn = qi[idist];
422: if ( tn == (ModNum)0 ) qo[0] = qo[istep] = t0;
423: else {
424: AxBmodP( tn, ModNum, tn, a, P, Pinv );
425: /*** qo[0] = AplusBmodP( t0, tn, P ); ***/
426: /*** qo[istep] = A_BmodP( t0, tn, P ); ***/
427: qo[0] = A_BmodP( t0, tn, P );
428: qo[istep] = AplusBmodP( t0, tn, P );
429: }
430: }
431: }
432: }
433: #if 0
434: final:
435: #endif
436: /*
437: * The final step of k = 0. i can take only the value 0.
438: *
439: * Out[j] <= t0 + \alpha^{-j} * tn;
440: * Out[j + 2^{K-1}] <= t0 - \alpha^{-j} * tn;
441: *
442: * where t0 = In[j] and tn = In[j + 2^{K-1}],
443: * for 0 <= j < 2^{K-1}.
444: *
445: */
446: qo = out, qi = &wk[inwk ? N : 0];
447: if ( (n = osi + nout) > N ) nout = (n = N) - osi; /* true nout */
448: if ( (k = nout - halfN) <= 0 ) { /* no overlap, i.e., no +- */
449: if ( osi <= 0 ) {
450: t0 = qi[0], tn = qi[idist];
451: if ( tn != (ModNum)0 ) t0 = AplusBmodP( t0, tn, P );
452: if ( t0 == (ModNum)0 ) qo[0] = t0;
453: else {
454: AxBmodP( qo[0], ModNum, t0, Ninv, P, Pinv );
455: }
456: #ifdef POWA_STRIDE
457: aj = a1;
458: #endif
459: j = 1, k = n, i = 0, qo += o1;
460: } else {
461: j = osi;
462: #ifdef POWA_STRIDE
463: aj = osi*a1;
464: #endif
465: if ( n <= halfN ) i = 0, k = n; /* lower only */
466: else if ( j == halfN ) goto L_halfN; /* start with j = N/2 */
467: else if ( j > halfN ) goto L_upper; /* upper only */
468: else i = n - halfN + 1, k = halfN; /* we have lower, and i upper */
469: }
470: } else {
471: os = o1*halfN;
472: #ifdef EBUG
473: fprintf( stderr, "::::: DFT^{-1} ::: os=%d*%d=%d\n", o1, halfN, os );
474: #endif
475: if ( osi <= 0 ) {
476: t0 = qi[0], tn = qi[idist];
477: if ( tn == (ModNum)0 ) {
478: if ( t0 == (ModNum)0 ) tt = t0;
479: else { AxBmodP( tt, ModNum, t0, Ninv, P, Pinv ); }
480: qo[0] = qo[os] = tt;
481: } else {
482: tt = AplusBmodP( t0, tn, P );
483: if ( tt == (ModNum)0 ) qo[0] = tt;
484: else { AxBmodP( qo[0], ModNum, tt, Ninv, P, Pinv ); }
485: tt = A_BmodP( t0, tn, P );
486: if ( tt == (ModNum)0 ) qo[os] = tt;
487: else { AxBmodP( qo[os], ModNum, tt, Ninv, P, Pinv ); }
488: }
489: #ifdef EBUG
490: fprintf( stderr, "::::: DFT^{-1} ::: [0]=%d, [%d]=%d\n", (int)qo[0], os, (int)qo[os] );
491: #endif
492: #ifdef POWA_STRIDE
493: j = 1, aj = a1, n = halfN, i = 0, qo += o1, k--;
494: } else j = osi, aj = osi*a1, i = (n -= k) - halfN;
495: /**/
496: for ( ; k-- > 0; aj += a1, j++, qo += o1 ) {
497: #else
498: j = 1, n = halfN, i = 0, qo += o1, k--;
499: } else j = osi, i = (n -= k) - halfN;
500: /**/
501: for ( ; k-- > 0; j++, qo += o1 ) {
502: #endif
503: t0 = qi[j], tn = qi[j+idist];
504: if ( tn == (ModNum)0 ) {
505: if ( t0 != (ModNum)0 ) {
506: AxBmodP( t0, ModNum, t0, Ninv, P, Pinv );
507: }
508: qo[0] = qo[os] = t0;
509: } else {
510: #ifdef POWA_STRIDE
511: /*** a = P - powa[halfNa - aj]; ***/
512: a = powa[halfNa - aj];
513: #else
514: /*** a = P - powa[halfN - j]; ***/
515: a = powa[halfN - j];
516: #endif
517: AxBmodP( tn, ModNum, tn, a, P, Pinv );
518: /*** tt = AplusBmodP( t0, tn, P ); ***/
519: tt = A_BmodP( t0, tn, P );
520: if ( tt == (ModNum)0 ) qo[0] = tt;
521: else { AxBmodP( qo[0], ModNum, tt, Ninv, P, Pinv ); }
522: /*** tt = A_BmodP( t0, tn, P ); ***/
523: tt = AplusBmodP( t0, tn, P );
524: if ( tt == (ModNum)0 ) qo[os] = tt;
525: else { AxBmodP( qo[os], ModNum, tt, Ninv, P, Pinv ); }
526: }
527: }
528: k = halfN;
529: }
530: /*
531: * At this point, k is an upper limit < N/2 for j, i is the number of j's
532: * > N/2, and n is the real upper limit of j.
533: */
534: #ifdef POWA_STRIDE
535: for ( ; j < k; aj += a1, j++, qo += o1 ) {
536: #else
537: for ( ; j < k; j++, qo += o1 ) {
538: #endif
539: t0 = qi[j], tn = qi[j+idist];
540: if ( tn == (ModNum)0 ) {
541: if ( t0 == (ModNum)0 ) qo[0] = t0;
542: else { AxBmodP( qo[0], ModNum, t0, Ninv, P, Pinv ); }
543: } else {
544: #ifdef POWA_STRIDE
545: a = P - powa[halfNa - aj];
546: #else
547: a = P - powa[halfN - j];
548: #endif
549: AxBplusCmodP( tt, ModNum, tn, a, t0, P, Pinv );
550: if ( tt == (ModNum)0 ) qo[0] = tt;
551: else { AxBmodP( qo[0], ModNum, tt, Ninv, P, Pinv ); }
552: }
553: }
554: if ( i <= 0 ) return;
555: /**/
556: L_halfN:
557: t0 = qi[0], tn = qi[idist];
558: tt = tn == (ModNum)0 ? t0 : A_BmodP( t0, tn, P );
559: if ( tt == (ModNum)0 ) qo[0] = tt;
560: else { AxBmodP( qo[0], ModNum, tt, Ninv, P, Pinv ); }
561: qo += o1, j++;
562: /**/
563: L_upper:
564: #ifdef POWA_STRIDE
565: for ( n -= halfN, aj = a1, j = 1; j < n; aj += a1, j++, qo += o1 ) {
566: #else
567: for ( n -= halfN, j = 1; j < n; j++, qo += o1 ) {
568: #endif
569: t0 = qi[j], tn = qi[j+idist];
570: if ( tn == (ModNum)0 ) {
571: if ( t0 == (ModNum)0 ) qo[0] = t0;
572: else { AxBmodP( qo[0], ModNum, t0, Ninv, P, Pinv ); }
573: } else {
574: #ifdef POWA_STRIDE
575: a = powa[halfNa - aj];
576: #else
577: a = powa[halfN - j];
578: #endif
579: AxBplusCmodP( tt, ModNum, tn, a, t0, P, Pinv );
580: if ( tt == (ModNum)0 ) qo[0] = tt;
581: else { AxBmodP( qo[0], ModNum, tt, Ninv, P, Pinv ); }
582: }
583: }
584: }
585:
586: #if defined(USE_FLOAT)
587: void C_PREP_ALPHA( r, log2ord, log2k, n, tbl, P, Pinv )
588: ModNum r, tbl[], P;
589: int log2ord, log2k, n;
590: double Pinv;
591: /*
592: * Let K and k denote log2ord and log2k, respectively.
593: * Let r be a primitive (2^K)-th root of unity in Z/(P), where P is a prime.
594: * Compute a = r^{2^{K-k}}, a primitive (2^k)-th root of unity, and
595: * its powers to the degrees upto n in tbl[*].
596: */
597: {
598: int i;
599: ModNum *q;
600: double t, w, dp = (double) P;
601:
602: for ( w = (double) r, i = log2ord - log2k; i > 0; i-- ) {
603: w *= w;
604: if ( w < dp ) continue;
605: w -= (dp*(double)((int) (w*Pinv)));
606: }
607: tbl[0] = 1;
608: tbl[1] = (ModNum)(t = w);
609: for ( q = &tbl[2], i = n - 1; i > 0; i-- ) {
610: t *= w;
611: if ( t >= dp ) t -= (dp*(double)((int) (t*Pinv)));
612: *q++ = t;
613: }
614: }
615: #else
616:
617: void C_PREP_ALPHA( r, log2ord, log2k, n, tbl, P, Pinv )
618: ModNum r, tbl[], P;
619: int log2ord, log2k, n;
620: double Pinv;
621: /*
622: * Let K and k denote log2ord and log2k, respectively.
623: * Let r be a primitive (2^K)-th root of unity in Z/(P), where P is a prime.
624: * Compute a = r^{2^{K-k}}, a primitive (2^k)-th root of unity, and
625: * its powers to the degrees upto n in tbl[*].
626: */
627: {
628: int i;
629: ModNum *q;
630: ModNum t, w, s;
631:
632: for ( w = r, i = log2ord - log2k; i > 0; i-- ) {
633: AxBmodP( t, ModNum, w, w, P, Pinv );
634: w = t;
635: }
636: tbl[0] = 1;
637: tbl[1] = (ModNum)(t = w);
638: for ( q = &tbl[2], i = n - 1; i > 0; i-- ) {
639: AxBmodP( s, ModNum, t, w, P, Pinv );
640: t = s;
641: *q++ = t;
642: }
643: }
644: #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 / asir2018 / fft