FFT/NTT模板 既 HDU1402 A * B Problem Plus

@(学习笔记)[FFT, NTT]

Problem Description

Calculate A * B.

Input

Each line will contain two integers A and B. Process to end of file.
Note: the length of each integer will not exceed 50000.

Output

For each case, output A * B in one line.

Sample Input

1
2
1000
2

Sample Output

2
2000

Solution

FFT和NTT都是可以的.
FFT: 没有特别要求.

const int N = 1 << 17;
char str1[N], str2[N];

#include<cstdio>
#include<cstring>

int len1, len2;

struct complex
{
    double real, imag;
    
    inline complex(){}

    inline complex(double _real, double _imag)
    {
        real = _real, imag = _imag;
    }
    
    inline friend complex operator *(complex a, complex b)
    {
        return complex(a.real * b.real - a.imag * b.imag, a.imag * b.real + b.imag * a.real);
    }
    
    inline friend complex operator +(complex a, complex b)
    {
        return complex(a.real + b.real, a.imag + b.imag);
    }
    
    inline friend complex operator -(complex a, complex b)
    {
        return complex(a.real - b.real, a.imag - b.imag);
    }
}a[N << 1], b[N << 1];

#include<algorithm>

int len;
int rev[N << 1];

inline void prepare(complex *a, complex *b)
{
    int mx = std::max(len1, len2) << 1;
    len = 1;
    int bit = 0;
    
    while(len < mx)
        len <<= 1, ++ bit;
        
    rev[0] = 0;
    
    for(int i = 1; i < len; ++ i)
        rev[i] = rev[i >> 1] >> 1 | (i & 1) << (bit - 1);
}

inline void reverse(complex *a)
{
    for(int i = 0; i < len; ++ i)
        if(rev[i] < i)
            std::swap(a[i], a[rev[i]]);
}

#include<cmath>

const double PI = acos(-1.0);

inline void fft(complex *a, int opt)
{
    reverse(a);
    
    for(int i = 2; i <= len; i <<= 1)
    {
        complex omega_n = complex(cos(2.0 * PI * (double)opt / (double)i), sin(2.0 * PI * (double)opt / (double)i));
        
        for(int j = 0; j < len; j += i)
        {
            complex omega = complex(1.0, 0.0); //e ^ 0
            
            for(int k = j; k < j + i / 2; ++ k)
            {
                complex x = a[k], y = omega * a[k + i / 2];
                a[k] = x + y, a[k + i / 2] = x - y;
                omega = omega * omega_n;
            }
        }
    }
    
    if(opt == -1)
        for(int i = 0; i < len; ++ i)
            a[i].real /= (double)len;
}

int ans[N << 1];

inline void convolute(complex *a, complex *b)
{
    prepare(a, b);
    fft(a, 1), fft(b, 1);
    
    for(int i = 0; i < len; ++ i)
        a[i] = a[i] * b[i];
        
    fft(a, -1);
    
    for(int i = 0; i < len; ++ i)
        ans[i] = (int)(a[i].real + 0.5);
}

int main()
{   
    #ifndef ONLINE_JUDGE
    freopen("HDU1402.in", "r", stdin);
    freopen("HDU1402.out", "w", stdout);
    #endif

    while(gets(str1) && gets(str2))
    {
        len1 = strlen(str1), len2 = strlen(str2);
        
        memset(a, 0, sizeof(a));
        memset(b, 0, sizeof(b));
        
        for(int i = 0; i < len1; ++ i)
            a[i] = complex(str1[len1 - i - 1] - '0', 0);
            
        for(int i = 0; i < len2; ++ i)
            b[i] = complex(str2[len2 - i - 1] - '0', 0);
            
        convolute(a, b);
        
        for(int i = 0; i < len; ++ i)
            ans[i + 1] += ans[i] / 10, ans[i] %= 10;
        
        for(; len && ! ans[len - 1]; -- len);
        
        for(int i = len - 1; ~ i; -- i)
            putchar(ans[i] + '0');
            
        if(! len)
            putchar('0');
            
        putchar('\n');
    }
}

