Арифметика и числовые алгоритмы

      Здесь приведён код быстрого преобразование Фурье (функция fft2) в поле вычетов по некоторому модулю MOD. Этот код может быть использован для

• Быстрого умножения длинных чисел
  
• Умножения многочленов с целыми коэффициентами.

Artem Voroztsov - 11 Mar 2005

#define LT unsigned long

const LT N   = 0x00000100;  // длина вектора
const LT NX  = 0x000000ff;  // N - 1
const LT MOD = 0x00000101;  // модуль

LT Z = 5;     // Z^n должно пробегать все не нулевые остатки по модулю MOD
LT *ZI;       // ZI[n] = (z^n) mod MOD

// Долгий вариант прямого преобразования Фурье
void ft(LT* y, LT* c)
{
    long i, j;
    for(i = 0; i < N ; i++)
        for(c[i] = 0, j = 0; j < N ; j++ )
            c[i] = (c[i] + ZI[( i * j) & NX ] * y[j] ) % MOD;
}

// Долгий вариант обратного преобразования Фурье
void ift(LT* c, LT* y)
{
    long i,j;
    for(i = 0 ; i < N ; i++)
        for(y[i] = 0, j = 0; j < N ; j++)
            y[i] = (y[i] + MOD * MOD - ZI[ (-i*j) & NX ] * c[j]) % MOD;
}

// Быстрое преобразование фурье
//   I ==  1 — прямое
//   I == -1 — обратное
void fft2(LT *x, int I)
{
register LT  t1,t2,u;
register long i,j,p,l,zl,L,BL;
    L = N; BL = I;
    while(L >= 2){
        l = L/2; u = ZI[0]; zl = 0;
        for(j = 0; j < l; j++, u = ZI[ (zl += BL) & NX])
            for(i = j; i < N ; i += L){
                p  = i+l;
                t1 = (x[i] + x[p]) % MOD;
                t2 = (MOD + x[i] - x[p]) % MOD;
                x[p] = (t2 * u) % MOD;
                x[i] = t1 % MOD;
            }
        L /= 2; BL *= 2;
    }

    // Далее делаем перестановку элементов массива x
    //  x[i] --> x[ REVERSEBITS(i) ]
    j = 0;
    for(i=0; i < NX ; i++){
        if(i > j){ t1=x[j]; x[j]=x[i]; x[i]=t1;}
        l = N/2;
        while(j >=l){ j-=l; l /= 2;}
        j+=l;
      }
}

int main()
{
    int i;
    ZI = (LT*) malloc(N * sizeof(LT));
    for(ZI[0]=1,i=1 ; i < N ; i++)
        ZI[i] = ( Z * ZI[i-1] ) % MOD;
    free(ZI);
}