多项式模板

最新推荐文章于 2021-02-04 10:35:40 发布
最新推荐文章于 2021-02-04 10:35:40 发布
阅读量124 收藏
点赞数
原文链接：http://www.cnblogs.com/izcat/p/11398764.html
版权

 // 洛谷多项式模板题，整理了自己的数论FFT，NTT的模板。

 

P3803 【模板】多项式乘法（FFT）

#include<cstdio>
#include<iostream>
#include<complex>
#include<cmath>
const double pi = acos(-1);
// #define cp complex<double>
using namespace std;

const int maxn = 1500100;

struct cp {
    double x, y;
    cp(double xx=0, double yy=0):x(xx), y(yy) {}
    cp operator+(const cp& b) {
        return cp(x+b.x, y+b.y);
    }
    cp operator-(const cp& b) {
        return cp(x-b.x, y-b.y);
    }
    cp operator*(const cp& b) {
        return cp(x*b.x-y*b.y, x*b.y+y*b.x);
    }
    cp operator*=(const cp& b) {
        double tmpx = x*b.x - y*b.y;
        y = x*b.y + y*b.x;
        x = tmpx;
        return *this;
    }
};


int n, m;
cp a[maxn*4], b[maxn*4];   // 一定要4倍内存！！！
int rev[maxn*4];

void FFT(cp a[], int n, int inv) {
    // if(!n) return;
    for(int i=0;i<n;i++) {
        // rev[i] = (rev[i>>1]>>1)|((i&1)<<(bit-1));
        if(i<rev[i])                //不加这条if会交换两次（就是没交换）
            swap(a[i], a[rev[i]]);          
    }
    for(int mid=1;mid<n;mid*=2) {   //mid是准备合并序列的长度的二分之一 
        cp temp(cos(pi/mid), inv*sin(pi/mid));
        for(int i=0;i<n;i+=mid*2) { //mid*2是准备合并序列的长度，i是合并到了哪一位
            cp w(1, 0);
            for(int j=0;j<mid;j++,w*=temp) {//只扫左半部分，得到右半部分的答案
                cp Nx = a[i+j];
                cp Ny = w*a[i+j+mid];
                a[i+j] = Nx + Ny;
                a[i+j+mid] = Nx - Ny;
            }
        }
    }
}


int main() {
    scanf("%d %d", &n, &m);
    for(int i=0;i<=n;i++) {
        scanf("%lf", &a[i].x);
    }
    for(int i=0;i<=m;i++) {
        scanf("%lf", &b[i].x);
    }

    int len = 1, bit = 0;
    while(len<=n+m) len <<= 1, ++bit;
    for(int i=0;i<len;i++) 
        rev[i] = (rev[i>>1]>>1)|((i&1)<<(bit-1)); 

    FFT(a, len, 1);
    FFT(b, len, 1);
    for(int i=0;i<=len;i++)
        a[i] *= b[i];
    FFT(a, len, -1);

    for(int i=0;i<=n+m;i++) 
        printf("%d ", (int)(a[i].x/len + 0.5));
    return 0;
}

 

 

P1919 【模板】A*B Problem升级版

利用上面的FFT模板，可以用来做大数乘法。

朴素的模拟手算乘法时间复杂度O(n^2)，n为A，B的位数。改用FFT时间复杂度将优化到 O(nlogn)。

#include<cstdio>
#include<iostream>
#include<cmath>
const double pi = acos(-1);
using namespace std;

const int maxn = 60100;
int n;
struct cp { // 复数类定义
    double x, y;
    cp(double xx=0, double yy=0):x(xx), y(yy) {}
    cp operator+(const cp& b) {
        return cp(x+b.x, y+b.y);
    }
    cp operator-(const cp& b) {
        return cp(x-b.x, y-b.y);
    }
    cp operator*(const cp& b) {
        return cp(x*b.x-y*b.y, x*b.y+y*b.x);
    }
    cp operator*=(const cp& b) {
        double tmpx = x*b.x - y*b.y;
        y = x*b.y + y*b.x;
        x = tmpx;
        return *this;
    }
};