NTT: 要求取模的数\(P\)是费马素数. 对于原根\(g\), 要确保有 \(g ^ 0 \ne g^1 \ne ... \ne g^{p - 2} (modP)\). NTT的过程相当于用 \(\omega_n = g^{ \frac{p - 1}{n}}\) 替代了 \(\omega_n = e^{\frac{2 \cdot \pi}{n} \cdot i}\)的FFT. \(\omega_i\)和\(inv(\omega_i)\)可以在预处理中求出.

const long long P = (479 << 21) + 1;
const long long G = 3;

inline long long quickPower(long long a, long long k, long long mod)
{
    if(! k)
        return 1;
        
    long long ret = quickPower(a, k >> 1, mod);
    ret = ret * ret % mod;
    
    if(k & 1)
        ret = ret * a % mod;
        
    return ret;
}

const long long QUAN = 1 << 5;
long long omega_[QUAN];

inline void getOmega_n()
{
    for(long long i = 0; i < QUAN; ++ i)
        omega_[i] = quickPower(G, (P - 1) / (1 << i), P);
}

const long long LEN = 1 << 16;
char str1[LEN], str2[LEN];

#include<cstdio>
#include<cstring>

long long len1, len2;
long long a[LEN << 1], b[LEN << 1];

#include<algorithm>

long long len;
long long rev[LEN << 1];

inline void prepare()
{
    long long mx = std::max(len1, len2) << 1;
    len = 1;
    long long bit = 0;
    
    while(len < mx)
        len <<= 1, ++ bit;
        
    rev[0] = 0;
        
    for(long long i = 0; i < len; ++ i)
        rev[i] = rev[i >> 1] >> 1 | (i & 1) << (bit - 1);
}

inline void reverse(long long *a)
{
    for(long long i = 0; i < len; ++ i)
        if(rev[i] < i)
            std::swap(a[i], a[rev[i]]);
}

inline void NTT(long long *a, long long opt)
{
    reverse(a);
    long long temp = 0;
    
    for(long long i = 2; i <= len; i <<= 1)
    {
        ++ temp;
        int omega_i = ~ opt ? omega_[temp] : quickPower(omega_[temp], P - 2, P);
        
        for(long long j = 0; j < len; j += i)
        {
            long long omega = 1;
            
            for(long long k = j; k < j + i / 2; ++ k)
            {
                long long u = a[k], t = omega * a[k + i / 2] % P;
                a[k] = (u + t) % P, a[k + i / 2] = (u - t + P) % P;
                
                omega = omega * omega_i % P;
            }
        }
    }
    
    if(opt == -1)
    {
        long long inv = quickPower(len, P - 2, P);
        
        for(long long i = 0; i < len; ++ i)
            a[i] = a[i] * inv % P;
    }
}

long long ans[LEN << 1];

inline void convolute(long long *a, long long *b)
{
    prepare();
    NTT(a, 1), NTT(b, 1);
    
    for(long long i = 0; i < len; ++ i)
        a[i] = a[i] * b[i] % P;
        
    NTT(a, -1);
    
    for(int i = 0; i < len; ++ i)
        ans[i] = a[i];
}

int main()
{
    #ifndef ONLINE_JUDGE
    freopen("HDU1402.in", "r", stdin);
    #endif
    
    getOmega_n();
    
    while(gets(str1) && gets(str2))
    {
        len1 = strlen(str1), len2 = strlen(str2);
        
        memset(a, 0, sizeof(a)), memset(b, 0, sizeof(b));
        
        for(long long i = 0; i < len1; ++ i)
            a[len1 - i - 1] = str1[i] - '0';
            
        for(long long i = 0; i < len2; ++ i)
            b[len2 - i - 1] = str2[i] - '0';
            
        convolute(a, b);
        
        for(int i = 0; i < len; ++ i)
            ans[i + 1] += ans[i] / 10, ans[i] %= 10;
        
        for(; len && ! ans[len - 1]; -- len);
        
        for(int i = len - 1; ~ i; -- i)
            putchar(ans[i] + '0');
            
        if(! len)
            putchar('0');
            
        putchar('\n');
    }
}

转载于:https://www.cnblogs.com/ZeonfaiHo/p/6552262.html