Annotation of OpenXM_contrib2/asir2000/fft/dft.c, Revision 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