// FFT模板
int rev[maxn*4];   
void FFT(cp a[], int n, int inv) {
    // if(!n) return;
    for(int i=0;i<n;i++) {
        if(i<rev[i])             
            swap(a[i], a[rev[i]]);          
    }
    for(int mid=1;mid<n;mid*=2) { 
        cp temp(cos(pi/mid), inv*sin(pi/mid));
        for(int i=0;i<n;i+=mid*2) { 
            cp w(1, 0);
            for(int j=0;j<mid;j++,w*=temp) {
                cp Nx = a[i+j];
                cp Ny = w*a[i+j+mid];
                a[i+j] = Nx + Ny;
                a[i+j+mid] = Nx - Ny;
            }
        }
    }
}
cp a[maxn*4], b[maxn*4];   // 一定要4倍内存！！！
int ans[maxn*4];
int main() {
    scanf("%d", &n);
    for(int i=n-1;i>=0;i--) {
        scanf("%1lf", &a[i].x);
    }
    for(int i=n-1;i>=0;i--) {
        scanf("%1lf", &b[i].x);
    }

    int len = 1, bit = 0;    
    n <<= 1;
    while(len<=n) len <<= 1, ++bit;
    for(int i=0;i<len;i++) 
        rev[i] = (rev[i>>1]>>1)|((i&1)<<(bit-1)); 

    FFT(a, len, 1);
    FFT(b, len, 1);
    for(int i=0;i<len;i++)
        a[i] *= b[i];
    FFT(a, len, -1);
    
    int carry = 0;
    for(int i=0;i<len;i++) {
        int now = (int)(a[i].x/len + 0.5);
        ans[i] = (carry + now) % 10;
        carry = (carry + now) / 10;
    }
    int s = len-1;
    while(s>=0 && ans[s]==0) --s;
    for(int i=s;i>=0;i--)
        printf("%d", ans[i]);
    puts("");
    return 0;
}
View Code

 

 

 

P4238 【模板】多项式求逆

题意：

给一个多项式 a[x] = ∑ ai*x^i ，求 F[x] 满足 a[x] * F[x] = 1 (mod x^n) 。

求解原理：

假设我们找到 a[x] * G[x] = 1 (mod x^(n/2)) 时的解 G[x]

又由于　　a[x] * F[x] = 1

所以　　　a[x] * ( F[x] - G[x] ) = 0 (mod x^(n/2))

则　　　　F[x] - G[x]  = 0 (mod x^(n/2))

平方后，得

　　　　　F[x] * F[x] + G[x] * G[x] - 2 * F[x] * G[x]  = 0 (mod x^n)

再乘上 a[x]

　　　　　F[x] + a[x] * G[x] * G[x]  - 2 * G[x] = 0 (mod x^n)

所以

　　　　　F[x] = 2 * G[x] - a[x] * G[x] * G[x]

我们找到了递推式，当 n = 1 时，G[0] 为 a[0] 的逆元。

模板适用条件：

模数 mod = 998244353 等原根为 3 的素数。

（代码模板见下一题。）

 

 

 

 

P4725 【模板】多项式对数函数（多项式 ln）

原理：

ln(A[x]) ==> B[x]

两边求导，得

　　　　　A'[x] / A[x] = B'[x]

积分，得

