FFT,NTT学习笔记

最新推荐文章于 2023-04-29 11:43:51 发布
最新推荐文章于 2023-04-29 11:43:51 发布
阅读量1.4k 收藏
点赞数
分类专栏： FFT && NTT
版权声明：本文为博主原创文章，遵循 CC 4.0 BY-SA 版权协议，转载请附上原文出处链接和本声明。
本文链接：https://blog.csdn.net/morejarphone/article/details/50899112
版权

快速傅里叶变换,可以将多项式相乘的时间复杂度从最简单的O(n^2)优化到O(nlgn),详细过程参考算法导论.

FFT的流程大致是:

1):构造多项式,复杂度O(n)

2):求两个多项式的DFT,复杂度O(nlgn)

3):构造多项式乘积的点值表达式,复杂度O(n)

4):求点值表达式的IDFT,复杂度O(nlgn).


下面是两道最简单的习题:

HDU 1402:点击打开链接

求两个大数乘积.

因为一个大数可以看成是一个多项式,每一位上的值都表示对应次数下的系数,因此可以用FFT加速.

本体的一个坑点就是

len = l1+l2-1;
这句代码,可能是精度问题在len更加高位的地方出现了非0值. 

#include <bits/stdc++.h>
using namespace std;
#define pi acos (-1)
#define maxn 200010

struct plex {
    double x, y;
    plex (double _x = 0.0, double _y = 0.0) : x (_x), y (_y) {}
    plex operator + (const plex &a) const {
        return plex (x+a.x, y+a.y);
    }
    plex operator - (const plex &a) const {
        return plex (x-a.x, y-a.y);
    }
    plex operator * (const plex &a) const {
        return plex (x*a.x-y*a.y, x*a.y+y*a.x);
    }
};

