一切皆整数的FFT

    可以知道，对于计算机的储存和处理，我们使用一个double数对来保存复数，同时运算的时候大量的运用了三角函数，浮点乘法等。撇开低效率不谈，光是其精确度都不会太高！

       对于数字的估计可以知道最终的结果可能非常大(浮点数),之后各种误差开始积累~而且速度越来越慢。

       于是对于n很大的情况，一般来说系数不会太大。

      下面介绍一种vb神教我的思想。

    这样就把 复数对应到了一个整数，之后一切都在mod P系统内考虑。显然如果最终得到的结果是小数的话，那么是无法等价的。

       其中，由于N经常是2的次幂，于是我们可以通过 构造 P = c * (2^k) + 1类的素数来达到目的。

       例如在整系数的多项式乘法中，大整数乘法中。。。这些的最终结果均为整数，于是可以用此来等价。

       下面来举一个例子:

       例子:

      我们按照数论替换的思路来YY看吧……

       首先我们找到一个素数P以及它的原根G:

       P = 1004535809 

       G = 3 

      此时G^(P-1)/8 =  395918948 (mod P)

      之后按照例子的做法带入，而已得到A,B序列:

       A = {21 648936765 369939602853529083 7 128750538 634596211 377855264}

       B = {9 107604863 445555772536309069 3 745698602 558980033 619459092}

       对比发现，对于原来是整数项的，得到的结果相同(想想为什么？)，同时纯复数项的各种不同(为什么，通过这个整数是否有可能还原原复数呢？)

       C = {189 956821336 1512323231000195688 21 603796229 853303436 452794142}

       接下来我们知道是需要将C的每个元素取共轭，然后再做一次FFT,并将结果/N得到最终答案。

       问题来了，复数都被我搞成一坨一坨的整数了，如何“共轭”?

    于是牛逼的事情发生了，只需要第0项不变，其余的通过 k -> N – k(1 <= k <N)的映射交换位置就完成了所谓的“共轭”。

之后就是裸的FFT了，并将得到的结果乘上N的逆元即可。

         http://61.187.179.132:8080/JudgeOnline/showproblem?problem_id=2179

         此题是高精度乘法的题目，使用以上算法能在2000~3000MS通过

         http://61.187.179.132:8080/JudgeOnline/showproblem?problem_id=2194

         此题是类多项式乘法的题目，使用以上算法能在6000ms+通过

         https://www.spoj.pl/problems/TSUM/

         此题为裸FFT + 容斥计数

 ** 注意要点: 当系数的结果过大，需要选择合适的P

【FFT高精度乘法(数论变换)代码】

#include<iostream>
#include<stdio.h>
#include<string.h>
#include<cmath>
#include<string>
#include<map>
using namespace std;
typedef long long LL;
#define Inv(n) powmod( n, P - 2, P)

/*
 * P = C * 2^k + 1 ， P是素数 
 * G 为 原根 
 * 对于 N = 2^w 的 FFT, 在 Zp 中 用 g = G^((P - 1) / N) (mod P) 来代替复 根 e^[ -j(2PI / N)]
 */
int P;
int _g[ 25 ]; 
int BIT_CNT;
int powmod( LL a, int b, int c){
    LL ret = 1 ;
    while(b){
         if(b & 0x1) ret = ret * a % c;
         a = a * a % c;
         b >>= 1;         
    }   
    return ret;
}

bool is_prime( int n){    // 小 心溢出！ 
     int i;
     for(i = 2; i * i <= n; ++ i) if(n % i == 0) return 0;
     return 1;     
}

int getP(int Lim){
     // P = C * 2^21 + 1, P >= Lim
     int c;
     for(c = 3; ; ++ c){
         int t = c << 21 | 1;
         if( is_prime( t ) && t >= Lim) return t;        
     }
     return - 1;     
}
bool is_g( int a, int P){
     int i, p0 = P - 1;
     for(i = 1; i * i <= p0; ++ i)
           if( p0 % i == 0) {
               if( powmod( a, i, P) == 1 && i < p0) return 0;
               if( powmod( a, p0 / i, P) == 1 && p0 / i < p0) return 0;     
           }  
     return 1;   
} 

int getG( int P){
    int g;
    for(g = 2; !is_g( g, P ); ++ g);
    return g;    
}

void get_g(int G, int P, int blim,int _g[]){
    int i, j;
    for(i = 0; i < blim; ++ i){
          j = 1 << i;
          _g[ i ] = powmod( G, (P - 1) / j, P );      
    }    
}
int reverse(int j){
  int i,k;
  k = 0;
  for(i = 0;i < BIT_CNT;++ i) if( j& (1 << i)) k |= 1 << (BIT_CNT - i - 1);
  return k; 
}

void FFT(int x[], int n){
  int i,j,m,t0, t1, i0, j0, tt;
  for(m = 1;m <= BIT_CNT;++ m){
        i0 = 1 << m;
        j0 = i0 >> 1;
    for(i = 0;i < n;i += i0 ) 
     for(j = 0, tt = 1;j < j0 ;++ j, tt = (LL)tt * _g[ m ] % P)  {
       t0 =  tt;
       t1 =  (LL) x[i + j + j0 ] * t0 % P;
      t0 = ( x[i + j] + t1) % P;
      t1 = (x[i + j] - t1) % P; 
            if( t1 < 0) t1 += P;
      x[i + j] = t0;
      x[ i + j + j0 ] = t1;
    } 
  }
  return ;
}

void conv( int a[], int b[], int n) {
     int i;
     FFT( a , n);
     FFT( b , n);
     for(i = 0; i < n; ++ i) b[ i ] = (LL) a[ i ] * b[ i ] % P;
     for(i = 0; i < n; ++ i) a[ reverse( i ) ] = b[ i == 0 ? 0 : n - i ];
     FFT( a , n);
     for(i = 0; i < n; ++ i) a[ i ] = (LL) a[i] * Inv( n ) % P ;      
}

const int maxn = 1 << 19;
char A[ maxn ], B[ maxn ];
int a[ maxn ], b[ maxn ], n;

void init(){
     P = getP( 1000000000 );
     get_g( getG( P) , P, 21, _g) ;    
}

void get(){
     int i,j;
     scanf("%d", &n);    
     scanf("%s%s", A, B);
     int v, c = 0,k = 0;
     int av, bv,t = 1;
     av = bv = 0;
     int on = n / 1 + (n % 1 != 0);
     for( BIT_CNT = 1; on + on > (1 << BIT_CNT); ++ BIT_CNT);
     for(i = n - 1; i >= 0; -- i) {
           av = av + t * (A[ i ] - '0');
           bv = bv + t * (B[ i ] - '0');
           ++ c;
           if(c == 1 || i == 0){
               j = reverse( k );
               a[ j ] = av;
               b[ j ] = bv;
               ++ k;
               c = av = bv = 0; t = 1;       
           }else t *= 10;       
     }
     n = 1 << BIT_CNT;
}

int ans[ maxn ];

void work(){
     int i, j = 0, k;
     conv( a, b, n );
     for(i = 0; i < n; ++ i){
           k = a[ i ] + j;
           ans[ i ] = k % 10;
           j = k / 10;      
     }
     for(i = n - 1; i >= 0 && ans[ i ] == 0; -- i);
     for( printf("%d", ans[i --]); i >= 0; -- i) printf("%d", ans[i]);
     puts("");  
}

int main(){
    init();
    get();
    work();
    //while( 1 ) ;
   return 0;
}