/*
* Copyright (c) 1994-2000 FUJITSU LABORATORIES LIMITED
* All rights reserved.
*
* FUJITSU LABORATORIES LIMITED ("FLL") hereby grants you a limited,
* non-exclusive and royalty-free license to use, copy, modify and
* 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/asir2018/fft/dft.c,v 1.2 2020/10/04 03:14:09 noro 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
*/
#include "dft.h"
/*
#define C_DFT_FORE C_DFTfore
#define C_DFT_BACK C_DFTback
#define C_PREP_ALPHA C_prep_alpha
*/
/*****************************************************/
void C_DFT_FORE( in, nin, i1, K, powa,
#ifdef POWA_STRIDE
a1,
#endif
out, o1, P, Pinv, wk )
ModNum in[], powa[], out[], wk[];
ModNum P;
int nin, i1, K, o1;
#ifdef POWA_STRIDE
int a1;
#endif
double Pinv;
/*
* Let N = 2^K, and \alpha be a primitive N-th root of 1.
* out[0:N-1][o1] = DFT of in[0:nin-1][i1] using \alpha modulo P.
* The powers \alpha (to the degrees upto N/2=2^{K-1} at least) are given in powa[*],
* (or powa[a1**]).
* nin <= N assumed.
* wk[] must have 2^K entries at least.
*/
{
int istep, ostep, i, j, k, idist, odist, n, Ndiv2n, K_k,
#ifndef OPT
K_k_1,
#endif
#ifdef POWA_STRIDE
aj,
#endif
inwk;
int idistwk, idistout, iostep, oostep;
ModNum *qi, *qo, *qis, *qos, a, t0, tn;
/*
* 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 + 2^{K-k-1}] <= t0 - \alpha^{2^k * j} * tn;
*
* where t0 = In[2^{K-k-1}*i + j] and
* tn = In[2^{K-k-1}*i + j + 2^{K-1}],
* for 0 <= i < 2^k and 0 <= j < 2^{K-k-1}.
*
*/
/* In the following, n = 2^k and Ndiv2n = 2^{K-k-1}.
*/
n = 1 << (k = K - 1);
idistwk = n, idistout = o1 << k;
/* when k = K-1, Out[2*i] <= In[i] + In[i+2^{K-1}],
* and Out[2*i+1] <= In[i] - In[i+2^{K-1}], for 0 <= i < 2^{K-1}.
*/
if ( k&1 ) qo = wk, ostep = 2, inwk = 1, odist = 1;
else qo = out, ostep = o1 << 1, inwk = 0, odist = o1;
qi = in;
/**/
if ( nin <= n ) {
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;
} else {
idist = i1 << k; /* distance = 2^{K-1} */
for ( j = i = nin - n; i > 0; i--, qi += i1, qo += ostep ) {
t0 = qi[0], tn = qi[idist];
if ( tn == (ModNum)0 ) qo[0] = qo[odist] = t0;
else {
qo[0] = AplusBmodP( t0, tn, P );
qo[odist] = A_BmodP( t0, tn, P );
}
}
for ( i = n - j; i > 0; i--, qi += i1, qo += ostep )
qo[0] = qo[odist] = qi[0];
}
/******/
#ifndef J_INNER
for ( K_k = 1, Ndiv2n = 1; --k > 0; inwk = 1 - inwk ) {
#ifndef OPT
K_k = (K_k_1 = K_k) + 1;
#endif
n >>= 1; /* == 2^k */
Ndiv2n <<= 1; /* == 2^{K-k-1} */
if ( inwk )
qi = wk, qo = out, iostep = 1, oostep = o1, idist = idistwk;
else
qi = out, qo = wk, iostep = o1, oostep = 1, idist = idistout;
qis = qi, qos = qo;
/**/
istep = ostep; /* = iostep << K_k_1; */
#ifndef OPT
ostep = oostep << K_k;
odist = oostep << K_k_1;
#else
odist = oostep << K_k;
++K_k;
ostep = oostep << K_k;
#endif /* OPT */
/* j = 0 */
for ( i = n; i-- > 0; qi += istep, qo += ostep ) {
t0 = qi[0], tn = qi[idist];
if ( tn == (ModNum)0 ) qo[0] = qo[odist] = t0;
else {
qo[0] = AplusBmodP( t0, tn, P );
qo[odist] = A_BmodP( t0, tn, P );
}
}
/**/
#ifdef POWA_STRIDE
for ( aj = a1, j = 1; j < Ndiv2n; aj += a1, j++ ) {
a = powa[aj << k];
#else
for ( j = 1; j < Ndiv2n; j++ ) {
a = powa[j << k];
#endif
qi = (qis += iostep), qo = (qos += oostep);
for ( i = 0; i < n; i++, qi += istep, qo += ostep ) {
t0 = qi[0], tn = qi[idist];
if ( tn == (ModNum)0 ) qo[0] = qo[odist] = t0;
else {
AxBmodP( tn, ModNum, tn, a, P, Pinv );
qo[0] = AplusBmodP( 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},
* 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}].
*/
qi = wk, qo = out, idist = idistwk, odist = idistout;
/* == 2^{K-1}, == o1 << (K-1) */
t0 = qi[0], tn = qi[idist];
if ( tn == (ModNum)0 ) qo[0] = qo[odist] = t0;
else {
qo[0] = AplusBmodP( t0, tn, P );
qo[odist] = A_BmodP( t0, tn, P );
}
#ifdef POWA_STRIDE
for ( k = o1, aj = a1, j = 1; j < idist; aj += a1, j++ ) {
#else
for ( k = o1, j = 1; j < idist; j++ ) {
#endif
qi++, qo += k;
t0 = qi[0], tn = qi[idist];
if ( tn == (ModNum)0 ) qo[0] = qo[odist] = t0;
else {
#ifdef POWA_STRIDE
a = powa[aj];
#else
a = powa[j];
#endif
AxBmodP( tn, ModNum, tn, a, P, Pinv );
qo[0] = AplusBmodP( t0, tn, P );
qo[odist] = A_BmodP( t0, tn, P );
}
}
#else
for ( K_k = 1, Ndiv2n = 1; --k >= 0; inwk = 1 - inwk ) {
#ifndef OPT
K_k = (K_k_1 = K_k) + 1;
#endif
n >>= 1; /* == 2^k */
Ndiv2n <<= 1; /* == 2^{K-k-1} */
if ( inwk ) qi = wk, qo = out, istep = 1, ostep = o1, idist = idistwk;
else qi = out, qo = wk, istep = o1, ostep = 1, idist = idistout;
qis = qi, qos = qo;
/**/
iostep = oostep; /* = istep << K_k_1; */
#ifndef OPT
oostep = ostep << K_k;
odist = ostep << K_k_1;
#else
odist = ostep << K_k;
++K_k;
oostep = ostep << K_k;
#endif /* OPT */
for ( i = n; i-- > 0; qis += iostep, qos += oostep ) {
/* j = 0 */
t0 = qis[0], tn = qis[idist];
if ( tn == (ModNum)0 ) qos[0] = qos[odist] = t0;
else {
qos[0] = AplusBmodP( t0, tn, P );
qos[odist] = A_BmodP( t0, tn, P );
}
#ifdef POWA_STRIDE
for ( aj = a1, j = 1, qi = qis, qo = qos; j < Ndiv2n; aj += a1, j++ ) {
#else
for ( j = 1, qi = qis, qo = qos; j < Ndiv2n; j++ ) {
#endif
qi += istep, qo += ostep;
t0 = qi[0], tn = qi[idist];
if ( tn == (ModNum)0 ) qo[0] = qo[odist] = t0;
else {
#ifdef POWA_STRIDE
a = powa[aj << k];
#else
a = powa[j << k];
#endif
AxBmodP( tn, ModNum, tn, a, P, Pinv );
qo[0] = AplusBmodP( t0, tn, P );
qo[odist] = A_BmodP( t0, tn, P );
}
}
}
}
#endif
}
/*********************************************************************
*
* NOTE:
* (1) Let \alpha be a primitive N-th root of unity modulo P,
* where N is even,
* 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.
*
* (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,
* because 1 \equiv - 2^n*Q \bmod P.
*
**********************************************************************/
void C_DFT_BACK( in, N, i1, log2N, powa,
#ifdef POWA_STRIDE
a1,
#endif
out, o1, osi, nout, Ninv, P, Pinv, wk )
ModNum in[], powa[], out[], wk[];
int N, log2N, osi, nout, i1, o1;
#ifdef POWA_STRIDE
int a1;
#endif
ModNum P, Ninv;
double Pinv;
/*
* Let K denote log2N. Let N = 2^K, and \alpha be a primitive N-th root of 1.
* This routine computes the inverse discrete-Fourier transform of in[0:N-1][i1]
* using \alpha^{-1} modulo P, and store the `nout' transformed data from the
* `osi'-th in out[0:nout-1][o1]. 0 <= osi < N-1 and 0 < nout <= N.
* The powers of \alpha (to the degrees upto N/2=2^{K-1} at least) are given in powa[*] (or powa[a1* *]).
* The area wk[] is used to contain the intermediate transformed data, and thus,
* must have 2*N entries at least. Notice that `out' cannot be used because
* its space amount (`nout') may not be sufficient for this purpose.
*/
{
int i, j, k, n, is, os, istep, ostep, idist, inwk, halfN;
#ifdef POWA_STRIDE
int aj, halfNa;
#endif
ModNum *qi, *qo, *qis, *qos, a, t0, tn, tt;
/*
* for k = K-1, ..., 0 do the following:
*
* Out[2^{K-k}*i + j] <= t0 + \alpha^{- 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;
* = t0 + \alpha^{N/2 - 2^k * j}*tn;
*
* where t0 = In[2^{K-k-1}*i + j] and
* tn = In[2^{K-k-1}*i + j + 2^{K-1}],
* for 0 <= i < 2^k and 0 <= j < 2^{K-k-1}.
*
*/
if ( log2N == 1 ) {
/* K = 1, k = 0, 0 <= i < 1, 0 <= j < 1.
* Out[0] <= t0 + tn, Out[1] <= t0 - tn,
* where t0 = In[0] and Tn = In[1].
*/
t0 = in[0], tn = in[i1];
if ( osi == 0 ) {
/* out[0] = (t0 + tn)*Ninv; */
tt = tn == (ModNum)0 ? t0 : AplusBmodP( t0, tn, P );
if ( tt != (ModNum)0 ) {
AxBmodP( tt, ModNum, tt, Ninv, P, Pinv );
}
out[0] = tt;
i = 1;
} else i = 0;
/**/
if ( osi + nout >= 2 ) {
/* out[osi == 0 ? 1 : 0] = (t0 - tn)*Ninv; */
tt = tn == (ModNum)0 ? t0 : A_BmodP( t0, tn, P );
if ( tt != (ModNum)0 ) {
AxBmodP( tt, ModNum, tt, Ninv, P, Pinv );
}
out[i] = tt;
}
return;
}
/****/
halfN = n = 1 << (k = log2N - 1); /* == 2^{K-1} */
#ifdef POWA_STRIDE
halfNa = a1 << k;
#endif
/*
* when k = K-1,
* Out[2*i] <= t0 + tn, and Out[2*i+1] <= t0 - tn,
* where t0 = In[i] and tn = In[i+2^{K-1}],
* for 0 <= i < 2^{K-1}.
*/
qi = in, istep = i1, idist = i1 << k, qo = wk, inwk = 0;
for ( i = n; i-- > 0; qi += istep, qo += 2 ) {
t0 = qi[0], tn = qi[idist];
if ( tn == (ModNum)0 ) qo[0] = qo[1] = t0;
else {
qo[0] = AplusBmodP( t0, tn, P );
qo[1] = A_BmodP( t0, tn, P );
}
}
#ifdef EBUG
fprintf( stderr, "::: DFT^{-1} ::: after the first step\n" );
for ( qo = wk, i = 0; i < N; i++ ) {
if ( (i%5) == 0 ) fprintf( stderr, "\t" );
fprintf( stderr, "%10d,", (int)wk[i] );
if ( (i%5) == 4 ) fprintf( stderr, "\n" );
}
if ( (i%5) != 0 ) fprintf( stderr, "\n" );
#endif
/**/
idist = halfN;
for ( ostep = 2; --k > 0; inwk = 1 - inwk ) {
n >>= 1;
istep = ostep; /* == 2^{K-k-1} */
ostep <<= 1; /* == 2^{K-k} */
if ( inwk ) qi = &wk[N], qo = wk;
else qi = wk, qo = &wk[N];
qis = qi, qos = qo;
/*
* for j = 0,
* 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}].
*/
for ( i = n, is = os = 0; i-- > 0; qi += istep, qo += ostep ) {
t0 = qi[0], tn = qi[idist];
if ( tn == (ModNum)0 ) qo[0] = qo[istep] = t0;
else {
qo[0] = AplusBmodP( t0, tn, P );
qo[istep] = A_BmodP( t0, tn, P );
}
}
#ifdef POWA_STRIDE
for ( aj = a1, j = 1; j < istep; aj += a1, j++ ) {
/*** a = P - powa[halfNa - (aj << k)]; ***/
a = powa[halfNa - (aj << k)];
#else
for ( j = 1; j < istep; j++ ) {
/*** a = P - powa[halfN - (j << k)]; ***/
a = powa[halfN - (j << k)];
#endif
qi = ++qis, qo = ++qos;
for ( i = n; i-- > 0; qi += istep, qo += ostep ) {
t0 = qi[0], tn = qi[idist];
if ( tn == (ModNum)0 ) qo[0] = qo[istep] = t0;
else {
AxBmodP( tn, ModNum, tn, a, P, Pinv );
/*** qo[0] = AplusBmodP( t0, tn, P ); ***/
/*** qo[istep] = A_BmodP( t0, tn, P ); ***/
qo[0] = A_BmodP( t0, tn, P );
qo[istep] = AplusBmodP( t0, tn, P );
}
}
}
}
#if 0
final:
#endif
/*
* The final step of k = 0. i can take only the value 0.
*
* Out[j] <= t0 + \alpha^{-j} * tn;
* Out[j + 2^{K-1}] <= t0 - \alpha^{-j} * tn;
*
* where t0 = In[j] and tn = In[j + 2^{K-1}],
* for 0 <= j < 2^{K-1}.
*
*/
qo = out, qi = &wk[inwk ? N : 0];
if ( (n = osi + nout) > N ) nout = (n = N) - osi; /* true nout */
if ( (k = nout - halfN) <= 0 ) { /* no overlap, i.e., no +- */
if ( osi <= 0 ) {
t0 = qi[0], tn = qi[idist];
if ( tn != (ModNum)0 ) t0 = AplusBmodP( t0, tn, P );
if ( t0 == (ModNum)0 ) qo[0] = t0;
else {
AxBmodP( qo[0], ModNum, t0, Ninv, P, Pinv );
}
#ifdef POWA_STRIDE
aj = a1;
#endif
j = 1, k = n, i = 0, qo += o1;
} else {
j = osi;
#ifdef POWA_STRIDE
aj = osi*a1;
#endif
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_upper; /* upper only */
else i = n - halfN + 1, k = halfN; /* we have lower, and i upper */
}
} else {
os = o1*halfN;
#ifdef EBUG
fprintf( stderr, "::::: DFT^{-1} ::: os=%d*%d=%d\n", o1, halfN, os );
#endif
if ( osi <= 0 ) {
t0 = qi[0], tn = qi[idist];
if ( tn == (ModNum)0 ) {
if ( t0 == (ModNum)0 ) tt = t0;
else { AxBmodP( tt, ModNum, t0, Ninv, P, Pinv ); }
qo[0] = qo[os] = tt;
} else {
tt = AplusBmodP( t0, tn, P );
if ( tt == (ModNum)0 ) qo[0] = tt;
else { AxBmodP( qo[0], ModNum, tt, Ninv, P, Pinv ); }
tt = A_BmodP( t0, tn, P );
if ( tt == (ModNum)0 ) qo[os] = tt;
else { AxBmodP( qo[os], ModNum, tt, Ninv, P, Pinv ); }
}
#ifdef EBUG
fprintf( stderr, "::::: DFT^{-1} ::: [0]=%d, [%d]=%d\n", (int)qo[0], os, (int)qo[os] );
#endif
#ifdef POWA_STRIDE
j = 1, aj = a1, n = halfN, i = 0, qo += o1, k--;
} else j = osi, aj = osi*a1, i = (n -= k) - halfN;
/**/
for ( ; k-- > 0; aj += a1, j++, qo += o1 ) {
#else
j = 1, n = halfN, i = 0, qo += o1, k--;
} else j = osi, i = (n -= k) - halfN;
/**/
for ( ; k-- > 0; j++, qo += o1 ) {
#endif
t0 = qi[j], tn = qi[j+idist];
if ( tn == (ModNum)0 ) {
if ( t0 != (ModNum)0 ) {
AxBmodP( t0, ModNum, t0, Ninv, P, Pinv );
}
qo[0] = qo[os] = t0;
} else {
#ifdef POWA_STRIDE
/*** a = P - powa[halfNa - aj]; ***/
a = powa[halfNa - aj];
#else
/*** a = P - powa[halfN - j]; ***/
a = powa[halfN - j];
#endif
AxBmodP( tn, ModNum, tn, a, P, Pinv );
/*** tt = AplusBmodP( t0, tn, P ); ***/
tt = A_BmodP( t0, tn, P );
if ( tt == (ModNum)0 ) qo[0] = tt;
else { AxBmodP( qo[0], ModNum, tt, Ninv, P, Pinv ); }
/*** tt = A_BmodP( t0, tn, P ); ***/
tt = AplusBmodP( t0, tn, P );
if ( tt == (ModNum)0 ) qo[os] = tt;
else { AxBmodP( qo[os], ModNum, tt, Ninv, P, Pinv ); }
}
}
k = halfN;
}
/*
* At this point, k is an upper limit < N/2 for j, i is the number of j's
* > N/2, and n is the real upper limit of j.
*/
#ifdef POWA_STRIDE
for ( ; j < k; aj += a1, j++, qo += o1 ) {
#else
for ( ; j < k; j++, qo += o1 ) {
#endif
t0 = qi[j], tn = qi[j+idist];
if ( tn == (ModNum)0 ) {
if ( t0 == (ModNum)0 ) qo[0] = t0;
else { AxBmodP( qo[0], ModNum, t0, Ninv, P, Pinv ); }
} else {
#ifdef POWA_STRIDE
a = P - powa[halfNa - aj];
#else
a = P - powa[halfN - j];
#endif
AxBplusCmodP( tt, ModNum, tn, a, t0, P, Pinv );
if ( tt == (ModNum)0 ) qo[0] = tt;
else { AxBmodP( qo[0], ModNum, tt, Ninv, P, Pinv ); }
}
}
if ( i <= 0 ) return;
/**/
L_halfN:
t0 = qi[0], tn = qi[idist];
tt = tn == (ModNum)0 ? t0 : A_BmodP( t0, tn, P );
if ( tt == (ModNum)0 ) qo[0] = tt;
else { AxBmodP( qo[0], ModNum, tt, Ninv, P, Pinv ); }
qo += o1, j++;
/**/
L_upper:
#ifdef POWA_STRIDE
for ( n -= halfN, aj = a1, j = 1; j < n; aj += a1, j++, qo += o1 ) {
#else
for ( n -= halfN, j = 1; j < n; j++, qo += o1 ) {
#endif
t0 = qi[j], tn = qi[j+idist];
if ( tn == (ModNum)0 ) {
if ( t0 == (ModNum)0 ) qo[0] = t0;
else { AxBmodP( qo[0], ModNum, t0, Ninv, P, Pinv ); }
} else {
#ifdef POWA_STRIDE
a = powa[halfNa - aj];
#else
a = powa[halfN - j];
#endif
AxBplusCmodP( tt, ModNum, tn, a, t0, P, Pinv );
if ( tt == (ModNum)0 ) qo[0] = tt;
else { AxBmodP( qo[0], ModNum, tt, Ninv, P, Pinv ); }
}
}
}
#if defined(USE_FLOAT)
void C_PREP_ALPHA( r, log2ord, log2k, n, tbl, P, Pinv )
ModNum r, tbl[], P;
int log2ord, log2k, n;
double Pinv;
/*
* Let K and k denote log2ord and log2k, respectively.
* Let r be a primitive (2^K)-th root of unity in Z/(P), where P is a prime.
* Compute a = r^{2^{K-k}}, a primitive (2^k)-th root of unity, and
* its powers to the degrees upto n in tbl[*].
*/
{
int i;
ModNum *q;
double t, w, dp = (double) P;
for ( w = (double) r, i = log2ord - log2k; i > 0; i-- ) {
w *= w;
if ( w < dp ) continue;
w -= (dp*(double)((int) (w*Pinv)));
}
tbl[0] = 1;
tbl[1] = (ModNum)(t = w);
for ( q = &tbl[2], i = n - 1; i > 0; i-- ) {
t *= w;
if ( t >= dp ) t -= (dp*(double)((int) (t*Pinv)));
*q++ = t;
}
}
#else
void C_PREP_ALPHA( r, log2ord, log2k, n, tbl, P, Pinv )
ModNum r, tbl[], P;
int log2ord, log2k, n;
double Pinv;
/*
* Let K and k denote log2ord and log2k, respectively.
* Let r be a primitive (2^K)-th root of unity in Z/(P), where P is a prime.
* Compute a = r^{2^{K-k}}, a primitive (2^k)-th root of unity, and
* its powers to the degrees upto n in tbl[*].
*/
{
int i;
ModNum *q;
ModNum t, w, s;
for ( w = r, i = log2ord - log2k; i > 0; i-- ) {
AxBmodP( t, ModNum, w, w, P, Pinv );
w = t;
}
tbl[0] = 1;
tbl[1] = (ModNum)(t = w);
for ( q = &tbl[2], i = n - 1; i > 0; i-- ) {
AxBmodP( s, ModNum, t, w, P, Pinv );
t = s;
*q++ = t;
}
}
#endif
![[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