version 1.1, 1999/12/03 07:39:08 |
version 1.5, 2018/03/29 01:32:53 |
|
|
/* $OpenXM: OpenXM/src/asir99/fft/dft.c,v 1.1.1.1 1999/11/10 08:12:27 noro Exp $ */ |
|
/* |
/* |
* store arith modulus |
* Copyright (c) 1994-2000 FUJITSU LABORATORIES LIMITED |
* ------------------------------------------------------- |
* All rights reserved. |
* IEEE: float (23+1) <-> double (52+1) 24 bits |
* |
* int (32) <-> double (52+1) 26 bits |
* FUJITSU LABORATORIES LIMITED ("FLL") hereby grants you a limited, |
* 370: float (24) <-> double (56) 24/25 bits |
* non-exclusive and royalty-free license to use, copy, modify and |
* int (32) <-> double (56) 28 bits |
* redistribute, solely for non-commercial and non-profit purposes, the |
|
* computer program, "Risa/Asir" ("SOFTWARE"), subject to the terms and |
|
* conditions of this Agreement. For the avoidance of doubt, you acquire |
|
* only a limited right to use the SOFTWARE hereunder, and FLL or any |
|
* third party developer retains all rights, including but not limited to |
|
* copyrights, in and to the SOFTWARE. |
|
* |
|
* (1) FLL does not grant you a license in any way for commercial |
|
* purposes. You may use the SOFTWARE only for non-commercial and |
|
* non-profit purposes only, such as academic, research and internal |
|
* business use. |
|
* (2) The SOFTWARE is protected by the Copyright Law of Japan and |
|
* international copyright treaties. If you make copies of the SOFTWARE, |
|
* with or without modification, as permitted hereunder, you shall affix |
|
* to all such copies of the SOFTWARE the above copyright notice. |
|
* (3) An explicit reference to this SOFTWARE and its copyright owner |
|
* shall be made on your publication or presentation in any form of the |
|
* results obtained by use of the SOFTWARE. |
|
* (4) In the event that you modify the SOFTWARE, you shall notify FLL by |
|
* e-mail at risa-admin@sec.flab.fujitsu.co.jp of the detailed specification |
|
* for such modification or the source code of the modified part of the |
|
* SOFTWARE. |
|
* |
|
* THE SOFTWARE IS PROVIDED AS IS WITHOUT ANY WARRANTY OF ANY KIND. FLL |
|
* MAKES ABSOLUTELY NO WARRANTIES, EXPRESSED, IMPLIED OR STATUTORY, AND |
|
* EXPRESSLY DISCLAIMS ANY IMPLIED WARRANTY OF MERCHANTABILITY, FITNESS |
|
* FOR A PARTICULAR PURPOSE OR NONINFRINGEMENT OF THIRD PARTIES' |
|
* RIGHTS. NO FLL DEALER, AGENT, EMPLOYEES IS AUTHORIZED TO MAKE ANY |
|
* MODIFICATIONS, EXTENSIONS, OR ADDITIONS TO THIS WARRANTY. |
|
* UNDER NO CIRCUMSTANCES AND UNDER NO LEGAL THEORY, TORT, CONTRACT, |
|
* OR OTHERWISE, SHALL FLL BE LIABLE TO YOU OR ANY OTHER PERSON FOR ANY |
|
* DIRECT, INDIRECT, SPECIAL, INCIDENTAL, PUNITIVE OR CONSEQUENTIAL |
|
* DAMAGES OF ANY CHARACTER, INCLUDING, WITHOUT LIMITATION, DAMAGES |
|
* ARISING OUT OF OR RELATING TO THE SOFTWARE OR THIS AGREEMENT, DAMAGES |
|
* FOR LOSS OF GOODWILL, WORK STOPPAGE, OR LOSS OF DATA, OR FOR ANY |
|
* DAMAGES, EVEN IF FLL SHALL HAVE BEEN INFORMED OF THE POSSIBILITY OF |
|
* SUCH DAMAGES, OR FOR ANY CLAIM BY ANY OTHER PARTY. EVEN IF A PART |
|
* OF THE SOFTWARE HAS BEEN DEVELOPED BY A THIRD PARTY, THE THIRD PARTY |
|
* DEVELOPER SHALL HAVE NO LIABILITY IN CONNECTION WITH THE USE, |
|
* PERFORMANCE OR NON-PERFORMANCE OF THE SOFTWARE. |
|
* |
|
* $OpenXM: OpenXM_contrib2/asir2000/fft/dft.c,v 1.4 2003/02/14 22:29:10 ohara Exp $ |
|
*/ |
|
/* |
|
* store arith modulus |
|
* ------------------------------------------------------- |
|
* IEEE: float (23+1) <-> double (52+1) 24 bits |
|
* int (32) <-> double (52+1) 26 bits |
|
* 370: float (24) <-> double (56) 24/25 bits |
|
* int (32) <-> double (56) 28 bits |
*/ |
*/ |
|
|
|
|
|
|
|
|
void C_DFT_FORE( in, nin, i1, K, powa, |
void C_DFT_FORE( in, nin, i1, K, powa, |
#ifdef POWA_STRIDE |
#ifdef POWA_STRIDE |
a1, |
a1, |
#endif |
#endif |
out, o1, P, Pinv, wk ) |
out, o1, P, Pinv, wk ) |
ModNum in[], powa[], out[], wk[]; |
ModNum in[], powa[], out[], wk[]; |
ModNum P; |
ModNum P; |
int nin, i1, K, o1; |
int nin, i1, K, o1; |
|
|
ModNum *qi, *qo, *qis, *qos, a, t0, tn; |
ModNum *qi, *qo, *qis, *qos, a, t0, tn; |
|
|
/* |
/* |
* for k = K-1, ..., 0 do the following: |
* for k = K-1, ..., 0 do the following: |
* |
* |
* Out[2^{K-k}*i + j] <= t0 + \alpha^{2^k * j} * tn; |
* Out[2^{K-k}*i + j] <= t0 + \alpha^{2^k * j} * tn; |
* Out[2^{K-k}*i + j + 2^{K-k-1}] <= t0 - \alpha^{2^k * j} * tn; |
* Out[2^{K-k}*i + j + 2^{K-k-1}] <= t0 - \alpha^{2^k * j} * tn; |
* |
* |
* where t0 = In[2^{K-k-1}*i + j] and |
* where t0 = In[2^{K-k-1}*i + j] and |
* tn = In[2^{K-k-1}*i + j + 2^{K-1}], |
* tn = In[2^{K-k-1}*i + j + 2^{K-1}], |
* for 0 <= i < 2^k and 0 <= j < 2^{K-k-1}. |
* for 0 <= i < 2^k and 0 <= j < 2^{K-k-1}. |
* |
* |
*/ |
*/ |
|
|
for ( i = nin; i > 0; i--, qi += i1, qo += ostep ) qo[0] = qo[odist] = qi[0]; |
for ( i = nin; i > 0; i--, qi += i1, qo += ostep ) qo[0] = qo[odist] = qi[0]; |
for ( i = n - nin; i > 0; i--, qo += ostep ) qo[0] = qo[odist] = 0; |
for ( i = n - nin; i > 0; i--, qo += ostep ) qo[0] = qo[odist] = 0; |
} else { |
} else { |
idist = i1 << k; /* distance = 2^{K-1} */ |
idist = i1 << k; /* distance = 2^{K-1} */ |
for ( j = i = nin - n; i > 0; i--, qi += i1, qo += ostep ) { |
for ( j = i = nin - n; i > 0; i--, qi += i1, qo += ostep ) { |
t0 = qi[0], tn = qi[idist]; |
t0 = qi[0], tn = qi[idist]; |
if ( tn == (ModNum)0 ) qo[0] = qo[odist] = t0; |
if ( tn == (ModNum)0 ) qo[0] = qo[odist] = t0; |
else { |
else { |
qo[0] = AplusBmodP( t0, tn, P ); |
qo[0] = AplusBmodP( t0, tn, P ); |
qo[odist] = A_BmodP( t0, tn, P ); |
qo[odist] = A_BmodP( t0, tn, P ); |
} |
} |
} |
} |
for ( i = n - j; i > 0; i--, qi += i1, qo += ostep ) |
for ( i = n - j; i > 0; i--, qi += i1, qo += ostep ) |
|
|
#ifndef OPT |
#ifndef OPT |
K_k = (K_k_1 = K_k) + 1; |
K_k = (K_k_1 = K_k) + 1; |
#endif |
#endif |
n >>= 1; /* == 2^k */ |
n >>= 1; /* == 2^k */ |
Ndiv2n <<= 1; /* == 2^{K-k-1} */ |
Ndiv2n <<= 1; /* == 2^{K-k-1} */ |
if ( inwk ) |
if ( inwk ) |
qi = wk, qo = out, iostep = 1, oostep = o1, idist = idistwk; |
qi = wk, qo = out, iostep = 1, oostep = o1, idist = idistwk; |
else |
else |
qi = out, qo = wk, iostep = o1, oostep = 1, idist = idistout; |
qi = out, qo = wk, iostep = o1, oostep = 1, idist = idistout; |
qis = qi, qos = qo; |
qis = qi, qos = qo; |
/**/ |
/**/ |
istep = ostep; /* = iostep << K_k_1; */ |
istep = ostep; /* = iostep << K_k_1; */ |
#ifndef OPT |
#ifndef OPT |
ostep = oostep << K_k; |
ostep = oostep << K_k; |
odist = oostep << K_k_1; |
odist = oostep << K_k_1; |
|
|
t0 = qi[0], tn = qi[idist]; |
t0 = qi[0], tn = qi[idist]; |
if ( tn == (ModNum)0 ) qo[0] = qo[odist] = t0; |
if ( tn == (ModNum)0 ) qo[0] = qo[odist] = t0; |
else { |
else { |
qo[0] = AplusBmodP( t0, tn, P ); |
qo[0] = AplusBmodP( t0, tn, P ); |
qo[odist] = A_BmodP( t0, tn, P ); |
qo[odist] = A_BmodP( t0, tn, P ); |
} |
} |
} |
} |
/**/ |
/**/ |
|
|
#endif |
#endif |
qi = (qis += iostep), qo = (qos += oostep); |
qi = (qis += iostep), qo = (qos += oostep); |
for ( i = 0; i < n; i++, qi += istep, qo += ostep ) { |
for ( i = 0; i < n; i++, qi += istep, qo += ostep ) { |
t0 = qi[0], tn = qi[idist]; |
t0 = qi[0], tn = qi[idist]; |
if ( tn == (ModNum)0 ) qo[0] = qo[odist] = t0; |
if ( tn == (ModNum)0 ) qo[0] = qo[odist] = t0; |
else { |
else { |
AxBmodP( tn, ModNum, tn, a, P, Pinv ); |
AxBmodP( tn, ModNum, tn, a, P, Pinv ); |
qo[0] = AplusBmodP( t0, tn, P ); |
qo[0] = AplusBmodP( t0, tn, P ); |
qo[odist] = A_BmodP( t0, tn, P ); |
qo[odist] = A_BmodP( t0, tn, P ); |
} |
} |
} |
} |
} |
} |
} |
} |
/* |
/* |
* When k = 0, for i = 0 (0 <= i < 2^0=1) and 0 <= j < 2^{K-1}, |
* When k = 0, for i = 0 (0 <= i < 2^0=1) and 0 <= j < 2^{K-1}, |
* Out[j] <= t0 + \alpha^{j}*tn and Out[j+2^{K-1}] <= t0 - \alpha^{j}*tn, |
* Out[j] <= t0 + \alpha^{j}*tn and Out[j+2^{K-1}] <= t0 - \alpha^{j}*tn, |
* where t0 = In[j] and tn = In[j + 2^{K-1}]. |
* where t0 = In[j] and tn = In[j + 2^{K-1}]. |
*/ |
*/ |
qi = wk, qo = out, idist = idistwk, odist = idistout; |
qi = wk, qo = out, idist = idistwk, odist = idistout; |
|
|
#ifndef OPT |
#ifndef OPT |
K_k = (K_k_1 = K_k) + 1; |
K_k = (K_k_1 = K_k) + 1; |
#endif |
#endif |
n >>= 1; /* == 2^k */ |
n >>= 1; /* == 2^k */ |
Ndiv2n <<= 1; /* == 2^{K-k-1} */ |
Ndiv2n <<= 1; /* == 2^{K-k-1} */ |
if ( inwk ) qi = wk, qo = out, istep = 1, ostep = o1, idist = idistwk; |
if ( inwk ) qi = wk, qo = out, istep = 1, ostep = o1, idist = idistwk; |
else qi = out, qo = wk, istep = o1, ostep = 1, idist = idistout; |
else qi = out, qo = wk, istep = o1, ostep = 1, idist = idistout; |
qis = qi, qos = qo; |
qis = qi, qos = qo; |
/**/ |
/**/ |
iostep = oostep; /* = istep << K_k_1; */ |
iostep = oostep; /* = istep << K_k_1; */ |
#ifndef OPT |
#ifndef OPT |
oostep = ostep << K_k; |
oostep = ostep << K_k; |
odist = ostep << K_k_1; |
odist = ostep << K_k_1; |
|
|
t0 = qis[0], tn = qis[idist]; |
t0 = qis[0], tn = qis[idist]; |
if ( tn == (ModNum)0 ) qos[0] = qos[odist] = t0; |
if ( tn == (ModNum)0 ) qos[0] = qos[odist] = t0; |
else { |
else { |
qos[0] = AplusBmodP( t0, tn, P ); |
qos[0] = AplusBmodP( t0, tn, P ); |
qos[odist] = A_BmodP( t0, tn, P ); |
qos[odist] = A_BmodP( t0, tn, P ); |
} |
} |
#ifdef POWA_STRIDE |
#ifdef POWA_STRIDE |
for ( aj = a1, j = 1, qi = qis, qo = qos; j < Ndiv2n; aj += a1, j++ ) { |
for ( aj = a1, j = 1, qi = qis, qo = qos; j < Ndiv2n; aj += a1, j++ ) { |
#else |
#else |
for ( j = 1, qi = qis, qo = qos; j < Ndiv2n; j++ ) { |
for ( j = 1, qi = qis, qo = qos; j < Ndiv2n; j++ ) { |
#endif |
#endif |
qi += istep, qo += ostep; |
qi += istep, qo += ostep; |
t0 = qi[0], tn = qi[idist]; |
t0 = qi[0], tn = qi[idist]; |
if ( tn == (ModNum)0 ) qo[0] = qo[odist] = t0; |
if ( tn == (ModNum)0 ) qo[0] = qo[odist] = t0; |
else { |
else { |
#ifdef POWA_STRIDE |
#ifdef POWA_STRIDE |
a = powa[aj << k]; |
a = powa[aj << k]; |
#else |
#else |
a = powa[j << k]; |
a = powa[j << k]; |
#endif |
#endif |
AxBmodP( tn, ModNum, tn, a, P, Pinv ); |
AxBmodP( tn, ModNum, tn, a, P, Pinv ); |
qo[0] = AplusBmodP( t0, tn, P ); |
qo[0] = AplusBmodP( t0, tn, P ); |
qo[odist] = A_BmodP( t0, tn, P ); |
qo[odist] = A_BmodP( t0, tn, P ); |
} |
} |
} |
} |
} |
} |
} |
} |
|
|
* NOTE: |
* NOTE: |
* (1) Let \alpha be a primitive N-th root of unity modulo P, |
* (1) Let \alpha be a primitive N-th root of unity modulo P, |
* where N is even, |
* where N is even, |
* and t be an integer s.t. 0 <= t < N/2. |
* and t be an integer s.t. 0 <= t < N/2. |
* Then, 1/\alpha^t is given by (- \alpha^s) with s = N/2 - t. |
* Then, 1/\alpha^t is given by (- \alpha^s) with s = N/2 - t. |
* |
* |
* (2) Let P be a prime s.t. P = 2^n*Q + 1, where Q is an odd integer. |
* (2) Let P be a prime s.t. P = 2^n*Q + 1, where Q is an odd integer. |
* Then, for 0 < s <= n, 2^{-s} \equiv - 2^{n-s}*Q \bmod P, |
* Then, for 0 < s <= n, 2^{-s} \equiv - 2^{n-s}*Q \bmod P, |
* because 1 \equiv - 2^n*Q \bmod P. |
* because 1 \equiv - 2^n*Q \bmod P. |
* |
* |
**********************************************************************/ |
**********************************************************************/ |
|
|
void C_DFT_BACK( in, N, i1, log2N, powa, |
void C_DFT_BACK( in, N, i1, log2N, powa, |
#ifdef POWA_STRIDE |
#ifdef POWA_STRIDE |
a1, |
a1, |
#endif |
#endif |
out, o1, osi, nout, Ninv, P, Pinv, wk ) |
out, o1, osi, nout, Ninv, P, Pinv, wk ) |
ModNum in[], powa[], out[], wk[]; |
ModNum in[], powa[], out[], wk[]; |
int N, log2N, osi, nout, i1, o1; |
int N, log2N, osi, nout, i1, o1; |
#ifdef POWA_STRIDE |
#ifdef POWA_STRIDE |
|
|
ModNum *qi, *qo, *qis, *qos, a, t0, tn, tt; |
ModNum *qi, *qo, *qis, *qos, a, t0, tn, tt; |
|
|
/* |
/* |
* for k = K-1, ..., 0 do the following: |
* for k = K-1, ..., 0 do the following: |
* |
* |
* Out[2^{K-k}*i + j] <= t0 + \alpha^{- 2^k * j} * tn; |
* Out[2^{K-k}*i + j] <= t0 + \alpha^{- 2^k * j} * tn; |
* = t0 - \alpha^{N/2 - 2^k * j}*tn; |
* = t0 - \alpha^{N/2 - 2^k * j}*tn; |
* Out[2^{K-k}*i + j + 2^{K-k-1}] <= t0 - \alpha^{- 2^k * j} * tn; |
* Out[2^{K-k}*i + j + 2^{K-k-1}] <= t0 - \alpha^{- 2^k * j} * tn; |
* = t0 + \alpha^{N/2 - 2^k * j}*tn; |
* = t0 + \alpha^{N/2 - 2^k * j}*tn; |
* |
* |
* where t0 = In[2^{K-k-1}*i + j] and |
* where t0 = In[2^{K-k-1}*i + j] and |
* tn = In[2^{K-k-1}*i + j + 2^{K-1}], |
* tn = In[2^{K-k-1}*i + j + 2^{K-1}], |
* for 0 <= i < 2^k and 0 <= j < 2^{K-k-1}. |
* for 0 <= i < 2^k and 0 <= j < 2^{K-k-1}. |
* |
* |
*/ |
*/ |
if ( log2N == 1 ) { |
if ( log2N == 1 ) { |
/* K = 1, k = 0, 0 <= i < 1, 0 <= j < 1. |
/* K = 1, k = 0, 0 <= i < 1, 0 <= j < 1. |
* Out[0] <= t0 + tn, Out[1] <= t0 - tn, |
* Out[0] <= t0 + tn, Out[1] <= t0 - tn, |
* where t0 = In[0] and Tn = In[1]. |
* where t0 = In[0] and Tn = In[1]. |
*/ |
*/ |
t0 = in[0], tn = in[i1]; |
t0 = in[0], tn = in[i1]; |
|
|
/* out[0] = (t0 + tn)*Ninv; */ |
/* out[0] = (t0 + tn)*Ninv; */ |
tt = tn == (ModNum)0 ? t0 : AplusBmodP( t0, tn, P ); |
tt = tn == (ModNum)0 ? t0 : AplusBmodP( t0, tn, P ); |
if ( tt != (ModNum)0 ) { |
if ( tt != (ModNum)0 ) { |
AxBmodP( tt, ModNum, tt, Ninv, P, Pinv ); |
AxBmodP( tt, ModNum, tt, Ninv, P, Pinv ); |
} |
} |
out[0] = tt; |
out[0] = tt; |
i = 1; |
i = 1; |
|
|
/* out[osi == 0 ? 1 : 0] = (t0 - tn)*Ninv; */ |
/* out[osi == 0 ? 1 : 0] = (t0 - tn)*Ninv; */ |
tt = tn == (ModNum)0 ? t0 : A_BmodP( t0, tn, P ); |
tt = tn == (ModNum)0 ? t0 : A_BmodP( t0, tn, P ); |
if ( tt != (ModNum)0 ) { |
if ( tt != (ModNum)0 ) { |
AxBmodP( tt, ModNum, tt, Ninv, P, Pinv ); |
AxBmodP( tt, ModNum, tt, Ninv, P, Pinv ); |
} |
} |
out[i] = tt; |
out[i] = tt; |
} |
} |
return; |
return; |
} |
} |
/****/ |
/****/ |
halfN = n = 1 << (k = log2N - 1); /* == 2^{K-1} */ |
halfN = n = 1 << (k = log2N - 1); /* == 2^{K-1} */ |
#ifdef POWA_STRIDE |
#ifdef POWA_STRIDE |
halfNa = a1 << k; |
halfNa = a1 << k; |
#endif |
#endif |
/* |
/* |
* when k = K-1, |
* when k = K-1, |
* Out[2*i] <= t0 + tn, and Out[2*i+1] <= t0 - tn, |
* Out[2*i] <= t0 + tn, and Out[2*i+1] <= t0 - tn, |
* where t0 = In[i] and tn = In[i+2^{K-1}], |
* where t0 = In[i] and tn = In[i+2^{K-1}], |
* for 0 <= i < 2^{K-1}. |
* for 0 <= i < 2^{K-1}. |
*/ |
*/ |
|
|
idist = halfN; |
idist = halfN; |
for ( ostep = 2; --k > 0; inwk = 1 - inwk ) { |
for ( ostep = 2; --k > 0; inwk = 1 - inwk ) { |
n >>= 1; |
n >>= 1; |
istep = ostep; /* == 2^{K-k-1} */ |
istep = ostep; /* == 2^{K-k-1} */ |
ostep <<= 1; /* == 2^{K-k} */ |
ostep <<= 1; /* == 2^{K-k} */ |
if ( inwk ) qi = &wk[N], qo = wk; |
if ( inwk ) qi = &wk[N], qo = wk; |
else qi = wk, qo = &wk[N]; |
else qi = wk, qo = &wk[N]; |
qis = qi, qos = qo; |
qis = qi, qos = qo; |
/* |
/* |
* for j = 0, |
* for j = 0, |
* Out[2^{K-k}*i] <= t0 + tn, and Out[2^{K-k}*i + 2^{K-k-1}] <= t0 - tn, |
* Out[2^{K-k}*i] <= t0 + tn, and Out[2^{K-k}*i + 2^{K-k-1}] <= t0 - tn, |
* where t0 = In[2^{K-k-1}*i] and tn = In[2^{K-k-1}*i + 2^{K-1}]. |
* where t0 = In[2^{K-k-1}*i] and tn = In[2^{K-k-1}*i + 2^{K-1}]. |
*/ |
*/ |
for ( i = n, is = os = 0; i-- > 0; qi += istep, qo += ostep ) { |
for ( i = n, is = os = 0; i-- > 0; qi += istep, qo += ostep ) { |
t0 = qi[0], tn = qi[idist]; |
t0 = qi[0], tn = qi[idist]; |
if ( tn == (ModNum)0 ) qo[0] = qo[istep] = t0; |
if ( tn == (ModNum)0 ) qo[0] = qo[istep] = t0; |
else { |
else { |
qo[0] = AplusBmodP( t0, tn, P ); |
qo[0] = AplusBmodP( t0, tn, P ); |
qo[istep] = A_BmodP( t0, tn, P ); |
qo[istep] = A_BmodP( t0, tn, P ); |
} |
} |
} |
} |
#ifdef POWA_STRIDE |
#ifdef POWA_STRIDE |
|
|
#endif |
#endif |
qi = ++qis, qo = ++qos; |
qi = ++qis, qo = ++qos; |
for ( i = n; i-- > 0; qi += istep, qo += ostep ) { |
for ( i = n; i-- > 0; qi += istep, qo += ostep ) { |
t0 = qi[0], tn = qi[idist]; |
t0 = qi[0], tn = qi[idist]; |
if ( tn == (ModNum)0 ) qo[0] = qo[istep] = t0; |
if ( tn == (ModNum)0 ) qo[0] = qo[istep] = t0; |
else { |
else { |
AxBmodP( tn, ModNum, tn, a, P, Pinv ); |
AxBmodP( tn, ModNum, tn, a, P, Pinv ); |
/*** qo[0] = AplusBmodP( t0, tn, P ); ***/ |
/*** qo[0] = AplusBmodP( t0, tn, P ); ***/ |
/*** qo[istep] = A_BmodP( t0, tn, P ); ***/ |
/*** qo[istep] = A_BmodP( t0, tn, P ); ***/ |
qo[0] = A_BmodP( t0, tn, P ); |
qo[0] = A_BmodP( t0, tn, P ); |
qo[istep] = AplusBmodP( t0, tn, P ); |
qo[istep] = AplusBmodP( t0, tn, P ); |
} |
} |
} |
} |
} |
} |
} |
} |
|
|
/* |
/* |
* The final step of k = 0. i can take only the value 0. |
* The final step of k = 0. i can take only the value 0. |
* |
* |
* Out[j] <= t0 + \alpha^{-j} * tn; |
* Out[j] <= t0 + \alpha^{-j} * tn; |
* Out[j + 2^{K-1}] <= t0 - \alpha^{-j} * tn; |
* Out[j + 2^{K-1}] <= t0 - \alpha^{-j} * tn; |
* |
* |
* where t0 = In[j] and tn = In[j + 2^{K-1}], |
* where t0 = In[j] and tn = In[j + 2^{K-1}], |
* for 0 <= j < 2^{K-1}. |
* for 0 <= j < 2^{K-1}. |
* |
* |
*/ |
*/ |
qo = out, qi = &wk[inwk ? N : 0]; |
qo = out, qi = &wk[inwk ? N : 0]; |
if ( (n = osi + nout) > N ) nout = (n = N) - osi; /* true nout */ |
if ( (n = osi + nout) > N ) nout = (n = N) - osi; /* true nout */ |
if ( (k = nout - halfN) <= 0 ) { /* no overlap, i.e., no +- */ |
if ( (k = nout - halfN) <= 0 ) { /* no overlap, i.e., no +- */ |
if ( osi <= 0 ) { |
if ( osi <= 0 ) { |
t0 = qi[0], tn = qi[idist]; |
t0 = qi[0], tn = qi[idist]; |
if ( tn != (ModNum)0 ) t0 = AplusBmodP( t0, tn, P ); |
if ( tn != (ModNum)0 ) t0 = AplusBmodP( t0, tn, P ); |
if ( t0 == (ModNum)0 ) qo[0] = t0; |
if ( t0 == (ModNum)0 ) qo[0] = t0; |
else { |
else { |
AxBmodP( qo[0], ModNum, t0, Ninv, P, Pinv ); |
AxBmodP( qo[0], ModNum, t0, Ninv, P, Pinv ); |
} |
} |
#ifdef POWA_STRIDE |
#ifdef POWA_STRIDE |
aj = a1; |
aj = a1; |
|
|
#ifdef POWA_STRIDE |
#ifdef POWA_STRIDE |
aj = osi*a1; |
aj = osi*a1; |
#endif |
#endif |
if ( n <= halfN ) i = 0, k = n; /* lower only */ |
if ( n <= halfN ) i = 0, k = n; /* lower only */ |
else if ( j == halfN ) goto L_halfN; /* start with j = N/2 */ |
else if ( j == halfN ) goto L_halfN; /* start with j = N/2 */ |
else if ( j > halfN ) goto L_upper; /* upper only */ |
else if ( j > halfN ) goto L_upper; /* upper only */ |
else i = n - halfN + 1, k = halfN; /* we have lower, and i upper */ |
else i = n - halfN + 1, k = halfN; /* we have lower, and i upper */ |
} |
} |
} else { |
} else { |
os = o1*halfN; |
os = o1*halfN; |
|
|
if ( osi <= 0 ) { |
if ( osi <= 0 ) { |
t0 = qi[0], tn = qi[idist]; |
t0 = qi[0], tn = qi[idist]; |
if ( tn == (ModNum)0 ) { |
if ( tn == (ModNum)0 ) { |
if ( t0 == (ModNum)0 ) tt = t0; |
if ( t0 == (ModNum)0 ) tt = t0; |
else { AxBmodP( tt, ModNum, t0, Ninv, P, Pinv ); } |
else { AxBmodP( tt, ModNum, t0, Ninv, P, Pinv ); } |
qo[0] = qo[os] = tt; |
qo[0] = qo[os] = tt; |
} else { |
} else { |
tt = AplusBmodP( t0, tn, P ); |
tt = AplusBmodP( t0, tn, P ); |
if ( tt == (ModNum)0 ) qo[0] = tt; |
if ( tt == (ModNum)0 ) qo[0] = tt; |
else { AxBmodP( qo[0], ModNum, tt, Ninv, P, Pinv ); } |
else { AxBmodP( qo[0], ModNum, tt, Ninv, P, Pinv ); } |
tt = A_BmodP( t0, tn, P ); |
tt = A_BmodP( t0, tn, P ); |
if ( tt == (ModNum)0 ) qo[os] = tt; |
if ( tt == (ModNum)0 ) qo[os] = tt; |
else { AxBmodP( qo[os], ModNum, tt, Ninv, P, Pinv ); } |
else { AxBmodP( qo[os], ModNum, tt, Ninv, P, Pinv ); } |
} |
} |
#ifdef EBUG |
#ifdef EBUG |
fprintf( stderr, "::::: DFT^{-1} ::: [0]=%d, [%d]=%d\n", (int)qo[0], os, (int)qo[os] ); |
fprintf( stderr, "::::: DFT^{-1} ::: [0]=%d, [%d]=%d\n", (int)qo[0], os, (int)qo[os] ); |
|
|
#endif |
#endif |
t0 = qi[j], tn = qi[j+idist]; |
t0 = qi[j], tn = qi[j+idist]; |
if ( tn == (ModNum)0 ) { |
if ( tn == (ModNum)0 ) { |
if ( t0 != (ModNum)0 ) { |
if ( t0 != (ModNum)0 ) { |
AxBmodP( t0, ModNum, t0, Ninv, P, Pinv ); |
AxBmodP( t0, ModNum, t0, Ninv, P, Pinv ); |
} |
} |
qo[0] = qo[os] = t0; |
qo[0] = qo[os] = t0; |
} else { |
} else { |
#ifdef POWA_STRIDE |
#ifdef POWA_STRIDE |
/*** a = P - powa[halfNa - aj]; ***/ |
/*** a = P - powa[halfNa - aj]; ***/ |
a = powa[halfNa - aj]; |
a = powa[halfNa - aj]; |
#else |
#else |
/*** a = P - powa[halfN - j]; ***/ |
/*** a = P - powa[halfN - j]; ***/ |
a = powa[halfN - j]; |
a = powa[halfN - j]; |
#endif |
#endif |
AxBmodP( tn, ModNum, tn, a, P, Pinv ); |
AxBmodP( tn, ModNum, tn, a, P, Pinv ); |
/*** tt = AplusBmodP( t0, tn, P ); ***/ |
/*** tt = AplusBmodP( t0, tn, P ); ***/ |
tt = A_BmodP( t0, tn, P ); |
tt = A_BmodP( t0, tn, P ); |
if ( tt == (ModNum)0 ) qo[0] = tt; |
if ( tt == (ModNum)0 ) qo[0] = tt; |
else { AxBmodP( qo[0], ModNum, tt, Ninv, P, Pinv ); } |
else { AxBmodP( qo[0], ModNum, tt, Ninv, P, Pinv ); } |
/*** tt = A_BmodP( t0, tn, P ); ***/ |
/*** tt = A_BmodP( t0, tn, P ); ***/ |
tt = AplusBmodP( t0, tn, P ); |
tt = AplusBmodP( t0, tn, P ); |
if ( tt == (ModNum)0 ) qo[os] = tt; |
if ( tt == (ModNum)0 ) qo[os] = tt; |
else { AxBmodP( qo[os], ModNum, tt, Ninv, P, Pinv ); } |
else { AxBmodP( qo[os], ModNum, tt, Ninv, P, Pinv ); } |
} |
} |
} |
} |
k = halfN; |
k = halfN; |
|
|
} |
} |
} |
} |
|
|
#if USE_FLOAT |
#if defined(USE_FLOAT) |
void C_PREP_ALPHA( r, log2ord, log2k, n, tbl, P, Pinv ) |
void C_PREP_ALPHA( r, log2ord, log2k, n, tbl, P, Pinv ) |
ModNum r, tbl[], P; |
ModNum r, tbl[], P; |
int log2ord, log2k, n; |
int log2ord, log2k, n; |