天天看點
傳回

FFT 快速傅裡葉變換 && NTT 快速數論變換

SDNU 1531 a*b III (FFT模闆)

Description

計算a乘b，多組輸入（50組以内）。

Input

輸入a b，資料範圍0 <= a,b <= 10^100000。

Output

輸出a與b的乘積。

Sample Input

2 2

4 4

Sample Output

4

16

Hint

FFT

FFT的原理是。。。（亂記的），大概原理懂了

FFT 快速傅裡葉變換 &amp;&amp; NTT 快速數論變換
FFT 快速傅裡葉變換 &amp;&amp; NTT 快速數論變換

AC代碼：

#include<bits/stdc++.h>

using namespace std;

const double PI = acos(-1.0);
struct Complex
{
    double x,y;
    Complex(double _x = 0.0,double _y = 0.0)
    {
        x = _x;
        y = _y;
    }
    Complex operator -(const Complex &b)const
    {
        return Complex(x-b.x,y-b.y);
    }
    Complex operator +(const Complex &b)const
    {
        return Complex(x+b.x,y+b.y);
    }
    Complex operator *(const Complex &b)const
    {
        return Complex(x*b.x-y*b.y,x*b.y+y*b.x);
    }
};

void change(Complex y[],int len)
{
    int i,j,k;
    for(i = 1, j = len/2; i <len-1; i++)
    {
        if(i < j)
            swap(y[i],y[j]);
        k = len/2;
        while(j >= k)
        {
            j -= k;
            k /= 2;
        }
        if(j < k)
            j += k;
    }
}

void fft(Complex y[],int len,int on)
{
    change(y,len);
    for(int h = 2; h <= len; h <<= 1)
    {
        Complex wn(cos(-on*2*PI/h), sin(-on*2*PI/h));
        for(int j = 0; j < len; j+=h)
        {
            Complex w(1,0);
            for(int k = j; k < j+h/2; k++)
            {
                Complex u = y[k];
                Complex t = w*y[k+h/2];
                y[k] = u+t;
                y[k+h/2] = u-t;
                w = w*wn;
            }
        }
    }
    if(on == -1)
        for(int i = 0; i < len; i++)
            y[i].x /= len;
}

const int MAXN = 200010;
Complex x1[MAXN],x2[MAXN];
char str1[MAXN/2],str2[MAXN/2];
int sum[MAXN];

int main()
{
    while(~scanf("%s%s",str1,str2))
    {
        int len1 = strlen(str1);
        int len2 = strlen(str2);
        int len = 1;
        while(len < len1*2 || len < len2*2) len<<=1;

        for(int i = 0; i < len1; i++) x1[i] = Complex(str1[len1-1-i]-'0',0);

        for(int i = len1; i < len; i++) x1[i] = Complex(0,0);

        for(int i = 0; i < len2; i++) x2[i] = Complex(str2[len2-1-i]-'0',0);

        for(int i = len2; i < len; i++) x2[i] = Complex(0,0);

        fft(x1,len,1);
        fft(x2,len,1);

        for(int i = 0; i < len; i++) x1[i] = x1[i]*x2[i];

        fft(x1,len,-1);
        for(int i = 0; i < len; i++) sum[i] = (int)(x1[i].x+0.5);

        for(int i = 0; i < len; i++)
        {
            sum[i+1]+=sum[i]/10;
            sum[i]%=10;
        }
        len = len1+len2-1;

        while(sum[len] <= 0 && len > 0) len--;

        for(int i = len; i >= 0; i--) printf("%c",sum[i]+'0');

        printf("\n");
    }
    return 0;
}

SDNU 1532 a*b IV（NTT模闆）

Description

計算a*b，多組輸入（50組以内）。

Input

兩個數a,b,資料範圍 0<=a,b<=10^100000。

Output

輸出a與b的乘積。

Sample Input

2 2

4 4

Sample Output

4

16

Hint

NTT

AC代碼：

#include <bits/stdc++.h>

using namespace std;
typedef long long ll;

//拓展歐幾裡得
void exgcd(ll a, ll b, ll &x, ll &y)
{
    if (b == 0)
    {
        x = 1;
        y = 0;
        return;
    }
    ll x0, y0;
    exgcd(b, a % b, x0, y0);
    x = y0;
    y = x0 - (ll)(a / b) * y0;
}