　　　　　B[x] = ∫ ( A'[x] / A[x] ) 

其中 1 / A[x] 可以用前面的多项式求逆模板。

所以

　　　　　B[x] = Interx( A'[x] / A[x] ) = Interx( Dx(A[x]) * Inv(A[x]) )

这两题均使用了NTT和多项式求逆，下面的代码为P4725，P4238将main函数里的调用POLYMONY::Inv即可。

 代码：

#include<cstdio>
#include<iostream>
#include<cmath>
#include<cstring>
using namespace std;
typedef long long ll;
const int maxn = 400100;  // 多项式相乘时需要开4倍大小

namespace POLYNOME {
    const int mod = 998244353;
    const int g = 3;    // 原根： 3 为 mod 的原根
    int rev[maxn];

    int myPow(ll a, int n) {
        ll res = 1;
        while(n) {
            if(n&1) res = res*a % mod;
            a = a*a % mod;
            n >>= 1;
        }
        return res;
    }

    // NTT模板
    // inv==1 快速数论变换
    // inv==-1 逆变换   注意最后点值除以长度为系数结果!!! 
    void NTT(ll a[], int n, int inv) {
        // if(!n) return;
        for(int i=0;i<n;i++) {
            // rev[i] = (rev[i>>1]>>1)|((i&1)<<(bit-1));
            if(i<rev[i]) 
                swap(a[i], a[rev[i]]);          
        }
        for(int mid=1;mid<n;mid*=2) { 
            ll tmp = myPow(g, (mod-1)/(mid*2));
            if(inv==-1) tmp = myPow(tmp, mod-2);
            for(int i=0;i<n;i+=mid*2) {
                ll w = 1;
                for(int j=0;j<mid;j++,w=w*tmp%mod) {
                    int x = a[i+j];
                    int y = w * a[i+j+mid] % mod;
                    a[i+j] = (x + y) % mod;
                    a[i+j+mid] = ((x - y)%mod + mod) % mod;
                }
            }
        }
        //  if(inv==-1) {
        //     for(int i=0;i<n;i++) // a[i]/n
        //         a[i] = a[i]* myPow(n, mod-2) % mod;
        // }
    }


    // 多项式求逆
    // Inv(a[x]) = 1/a[x] => G[x]
    ll F[maxn];
    void Inv(ll a[], ll G[], int n) {
        if(n==1) {
            G[0] = myPow(a[0], mod-2);
            return;
        }

        Inv(a, G, n+1>>1);

        ll lim = 1, bit = 0;
        while(lim<(n<<1)) lim<<=1,   ++bit;

        for(int i=1;i<lim;i++)
            rev[i] = (rev[i>>1]>>1)|((i&1)<<(bit-1)); 

        memset(F, 0, sizeof(0));
        for(int i=0;i<n;i++)
            F[i] = a[i];

        NTT(G, lim, 1);
        NTT(F, lim, 1);
        for(int i=0;i<lim;i++)
            G[i] = (2- F[i]*G[i]%mod + mod)%mod * G[i] % mod;

        NTT(G, lim, -1);
        for(int i=0;i<lim;i++) // G[i]/lim
            G[i] = G[i]* myPow(lim, mod-2) % mod;

        for(int i=n;i<lim;i++) G[i] = 0;
    } 
    
    // 多项式求导
    // Dx(A[x]) => B[x]
    ll Da[maxn];
    void Dx(ll A[], ll B[], int n) {
        for(int i=1;i<n;i++)
            B[i-1] = A[i]*i % mod;
        B[n-1] = 0;
    }

    // 多项式积分
    // Inter(A[x]) => B[x]
    ll Sa[maxn];
    void Interx(ll A[], ll B[], int n) {
        for(int i=1;i<n;i++)
            B[i] = A[i-1]*myPow(i, mod-2) % mod;
        B[0] = 0;
    }

    // 多项式求对数
    // ln(A[x]) => B[x]
    // B[x] = Interx(A'(x)/A(x)) = Interx(Dx(A[x])*Inv(A[x]))
    //                                         Da * Sa
    void lnx(ll a[], ll res[], int n) {
        Inv(a, Sa, n);
        Dx(a, Da, n);

        int lim = 1;
        int bit = 0;
        while(lim<(n<<1)) lim<<=1, ++bit;
        for(int i=0;i<lim;i++)
            rev[i] = (rev[i>>1]>>1)|((i&1)<<(bit-1)); 

        NTT(Da, lim, 1);
        NTT(Sa, lim, 1);
        
        for(int i=0;i<lim;i++)
            Da[i] = Da[i] * Sa[i] % mod;
        NTT(Da, lim, -1);
        for(int i=0;i<lim;i++) // Da[i]/lim
            Da[i] = Da[i]* myPow(lim, mod-2) % mod;

        Interx(Da, res, n);
    }
}


int n;
ll a[maxn], ans[maxn];

int main() {
    scanf("%d", &n);
    for(int i=0;i<n;i++) {
        scanf("%lld", &a[i]);
    }

    POLYNOME::lnx(a, ans, n);

    for(int i=0;i<n;i++) 
        printf("%lld ", ans[i]);
    return 0;
}

 

转载于:https://www.cnblogs.com/izcat/p/11398764.html

weixin_30879833
关注
【多项式n连】各种多项式模板（从入门到入土）
shiveringkonnyaku的博客
03-04 1706
qwq。。。好像很久没有写博客了的样子。。。因为最近一段时间一直处在觉得自己太蒻了的焦虑当中。。。导致啥也没学进去。。。这样的。。。本来定下的肝高数和概率论的计划也因为各种各样的原因没能完成。。。而降低到了只是零散地又复习了一些数学。。。然后发现luogu上有很多多项式的板子。。。就去做了一点。。。这样子。。。然后发现好久没更博客了。。。就准备放上来一些。。。也算是水一水博客了。。。反正老年选手已...
【学习总结】生成函数题目，多项式模板
执子之手,与子偕老
04-22 764
学习资料： 王乐平、策爷冬令营讲义。 多项式技巧 牛顿迭代（泰勒展开） 这是一篇非常好的博客。特别是把所有多项式操作都用泰勒展开推导，以后就不用再死记，也不用担心推错了！from yyb 复合逆 拉格朗日反演 注意当f(x)本身不存在逆元的时候，可以求f(x) / x的逆元 否则（x / f（x））n 的第n - 1 项为0 例题 直接推生成函数 【BZOJ3625】小朋友和二叉树 from ...
多项式模板&题目整理
banchangyue3928的博客
11-27 260
注：多项式的题目，数组应开：N的最近2的整数次幂的4倍。 多项式乘法 FFT模板 时间复杂度\(O(n\log n)\)。 模板： void FFT(Z *a,int x,int K){ static int rev[N],lst; int n=(1<<x); if(n!=lst){ for(int i=0;i<n;i++)rev...
大景深系统多项式相位模板的优化设计
01-11
本篇论文的作者王晋疆、韩瑞雨以及莫绪涛三人，提出了一种基于波前编码的大景深成像系统优化设计，该设计基于立方相位模板，并进一步提出了一种新的奇对称多项式相位模板，以改善成像系统在不同离焦条件下的成像性能...
将不同的边缘检测算子的实验结果进行比较，如下图所示：从图中可以看出，Sobel算子和Prewitt算子能够检测出较为清晰的边缘线，而Roberts算子和Laplacian算子的效果相对较差。这是因为Sobel算子和Prewitt算子采用的是多项式模板，能够更加准确地计算图像梯度，而Roberts算子和Laplacian算子采用的是线性模板，对于图像中的噪声和平滑区域容易产生误检。
最新发布
05-30
相比之下，Sobel算子和Prewitt算子的效果确实要比Roberts算子和Laplacian算子更好，这也验证了多项式模板的优越性。另外，对于Laplacian算子的效果相对较差，我们可以使用高斯滤波器来平滑图像，以提高Laplacian算子...
吉林大学ACM模板.pdf
07-12
目录 目录 .............................................. 1 Graph 图论 ........................................ 3 | DAG的深度优先搜索标记 ................................................
多项式总结及模板
caoyang1123的博客
03-30 490
1.求逆 设要对多项式A求逆，逆为B： 求在模ceil（n/2）的逆，将原式与之相减后平方，发现可得在模n意义下也是0了，然后两边同乘A，移项即可倍增了（注意此时A也是在模当前长度意义下的）。 最后倍增形式为B'=2*B-B*B*A （b’等于2b减b方a），边界为常数取逆元，也可发现多项式有无逆取决于常数有无逆元 2.开根 与求逆类似，先求ceil（n/2）,然后将A移到左边使等式右...
c++模板实现多项式相加
03-03
用c++中的模板实现了链表多项式相加，做的不是很好。仅供参考
多项式模板QAQ
cgh_Andy's Blog
03-06 438
失踪人口回归？ 似乎抄了挺多天的多项式题目代码。。留几个模板吧。。 （都不知道是不是好模板。。 反正也是自己看 最长不超过80行我去 UOJ #34. 多项式乘法 #include &lt;bits/stdc++.h&gt; #define me(a,x) memset(a,x,sizeof a) using namespace std; const int N=3e5+2,inf=1...
多项式系列模板
Never give in.
04-23 719
FFT struct Z { double x,y; Z(double _x=0,double _y=0):x(_x),y(_y){}; }; Z operator +(Z x,Z y) {return Z(x.x+y.x,x.y+y.y);} Z operator -(Z x,Z y) {return Z(x.x-y.x,x.y-y.y);} Z operator *(Z ...
多项式模板（inv, sqrt, ln, exp）
lifehappy
02-04 1376
#include <bits/stdc++.h> using namespace std; const int mod = 998244353, inv2 = mod + 1 >> 1; namespace Quadratic_residue { struct Complex { int r, i; Complex(int _r = 0, int _i = 0) : r(_r), i(_i) {} }; int I2; Complex op
NTT和多项式各种基础运算模板
wzf_2000的博客
06-03 686
#include #define mod 998244353 #define N 400009 #define ll long long using namespace std; ll w[2][N],inv[N]; ll A[N],B[N],C[N],D[N],E[N],F[N]; ll W,G; ll pow(ll a,ll b,ll c) { ll ans=1; for (;b;a=a
【模板】多项式全家桶
weixin_33912638的博客
03-21 145
FFT struct Z { double x,y; Z(double _x=0,double _y=0):x(_x),y(_y){}; }; Z operator +(Z x,Z y) {return Z(x.x+y.x,x.y+y.y);} Z operator -(Z x,Z y) {return Z(x.x-y.x,x.y-y.y);} Z operator *(Z x,Z...
[多项式ln][多项式exp][多项式求幂][生成函数][DP][FNT] BZOJ 3684: 大朋友和多叉树
Vectorxj的博客
03-06 806
SolutionSolutionSolution 把DP写成生成函数的形式。f(x)=x+∑d∈Dfd(x)f(x)=x+∑d∈Dfd(x)f(x)=x+\sum_{d\in D}f^d(x)设g(f(x))=xg(f(x))=xg\big(f(x)\big)=x，有g(f(x))g(x)==f(x)−∑d∈Dfd(x)x−∑d∈Dxdg(f(x))=f(x)−∑d∈Dfd(x)g(x)=x−∑...
【系列】【模板】 多项式
weixin_30505225的博客
12-11 82
【模板】 多项式乘法 学习一波NTT 放弃FFT 1 #include<iostream> 2 #include<cstdio> 3 #include<cstdlib> 4 #include<cmath> 5 #include<algorithm> 6 #include<cstring> ...
评论
添加红包

请填写红包祝福语或标题

个

红包个数最小为10个

元

红包金额最低5元

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

抵扣说明：

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

余额充值