void change (plex *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(plex y[],int len,int on)
{
    change(y,len);
    for(int h = 2; h <= len; h <<= 1)
    {
        plex wn(cos(-on*2*pi/h),sin(-on*2*pi/h));
        for(int j = 0;j < len;j+=h)
        {
            plex w(1,0);
            for(int k = j;k < j+h/2;k++)
            {
                plex u = y[k];
                plex 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;
}
char a[maxn], b[maxn];
plex x1[maxn], x2[maxn];
int ans[maxn];

int main () {
    while (scanf ("%s%s", a, b) == 2) {
        int len = 2, l1 = strlen (a), l2 = strlen (b);
        while (len < l1*2 || len < l2*2)
            len <<= 1;
        for (int i = 0; i < l1; i++) {
            x1[i] = plex (a[l1-1-i]-'0', 0);
        }
        for (int i = l1; i < len; i++)
            x1[i] = plex (0, 0);
        for (int i = 0; i < l2; i++) {
            x2[i] = plex (b[l2-1-i]-'0', 0);
        }
        for (int i = l2; i < len; i++)
            x2[i] = plex (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++) {
            ans[i] = (int)(x1[i].x+0.5);
        }
        for (int i = 0; i < len; i++) {
            if (ans[i] >= 10) {
                ans[i+1] += ans[i]/10;
                ans[i] %= 10;
            }
        }
        len = l1+l2-1;
        while (ans[len] <= 0 && len > 0)
            len--;
        for (int i = len; i >= 0; i--) {
            printf ("%d", ans[i]);
        }
        printf ("\n");
    }
    return 0;
}

HDU 4609: 点击打开链接

题意是给出n个长度,任取3个求能组成三角形的概率.

首先记录下每个长度的数量,然后用FFT加速求出任取两个长度下的情况,这里面有重复:

首先减去两次都取同一根的情况,减完之后的结果都/2.

最后只需要所有的情况减去不能组成三角形的情况,将最初的长度序列排序后从小到大枚举下标,假设这条边是最长边,那么如果所有两条边长度小于等于这条边的情况就应该减去,这里用前缀和统计下就好了.

#include <bits/stdc++.h>
using namespace std;
#define pi acos (-1)
#define maxn 611111

struct plex {
    double x, y;
    plex (double _x = 0.0, double _y = 0.0) : x (_x), y (_y) {}
    plex operator + (const plex &a) const {
        return plex (x+a.x, y+a.y);
    }
    plex operator - (const plex &a) const {
        return plex (x-a.x, y-a.y);
    }
    plex operator * (const plex &a) const {
        return plex (x*a.x-y*a.y, x*a.y+y*a.x);
    }
};

void change (plex y[], int len) {
    if (len == 1)
        return ;
    plex a1[len/2], a2[len/2];
    for (int i = 0; i < len; i += 2) {
        a1[i/2] = y[i];
        a2[i/2] = y[i+1];
    }
    change (a1, len>>1);
    change (a2, len>>1);
    for (int i = 0; i < len/2; i++) {
        y[i] = a1[i];
        y[i+len/2] = a2[i];
    }
    return  ;
}

void fft(plex y[],int len,int on)
{
    change(y,len);
    for(int h = 2; h <= len; h <<= 1)
    {
        plex wn(cos(on*2*pi/h),sin(on*2*pi/h));
        for(int j = 0;j < len;j+=h)
        {
            plex w(1,0);
            for(int k = j;k < j+h/2;k++)
            {
                plex u = y[k];
                plex 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;
}
long long num[maxn], sum[maxn];
int a[maxn];
plex x[maxn];
long long n;

int main () {
    //freopen ("in.txt", "r", stdin);
    int t;
    scanf ("%d", &t);
    while (t--) {
        scanf ("%lld", &n);
        long long Max = 0;
        memset (num, 0, sizeof num);
        for (int i = 1; i <= n; i++) {
            scanf ("%d", &a[i]);
            num[a[i]]++;
            Max = max (Max, (long long)a[i]);
        }
        Max++;
        int len = 2;
        while (len < Max*2)
            len <<= 1;
        for (int i = 0; i < len; i++) {
            x[i] = plex (num[i], 0);
        }
        fft (x, len, 1);
        for (int i = 0; i < len; i++) {
            x[i] = x[i]*x[i];
        }
        fft (x, len, -1);
        for (int i = 0; i < len; i++) {
            num[i] = (long long) (x[i].x+0.5);
        }
        for (int i = 1; i <= n; i++) {//两次取同一个
            num[a[i]+a[i]]--;
        }
        for (int i = 0; i < len; i++) {//重复计算
            num[i] /= 2;
        }
        sum[0] = 0;
        for (int i = 1; i < len; i++) {
            sum[i] = sum[i-1]+num[i];
        }
        sort (a+1, a+1+n);
        long long tot = n*(n-1)*(n-2)/6, ans = tot;
        for (int i = 3; i <= n; i++) {
            ans -= sum[a[i]];
        }
        printf ("%.7f\n", ans*1.0/tot);
    }
    return 0;
}


但是FFT有一个很致命的弱点就是会产生精度误差，在换成long double都不行的时候就需要用到NTT。

NTT就是用数论域中的原根代替FFT中的单位负根，其他的代码完全相同。


求原根的代码：

#include <cstdio>
#include <cmath>
#include <algorithm>
#include <iostream>
#include <vector>
#include <cstring>
using namespace std;
int P;
const int NUM = 32170;
int prime[NUM/4];
bool f[NUM];
int pNum = 0;
void getPrime () {//线性筛选素数
    for (int i = 2; i < NUM; ++ i) {
        if (!f[i]) {
            f[i] = 1;
            prime[pNum++] = i;
        }
        for (int j = 0; j < pNum && i*prime[j] < NUM; ++ j) {
            f[i*prime[j]] = 1;
            if (i%prime[j] == 0) {
                break;
            }
        }
    }
}
long long getProduct(int a,int b,int P) {//快速求次幂mod
    long long ans = 1;
    long long tmp = a;
    while (b)
    {
        if (b&1)
        {
            ans = ans*tmp%P;
        }
        tmp = tmp*tmp%P;
        b>>=1;
    }
    return ans;
}

bool judge (int num) {//求num的所有的质因子
    int elem[1000];
    int elemNum = 0;
    int k = P - 1;
    for (int i = 0; i < pNum; ++ i) {
        bool flag = false;
        while (!(k%prime[i])) {
            flag = true;
            k /= prime[i];
        }
        if (flag) {
            elem[elemNum ++] = prime[i];
        }
        if (k == 1) {
            break;
        }
        if (k/prime[i]<prime[i]) {
            elem[elemNum ++] = prime[i];
            break;
        }
    }
    bool flag = true;
    for (int i = 0; i < elemNum; ++ i) {
        if (getProduct (num, (P-1)/elem[i], P) == 1) {
            flag = false;
            break;
        }
    }
    return flag;
}

int main()  
{  
      
    getPrime();  
    while (cin >> P)  
    {  
        for (int i = 2;;++i)  
        {  
            if (judge(i))  
            {  
                cout << i<< endl;  
                break;  
            }  
        }  
    }  
    return 0;  
}



HDU 1402：

随便选一个不太大的模数和他的原根就好了。

#include <cstdio>
#include <cstring>
#include <algorithm>
#include <iostream>
#include <cmath>
#include <map>
#include <vector>
#include <stack>
using namespace std;
#define mod 1004535809LL
#define G 3
#define maxn 400005

long long qpow (long long a, long long b) {  
    long long ret=1;  
    while (b) {  
        if (b&1) ret = (ret*a)%mod;  
        a = (a*a)%mod;  
        b >>= 1;  
    }  
    return ret;  
}  

void change (long long y[], int len) {
    for(int i = 1, j = len / 2; i < len - 1; i++) {
        if(i < j) swap(y[i], y[j]);
        int k = len / 2;
        while(j >= k) {
            j -= k;
            k /= 2;
        }
        if(j < k) j += k;
    }
}
void ntt(long long y[], int len, int on) {
    change (y, len);
    for(int h = 2; h <= len; h <<= 1) {
        long long wn = qpow(G, (mod-1)/h);
        if(on == -1) wn = qpow(wn, mod-2);
        for(int j = 0; j < len; j += h) {
            long long w = 1;
            for(int k = j; k < j + h / 2; k++) {
                long long u = y[k];
                long long t = w * y[k + h / 2] % mod;
                y[k] = (u + t) % mod;
                y[k+h/2] = (u - t + mod) % mod;
                w = w * wn % mod;
            }
        }
    }
    if(on == -1) {
        long long t = qpow (len, mod - 2);
        for(int i = 0; i < len; i++)
            y[i] = y[i] * t % mod;
    }
}

char a[maxn], b[maxn];
long long x1[maxn], x2[maxn];
long long ans[maxn];

int main () {
    while (scanf ("%s%s", a, b) == 2) {
        int len = 2, l1 = strlen (a), l2 = strlen (b);
        while (len < l1*2 || len < l2*2)
            len <<= 1;
        //cout << len << endl;
        for (int i = 0; i < l1; i++) {
            x1[i] = a[l1-1-i]-'0';
        }
        for (int i = l1; i < len; i++)
            x1[i] = 0;
        for (int i = 0; i < l2; i++) {
            x2[i] = b[l2-1-i]-'0';
        }
        for (int i = l2; i < len; i++)
            x2[i] = 0;
        ntt(x1, len, 1);
        ntt(x2, len, 1);
        for (int i = 0; i < len; i++)
            x1[i] = x1[i]*x2[i]%mod;
        ntt(x1, len, -1);
        for (int i = 0; i < len; i++) {
            ans[i] = x1[i];
        }
        for (int i = 0; i < len; i++) {
            if (ans[i] >= 10) {
                ans[i+1] += ans[i]/10;
                ans[i] %= 10;
            }
        }
        len = l1+l2-1;
        while (ans[len] <= 0 && len > 0)
            len--;
        for (int i = len; i >= 0; i--) {
            printf ("%lld", ans[i]);
        }
        printf ("\n");
    }
    return 0;
}


morejarphone
关注
专栏目录
NTT学习笔记
hipamp的博客
06-17 533
前言 NTTNTTNTT 其实和 FFTFFTFFT 差不多，只不过 NTTNTTNTT 可以在某些模意义下进行运算。 学习 NTTNTTNTT 之前，你得先了解了解 原根 关于模数 NTTNTTNTT 模数是 p=qn+1p=qn+1p=qn+1（满足 ppp 是个质数），其中 n=2kn=2^kn=2k。 为什么呢？ 这样，当长度为 2l(l<k)2^l(l<k)2l(l<k) 时，存在一个单位根 gp−12lg^{\frac{p-1}{2^l}}g2lp−1​，使其满足与 nnn 次
任意模数ntt_再探快速傅里叶变换(FFT)学习笔记(其三)(循环卷积的Bluestein算法+分治FFT+FFT的优化+任意模数NTT)...
weixin_39864261的博客
01-17 820
再探快速傅里叶变换(FFT)学习笔记(其三)(循环卷积的Bluestein算法+分治FFT+FFT的优化+任意模数NTT)写在前面为了不使篇幅过长，预计将把学习笔记分为四部分:DFT,IDFT,FFT的定义,实现与证明:快速傅里叶变换(FFT)学习笔记(其一)任意模数NTT与FFT的优化技巧多项式相关操作一些约定\([p(x)]=\begin{cases}1,p(x)为真 \\ 0,p(x)为假 ...
FFT&NTT
weixin_34150503的博客
08-30 126
https://winniechen.cn/?p=377 内容有些多，我就不复制了... 转载于:https://www.cnblogs.com/Winniechen/p/9562408.html
FFT和NTT
juruo_hejiarui的博客
09-13 1450
FFT 这个东西比较复杂，可能需要看得久一点，而且本人写的比较草率，极有可能让大家看不懂。 问题 现在我们需要解决一个加法卷积问题： 已知两个数组f,gf,gf,g，每个值为实数，大小都为nnn,现在需要求出一个新的数组hhh，大小为2n2n2n，其中hih_ihi​用如下方式计算： hi=∑j=0igjfi−j h_i=\sum\limits_{j=0}^i g_j f_{i-j} hi​=j=0∑i​gj​fi−j​ 显然这时时间复杂度是O(n2)O(n^2)O(n2)的，非常不优秀。需要我们进行优化。
【FFT&NTT入门】大整数乘法
solego’s blog
04-02 735
【FFT&NTT入门】大整数乘法
学习笔记 - 位逆序置换 / 非递归FFT / 快速数论变换(NTT)
qq_26087481的博客
03-28 2265
文章讨论了快速傅立叶变换的非递归实现以及快速数论变换的实现.
快速傅里叶变换（FFT）及快速数论变换（NTT）入门笔记
qq_39561693的博客
11-28 735
个人博客同步食用~ 本博客学习于洛谷P3803各路大佬题解，所以你肯定会找到很多相似的地方QWQ 什么是FFT 傅里叶变换 傅里叶说明了：一切函数/波形都可以用有限或无限个弦函数/波形叠加形成。 例如，现在有一个由三个不同正弦波组成的近似矩形波，将这些波形关于 xxx 轴平行地放在一起，就可以得到下图： 其中红色的是合成出来的近似矩形波。 （当然，图肯定是不太标准的，意思一下就行） 从 xOzxOzxOz 平面看，我们可以获得各个波的时域信号，也就是 t−f(t)t-f(t)t−f(t) 图像。而从 zO
[学习笔记][省选算法]多项式乘法之 FFT & NTT 及卷积
女装的你，如此好看！
06-29 1446
一、开头 （Place ：山东省神犇协会第 998244353 会议厅） SD 神犇 998244353 号（会长）：所有神犇，洛谷黑题都切完了吗？ 全体神犇：切完了！ 神犇 1004535809 号： 998244353 ，我昨天看到一个弱鸡，洛谷名叫 xyz32768 ，他一道黑题都没做，您认为应该怎样呢？ 神犇 998244353 号：这个问题您们考虑一下，应该用怎样的方式把 xy...
FFT与NTT
05-23 158
讲解:http://www.cnblogs.com/poorpool/p/8760748.html 递归版FFT #include <iostream> #include <cstdio> #include <cstring> #include <algorithm> #include <cstdlib> #include &l...
ntt学习笔记
DJ的博客
03-07 241
差不多就是模意义下的fft各个注意事项 各个注意事项 仔细检查爆long long!!加法也会爆 总项数是两个乘数多项式之和的2倍 注意输出的项数！ 注意len的意义！ ...
fft & ntt
swt0_0
08-22 360
FFT&NTT
浅谈FFT&NTT
weixin_30385925的博客
01-27 316
复数及单位根 复数的定义大概就是：\(i^2=-1\)，其中\(i\)就是虚数单位。 那么，在复数意义下，对于方程： \[ x^n=1 \] 就必定有\(n\)个解，这\(n\)个解的分布一定是在复平面上，以圆点为圆心，半径为\(1\)的圆的\(n\)等分点。 由于欧拉公式： \[ e^{i\theta}=\cos\theta+i\cdot \sin\theta \] 把\(2\pi\)带入： \...
NTT学习笔记（快速数论变换）
tanjunming2020的博客
04-29 3588
在$FFT$中，我们选择$n$次单位复根作为$x$的值，因为他们满足消去引理、折半引理、求和引理，所以可以用分治的方法将时间复杂度变为$O(n\log n)$。但因为用了复数，所以在精度上有误差。而在一些题目要求的是带模的多项式乘法，所以$FFT$就用不了了。那么，我们能不能让$x$取整数值，使得这个值也能满足消去引理、折半引理、求和引理呢？
fft/ntt
fighting
03-18 817
题bzoj3160 题意：在一个仅仅含有a,b的字符串里选取一个子序列，使得： 1.位置和字符都关于某条对称轴对称； 2.不能是连续的一段。 分析： 对于每一条对称轴只要求在对称轴两边的对应位置上有相同字符的位置有多少个。 答案就是2len−1−len2^{len}-1-len(减去空串和连续的串) 但这东东它可以是不连续的….. 我们可以将 a b a a b a 1 0
（原理篇）FFT与NTT （快速傅里叶变换）摘自oiwi，仅方便理解
weixin_45666296的博客
08-29 2061
来自oiwi整理加上自己理解，采用jls板子
FFT&NTT数学解释
weixin_30412577的博客
06-01 346
FFT和NTT真是噩梦呢 既然被FFT和NTT坑够了，坑一下其他的人也未尝不可呢 前置知识 多项式基础知识 矩阵基础知识（之后会一直用矩阵表达） FFT：复数基础知识 NTT：模运算基础知识 单位根介绍 设有一个数a，使得an=1，其中n为满足an=1的最小正整数 满足条件的a有哪些呢？ 复数域上的(cos(2π/n)+sin(2π/n)*i)（一般用ωn表...
FFT/NTT高精度
最新发布
08-22
对于使用FFT（快速傅里叶变换）或NTT（快速数论变换）来实现高精度计算，可以采以下步骤： 1. 将数字表示为多项式形式：将需要进行高精度计算的数字转换为多项式的系数形式。例如，将一个整数表示为多项式的系数。 ...
评论
添加红包

请填写红包祝福语或标题

个

红包个数最小为10个

元

红包金额最低5元

当前余额3.43元 前往充值 >
需支付：10.00元
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付元
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明：

1.余额是钱包充值的虚拟货币，按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载，可以购买VIP、付费专栏及课程。

余额充值