//求逆元
ll Inv(ll a, ll p)
{
    ll x, y;
    exgcd(a, p, x, y);
    x %= p;
    while (x < 0) x += p;
    return x;
}

//快速幂取模
ll qpow(ll a, ll b, ll p)
{
    if (b < 0)
    {
        b = -b;
        a = Inv(a, p);
    }
    ll ans = 1, mul = a % p;
    while (b)
    {
        if (b & 1) ans = ans * mul % p;
        mul = mul * mul % p;
        b >>= 1;
    }
    return ans;
}

//在模p意義下的計算
#define maxn (65537*2)
const int MOD = 479 * (1 << 21) + 1, G = 3;
//const ll MOD = 15 * (1 << 27) + 1, G = 31;

//翻轉數組
ll rev[maxn];

void get_rev(ll bit)
{
    for (ll i = 0; i < (1 << bit); i++)
    {
        rev[i] = (rev[i >> 1] >> 1) | ((i & 1) << (bit - 1));
    }
}

//存儲數組
ll ar[maxn], br[maxn];

//快速數論變換
void ntt(ll *a, ll n, ll dft)
{
    //翻轉
    for (ll i = 0; i < n; i++)
    {
        if (i < rev[i]) swap(a[i], a[rev[i]]);
    }
    //蝴蝶操作模拟
    for (ll step = 1; step < n; step <<= 1)
    {
        ll wn;
        wn = qpow(G, dft * (MOD - 1) / (step * 2), MOD);
        for (ll j = 0; j < n; j += (step << 1))
        {
            ll wnk = 1;//這裡一定要用long long不然會迷之溢出
            for (ll k = j; k < j + step; k++)
            {
                ll x = a[k] % MOD, y = (wnk * a[k + step]) % MOD;//這裡也要用long long
                a[k] = (x + y) % MOD;
                a[k + step] = ((x - y) % MOD + MOD) % MOD;
                wnk = (wnk * wn) % MOD;
            }
        }
    }
    if (dft == -1)
    {
        ll nI = Inv(n, MOD);
        for (ll i = 0; i < n; i++) a[i] = a[i] * nI % MOD;
    }
}

//輸入數組
char s1[maxn], s2[maxn];

int main()
{
    while (~scanf("%s%s", s1, s2)) {
        ll l1 = strlen(s1), l2 = strlen(s2);
        for (ll i = 0; i < l1; i++) ar[i] = s1[l1 - i - 1] - '0';
        for (ll i = 0; i < l2; i++) br[i] = s2[l2 - i - 1] - '0';

        ll bit, s = 2;
        for (bit = 1; (1 << bit) < (l1 + l2 - 1); bit++) s <<= 1;

        get_rev(bit);
        ntt(ar, s, 1);
        ntt(br, s, 1);

        for (ll i = 0; i < s; i++) ar[i] = ar[i] * br[i] % MOD;

        ntt(ar, s, -1);
        for (ll i = 0; i < s; i++)
        {
            ar[i + 1] += ar[i] / 10;
            ar[i] %= 10;
        }
        ll cnt = s;

        while (cnt >= 0 && ar[cnt] == 0) cnt--;

        for (ll i = cnt; i >= 0; i--) printf("%lld", ar[i]);

        if (cnt == -1) putchar('0');

        putchar('\n');

        memset(ar, 0, sizeof(ar));
        memset(br, 0, sizeof(br));
    }
    return 0;
}

21牛客暑假多校訓練營第一場H題Hash Function

FFT 快速傅裡葉變換 &amp;&amp; NTT 快速數論變換
FFT 快速傅裡葉變換 &amp;&amp; NTT 快速數論變換
FFT 快速傅裡葉變換 &amp;&amp; NTT 快速數論變換

題意：

輸入n，之後n個數，找最小的一個正整數p滿足這n個數模p取得數字都不相同。

分析：

對于任意

1<=i，j<=n&&i!=j

滿足

ai%p!=aj%p

,即

|ai-aj|%p != 0

,一旦有了任意兩個數的內插補點，之後暴力枚舉p即可。

暴力求任意兩個數的內插補點是

O(n^2)

,比較快的方法是轉化為多項式乘法，構造兩個多項式 ∑ 1 n x a i \sum_{1}^{n}{x^{a_i}} ∑1n​xai​和 ∑ 1 n x − a i \sum_{1}^{n}{x^{-a_i}} ∑1n​x−ai​,讓兩個多項式相乘，得到多項式 ∑ x b i \sum{x^{b_i}} ∑xbi​, ∣ b i ∣ |b_i| ∣bi​∣的所有取值就是任意兩個數字的內插補點（0除外，因為保證了 a i ! = b i a_i!=b_i ai​!=bi​）。多項式相乘的過程可以通過FFT或者NTT實作，由于數組下标無法是0，是以第二個多項式可以構造為 ∑ 1 n x M X − a i \sum_{1}^{n}{x^{MX-a_i}} ∑1n​xMX−ai​即将指數加上一個偏移量MX，偏移量MX應該設定為一個稍大于5e5的資料（不懂原因，我覺得ai最大值為5e5，MX設定為5e5即可，但是這樣确實被一個樣例卡了，duipai了很久也沒找出來）。

之後暴力枚舉每個數字p，判斷p的倍數是否出現即可，這個過程是 O ( n l o n g n ) O(nlongn) O(nlongn)的。

AC代碼：

FFT

#include<bits/stdc++.h>

using namespace std;

const double PI = acos(-1.0);
struct Complex
{
    double x,y;
    Complex(double _x = 0.0,double _y = 0.0)
    {
        x = _x;
        y = _y;
    }
    Complex operator -(const Complex &b)const
    {
        return Complex(x-b.x,y-b.y);
    }
    Complex operator +(const Complex &b)const
    {
        return Complex(x+b.x,y+b.y);
    }
    Complex operator *(const Complex &b)const
    {
        return Complex(x*b.x-y*b.y,x*b.y+y*b.x);
    }
};

void change(Complex y[],int len)
{
    int i,j,k;
    for(i = 1, j = len/2; i <len-1; i++)
    {
        if(i < j)
            swap(y[i],y[j]);
        k = len/2;
        while(j >= k)
        {
            j -= k;
            k /= 2;
        }
        if(j < k)
            j += k;
    }
}

void fft(Complex y[],int len,int on)
{
    change(y,len);
    for(int h = 2; h <= len; h <<= 1)
    {
        Complex wn(cos(-on*2*PI/h), sin(-on*2*PI/h));
        for(int j = 0; j < len; j+=h)
        {
            Complex w(1,0);
            for(int k = j; k < j+h/2; k++)
            {
                Complex u = y[k];
                Complex t = w*y[k+h/2];
                y[k] = u+t;
                y[k+h/2] = u-t;
                w = w*wn;
            }
        }
    }
    if(on == -1)
        for(int i = 0; i < len; i++)
            y[i].x /= len;
}

const int MAXN = 1050010;
const int MX = 500005;
Complex x1[MAXN],x2[MAXN];
int n;
int ar[MAXN>>1];
int sum[MAXN];
int len1, len2;
bool boo1[MAXN>>1], boo2[MAXN>>1];
bool vis[MAXN>>1];
bool flag;

int main()
{
    scanf("%d", &n);
    for(int i = 1; i <= n; ++i)
    {
        scanf("%d", &ar[i]);
        boo1[ar[i]] = true;
        boo2[MX-ar[i]] = true;
    }

    if(n == 1) 
    {
        printf("1\n");
        return 0;
    }
    
    len1 = len2 = MX;
    int len = 1;
    while(len < len1*2 || len < len2*2) len<<=1;

    for(int i = 1; i <= n; ++i) x1[ar[i]] = Complex(1, 0);
    for(int i = 1; i <= n; ++i) x2[MX-ar[i]] = Complex(1, 0);

    for(int i = 0; i < len; i++)
    {
        if(!boo1[i]) x1[i] = Complex(0, 0);
        if(!boo2[i]) x2[i] = Complex(0, 0);
    }

    fft(x1,len,1);
    fft(x2,len,1);

    for(int i = 0; i < len; i++) x1[i] = x1[i]*x2[i];

    fft(x1,len,-1);
    for(int i = 0; i < len; i++) sum[i] = (int)(x1[i].x+0.5);

    for(int i = 1; i <= MX; ++i) if(sum[i+MX] != 0) vis[i] = true;

    for(int i = 1; i <= MX; ++i)
    {
        flag = false;
        for(int j = i; j <= MX; j += i)
        {
            if(vis[j])
            {
                flag = true;
                break;
            }
        }
        if(!flag)
        {
            printf("%d\n", i);
            break;
        }
    }

    return 0;
}

NTT

#include <bits/stdc++.h>

using namespace std;
typedef long long ll;

//拓展歐幾裡得
void exgcd(ll a, ll b, ll &x, ll &y)
{
    if (b == 0)
    {
        x = 1;
        y = 0;
        return;
    }
    ll x0, y0;
    exgcd(b, a % b, x0, y0);
    x = y0;
    y = x0 - (ll)(a / b) * y0;
}

//求逆元
ll Inv(ll a, ll p)
{
    ll x, y;
    exgcd(a, p, x, y);
    x %= p;
    while (x < 0) x += p;
    return x;
}

//快速幂取模
ll qpow(ll a, ll b, ll p)
{
    if (b < 0)
    {
        b = -b;
        a = Inv(a, p);
    }
    ll ans = 1, mul = a % p;
    while (b)
    {
        if (b & 1) ans = ans * mul % p;
        mul = mul * mul % p;
        b >>= 1;
    }
    return ans;
}

//在模p意義下的計算
#define maxn (2000005)
const int MOD = 479 * (1 << 21) + 1, G = 3;
//const ll MOD = 15 * (1 << 27) + 1, G = 31;

//翻轉數組
ll rev[maxn];

void get_rev(ll bit)
{
    for (ll i = 0; i < (1 << bit); i++)
    {
        rev[i] = (rev[i >> 1] >> 1) | ((i & 1) << (bit - 1));
    }
}

//存儲數組
ll ar[maxn], br[maxn];

//快速數論變換
void ntt(ll *a, ll n, ll dft)
{
    //翻轉
    for (ll i = 0; i < n; i++)
    {
        if (i < rev[i]) swap(a[i], a[rev[i]]);
    }
    //蝴蝶操作模拟
    for (ll step = 1; step < n; step <<= 1)
    {
        ll wn;
        wn = qpow(G, dft * (MOD - 1) / (step * 2), MOD);
        for (ll j = 0; j < n; j += (step << 1))
        {
            ll wnk = 1;//這裡一定要用long long不然會迷之溢出
            for (ll k = j; k < j + step; k++)
            {
                ll x = a[k] % MOD, y = (wnk * a[k + step]) % MOD;//這裡也要用long long
                a[k] = (x + y) % MOD;
                a[k + step] = ((x - y) % MOD + MOD) % MOD;
                wnk = (wnk * wn) % MOD;
            }
        }
    }
    if (dft == -1)
    {
        ll nI = Inv(n, MOD);
        for (ll i = 0; i < n; i++) a[i] = a[i] * nI % MOD;
    }
}

//輸入數組
const int MX = 500005;
char s1[maxn], s2[maxn];
int n;
int arr[maxn];
bool boo1[MX + 5], boo2[MX + 5];
bool vis[MX + 5];
bool flag;

int main()
{
    scanf("%d", &n);
    for(int i = 1; i <= n; ++i)
    {
        scanf("%d", &arr[i]);
        boo1[arr[i]] = true;
        boo2[MX-arr[i]] = true;
    }

    ll l1 = MX + 1, l2 = MX + 1;
    for(int i = 1; i <= n; ++i)
    {
        ar[arr[i]] = 1;
        br[MX-arr[i]] = 1;
    }

    ll bit, s = 2;
    for (bit = 1; (1 << bit) < (l1 + l2 - 1); bit++) s <<= 1;

    get_rev(bit);
    ntt(ar, s, 1);
    ntt(br, s, 1);

    for (ll i = 0; i < s; i++) ar[i] = ar[i] * br[i] % MOD;

    ntt(ar, s, -1);

    for(int i = 1; i <= MX; ++i) if(ar[i+MX] != 0) vis[i] = true;

    for(int i = 1; i <= MX; ++i)
    {
        flag = false;
        for(int j = i; j <= MX; j += i)
        {
            if(vis[j])
            {
                flag = true;
                break;
            }
        }
        if(!flag)
        {
            printf("%d\n", i);
            break;
        }
    }

    return 0;
}
SDNUOJ 算法 傅裡葉變換 fft NTT ACM 牛客暑假多校
上一篇: 愛麗絲和鮑勃正在玩遊戲。 有兩堆牌。 每一堆都有N張牌，每張牌都有一張分數。 他們輪流從任意一堆拿起頂部或底部牌，并且卡的分數将被添加到他的總分中。 愛麗絲和鮑勃都非常聰明，并且會拿起卡片以獲得盡可能
下一篇: native bayes

繼續閱讀