精华内容
下载资源
问答
  • 多项式乘法
    2021-03-03 11:14:36

    给定一个多项式 \(F(x)\) ,请求出一个多项式 \(G(x)\), 满足 \(F(x) * G(x) \equiv 1 ( \mathrm{mod\:} x^n )\)。系数对 \(998244353\) 取模。

    考虑倍增

    多项式只有一项时就是乘法逆元

    假设我们现在得到了

    \[G^{'}(x) \equiv F(x) (mod\ x^{\frac{n}{2}})

    \]

    我们需要求

    \[G(x)\equiv F(x) (mod\ x^n)

    \]

    很明显

    \[G(x)-G^{'}(x)\equiv 0 (mod\ x^{\frac{n}{2}})

    \]

    两边同时平方

    \[G(x)^2-2*G(x)G^{'}(x)+G^{'}(x)\equiv 0 (mod\ x^n)

    \]

    两边同时乘\(F(x)\)得到

    \[F(x)(G(x)^2-2*G(x)G^{'}(x)+G^{'}(x))\equiv 0 (mod\ x^n)

    \]

    由于

    \[F(x)G(x)\equiv 1(mod\ x^n)

    \]

    所以

    \[G(x)-2G(x)^{'}+F(x)G^{'2}(x)\equiv (mod\ x^n)

    \]

    移下项

    \[G(x)\equiv 2G(x)^{'}-F(x)G^{'2}(x)(mod\ x^n)

    \]

    只要初始设置\(G(0)=F(0)^{-1}\)就可以倍增啦~

    额额额不知道有没有人和我有一样的疑惑:为啥递归一次多项式次数少一半,原因是模掉了……

    #include

    using namespace std;

    namespace red{

    #define int long long

    #define eps (1e-8)

    inline int read()

    {

    int x=0;char ch,f=1;

    for(ch=getchar();(ch'9')&&ch!='-';ch=getchar());

    if(ch=='-') f=0,ch=getchar();

    while(ch>='0'&&ch<='9'){x=(x<<1)+(x<<3)+ch-'0';ch=getchar();}

    return f?x:-x;

    }

    const int N=5e5+10,p=998244353,g=3,gi=332748118;

    int n;

    int a[N],b[N],c[N];

    int pos[N];

    inline int fast(int x,int k)

    {

    int ret=1;

    while(k)

    {

    if(k&1) ret=ret*x%p;

    x=x*x%p;

    k>>=1;

    }

    return ret;

    }

    inline void ntt(int limit,int *a,int inv)

    {

    for(int i=0;i

    if(i

    for(int mid=1;mid

    {

    int Wn=fast(g,(p-1)/(mid<<1));

    for(int r=mid<<1,j=0;j

    {

    int w=1;

    for(int k=0;k

    {

    int x=a[j+k],y=w*a[j+k+mid]%p;

    a[j+k]=x+y;

    if(a[j+k]>=p) a[j+k]-=p;

    a[j+k+mid]=x-y;

    if(a[j+k+mid]<0) a[j+k+mid]+=p;

    }

    }

    }

    if(inv) return;

    inv=fast(limit,p-2);reverse(a+1,a+limit);

    for(int i=0;i

    }

    inline void poly_inv(int pw,int *a,int *b)

    {

    if(pw==1) {b[0]=fast(a[0],p-2);return;}

    poly_inv((pw+1)>>1,a,b);

    int len=0,limit=1;

    while(limit

    for(int i=1;i>1]>>1)|((i&1)<

    for(int i=0;i

    for(int i=pw;i

    ntt(limit,c,1);ntt(limit,b,1);

    for(int i=0;i

    ntt(limit,b,0);

    for(int i=pw;i

    }

    inline void main()

    {

    n=read();

    for(int i=0;i

    poly_inv(n,a,b);

    for(int i=0;i

    }

    }

    signed main()

    {

    red::main();

    return 0;

    }

    更多相关内容
  • (2)设计算法实现一元多项式乘法; (3)分析算法的时间复杂度和空间复杂度 一、总体设计 1 二、详细设计 1 2.1存储结构 1 2.2建立链表 1 2.3遍历操作 1 2.4多项式相乘算法 2 三、调试与测试 2 3.1方案一 2 3.2...
  • 基于Ring-LWE的密码系统最关键且计算量大的操作是环上的多项式乘法。 在本文中,我们介绍了几种优化技术... 借助这些技术,我们的多项式乘法器能够在Spartan-6 FPGA上针对维度256/512每秒执行57304/26913多项式乘法
  • ifft实现代码 matlab
  • 多项式乘法实现

    2017-12-22 16:48:14
    数据结构作业,适合初学者,小白,自学者,和大家共享。
  • 今天小编就为大家分享一篇关于Java实现多项式乘法代码实例,小编觉得内容挺不错的,现在分享给大家,具有很好的参考价值,需要的朋友一起跟随小编来看看吧
  • 为降低基于重模多项式剩余类环矩阵的密码算法中乘法运算占用的硬件资源量,提出了一种剩余类环上乘法的流水线实现方法。该方法选用数模为216,多项式模为4次首一多项式的重模多项式剩余类环,对流水线设计进行了数学...
  • 多项式乘法

    2022-03-01 21:27:47
    多项式乘法问题: 给定两个多项式: A(x) = a0 + a1… + anXn B(x) = b0 + b1… + bnXn 求A(x)B(x): 假设A,B的系数ai,bi存储在数组A[0…n],B[0…m]中。 第一个分治算法法: 分:四次乘法 B(x)也以相同的方法分。...

    多项式乘法问题:
    给定两个多项式:
    A(x) = a0 + a1… + anXn
    B(x) = b0 + b1… + bnXn
    求A(x)B(x):

    假设A,B的系数ai,bi存储在数组A[0…n],B[0…m]中。
    第一个分治算法法
    四次乘法
    在这里插入图片描述
    B(x)也以相同的方法分。
    在这里插入图片描述
    把原问题(大小为n)分为4个输入大小为n/2的问题

    治:解决这四个子问题
    通过递归调用算法4次
    在这里插入图片描述
    合;添加以下四个多项式
    在这里插入图片描述
    时间复杂度
    在这里插入图片描述
    和暴力解法的时间复杂度一样,并不是最好的分治法。

    改进的分治算法
    分:三次乘法
    在这里插入图片描述
    在这里插入图片描述
    时间复杂度
    在这里插入图片描述在这里插入图片描述

    展开全文
  • 还使用多项式乘法算法来测试基本的DFT算法。 在O(n log n)多项式乘法中的应用 DFT算法可以方便地应用于多项式乘法。 结果是O(n log n)多项式乘法算法。 实现是。 此应用程序还用于测试FFT的实现,测试用例为。 ...
  • 一元多项式乘法算法

    2019-03-18 02:37:38
    NULL 博文链接:https://qingdaoguy.iteye.com/blog/694445
  • % 使用多项式乘法进行计算机双线性变换% 输入是 H(s) 降序多项式% 输出是降序多项式中的 H(z) % % 用法:[c,d] = mBilinear(a,b,Fs); % with a: 按 H(s) 分子降序排列的向量% b:按 H(s) 分子降序排列的向量% c: H(z...
  • 多项式乘法(FFT)

    万次阅读 多人点赞 2018-11-21 19:10:04
    多项式乘法(FFT) 前言 作为一名OI选手,至今未写过fft相关的博客,真是一大遗憾,这也导致我并没有真正推过fft的所有式子 这一篇fft的博客我将详细介绍多项式乘法,易于理解,当然,一些公式会有些枯燥,如果是...

    1 前言

    作为一名OI选手,至今未写过fft相关的博客,真是一大遗憾,这也导致我并没有真正推过fft的所有式子
    这一篇fft的博客我将详细介绍多项式乘法,易于理解,主要是为了等我啥时候忘了回来看,当然,一些公式会有些枯燥,如果是初学者请耐心看完哦,还有,毕竟这是手写出来的,如果有错误,欢迎指正!

    2 介绍

    本栏用来普及一些知识和对FFT的思路进行描述
    多项式乘法,顾名思义,首先是讲到多项式,那么什么是多项式呢?

    2.1 多项式

    首先是多项式的定义,想必大家都知道(你上过初中吧),而在这里,我们所说的多项式都是单个未知数x的
    所以,在我们正常人眼中的一个次多项式就是形如 f ( x ) = a 0 x 0 + a 1 x 1 + ⋅ ⋅ ⋅ + a n − 1 x n − 1 f(x)=a_0x^0+a_1x^1+···+a_{n-1}x^{n-1} f(x)=a0x0+a1x1++an1xn1
    没错,这就是大名鼎鼎的系数表示法
    然后呢,由于在后面要用到,所以我在这里再介绍一种点值表示法
    就是将n个不同的值 x 0 , x 1 ⋅ ⋅ ⋅ x n − 1 x_0,x_1···x_{n-1} x0,x1xn1分别带入 f ( x ) f(x) f(x),获得n个结果 y 0 , y 1 ⋅ ⋅ ⋅ y n − 1 y_0,y_1···y_{n-1} y0,y1yn1,这n对数 ( x 0 , y 0 ) , ( x 1 , y 1 ) ⋅ ⋅ ⋅ ( x n − 1 , y n − 1 ) (x_0,y_0),(x_1,y_1)···(x_{n-1},y_{n-1}) (x0,y0),(x1,y1)(xn1,yn1)就可以表示出这个多项式
    看到这里,如果你是初学者,你一定会感到非常迷茫,这为什么对呢?
    看到这里,如果你是个FFT高手,你可能会感到迷茫,这为什么对呢?
    (dalao勿喷)
    在这里插入图片描述
    这张图片里系数表示法相当于是最右侧的那个矩阵
    而点值表示法则包含了左边的两个矩阵,可以通过这两个矩阵计算出最右侧的那个矩阵,所以两种表示法是等价的
    撒花
    注:另外要说的是,由于算法需要,本博客所说的n次多项式都默认n是2的幂次(如果不足可以添加0来补)

    2.2多项式的乘法

    在做了那么久的各种数学题后,我对多项式乘法有了有了的理解
    对于一个一般的给定系数表示法的多项式乘法问题
    比如两个n次多项式A(x),B(x),给出系数表示法,求它们的乘积C(x)
    分别枚举两个多项式中的每一项,分别是 O ( n ) 的 O(n)的 O(n),所以总复杂度为 O ( n 2 ) O(n^2) O(n2)
    这是一个很方便的做法

    你可以发现一件很有趣的事情,那就是如果给出的是点值表示法,并且两个多项式的x分别对应相等,那么把y对应相乘,就能 O ( n ) O(n) O(n)的获取乘积的点值表示法

    2.3 快速傅立叶变换(FFT)

    那么,FFT是用来干什么的呢?
    对于一个多项式乘法问题,当给出系数表示法的时候, O ( n 2 ) O(n^2) O(n2)的复杂度有时候并不足够优越,而FFT就是一个能使多项式乘法做到 O ( n l o g n ) O(nlogn) Onlogn的一个算法,具体的原理其实非常清晰

    • 两个多项式的系数表示法
      求值,O(nlogn)
    • 两个多项式的点值表示法
      点值乘法,O(n)
    • 两个多项式乘积的点值表示法
      插值,O(nlogn)
    • 两个多项式乘积的系数表示法

    是不是一目了然呢?当然,要具体实现,还需要细细说来

    3 实现

    现在你已经大致知道FFT要干什么了,现在你已经会在点值情况下 O ( n ) O(n) O(n)进行多项式乘法,剩下的就是要解决两个问题——求值与插值了

    3.1 暴力算法( O ( n 3 ) O(n^3) O(n3)

    要先做题,必先暴力
    首先是求值,加入你现在随便找了n个互不相同的x,带入其中,是什么复杂度呢 O ( n 2 ) O(n^2) O(n2)
    然后是插值,有一个非常妙的方法,假设所有的a都是未知数,那么这个问题就变成了经典的高斯消元问题,复杂度 O ( n 3 ) O(n^3) O(n3)
    不好意思,这两个操作的复杂度都光荣的在 O ( n 2 ) O(n^2) O(n2)以上,使得当前这个算法的总复杂度为 O ( n 3 ) O(n^3) O(n3)比文章开始的那个 O ( n 2 ) O(n^2) On2都要差,不要灰心,既然复杂度不优,那就循序渐进的优化

    3.2 离散傅里叶变换(通过优化使上面算法复杂度降到 O ( n 2 ) O(n^2) O(n2),请仔细看完,这是基础)

    你会发现,点值表示法有一个很好的特性,就是那个代入的x可以自己选择
    离散傅里叶变换的思路是将n个x的值取n个单位根(模长为一的复数)

    复数(这是一个知识拓展框)

    − 1 \sqrt{-1} 1 这个数,在实数范围内是不存在的,所以拓展出复数这一概念, i = − 1 i=\sqrt{-1} i=1 ,复数就是能够被表示为 z = x + y ∗ i z=x+y*i z=x+yi的数。所以对一个复数,可以用有序数对(x,y)表示,在坐标轴上有对应的点,而这个复数就是从(0,0)到(x,y)的一条有向线段(只会向量的同学可以把它看成向量),而这个复数的模长就等于(0,0)到(x,y)的距离
    由于复数是数,所以也有各种运算
    加法:(a+bi)+(c+di)=(a+c)+(b+d)i
    减法:(a+bi)-(c+di)=(a-c)+(b-d)i
    乘法:(a+bi)*(c+di)=(ac-bd)+(ad+bc)i
    当然,C++有专门的complex变量可以声明,但是

    不推荐使用!!!

    为什么呢?因为FFT本身就有一定的常数,如果再用系统complex常数会更大,所以推荐自己手写struct

    那么什么是单位根呢?

    3.2.1 单位根

    单位根所在的点是把单位圆(以原点为圆心,半径为1的圆)从(0,1)开始平均分成n份的分割点
    如下图,这就是n=8时的单位圆,绿色圆上的红点就是单位根所在的点
    在这里插入图片描述
    从(0,1)开始逆时针将这n个点编号,所表示的单位根分别为 w n 1 , w n 1 ⋅ ⋅ ⋅ , w n n − 1 w_n^1,w_n^1···,w_n^{n-1} wn1,wn1,wnn1,特殊的, w n 1 w_n^1 wn1被称为n次单位根。容易发现每个单位根都非常好算,即 w n k = ( c o s k n 2 π , s i n k n 2 π ) w_n^k=(cos\frac{k}{n} 2π,sin\frac{k}{n} 2π) wnk=(cosnk2π,sinnk2π)
    这个用三角函数的想法非常好证
    知道了这个之后,你会发现很多性质

    性质1: w n k = ( w n 1 ) k w_n^k=(w_n^1)^k wnk=(wn1)k

    证明:
    w n k ∗ w n 1 = ( c o s k n 2 π , s i n k n 2 π ) ∗ ( c o s 1 n 2 π , s i n 1 n 2 π ) = ( c o s k n 2 π ∗ c o s 1 n 2 π − s i n k n 2 π ∗ s i n 1 n 2 π , s i n k n 2 π ∗ c o s 1 n 2 π + c o s k n 2 π ∗ s i n 1 n 2 π ) = ( c o s k + 1 n 2 π , s i n k + 1 n 2 π ) = w n k + 1 \begin{aligned} w_n^k*w_n^1&amp;=(cos\frac{k}{n} 2π,sin\frac{k}{n} 2π)*(cos\frac{1}{n} 2π,sin\frac{1}{n} 2π)\\ &amp;=(cos\frac{k}{n} 2π*cos\frac{1}{n} 2π-sin\frac{k}{n} 2π*sin\frac{1}{n} 2π, sin\frac{k}{n} 2π*cos\frac{1}{n} 2π+cos\frac{k}{n} 2π*sin\frac{1}{n} 2π)\\ &amp;=(cos\frac{k+1}{n} 2π,sin\frac{k+1}{n} 2π)\\ &amp;=w_n^{k+1} \end{aligned} wnkwn1=(cosnk2π,sinnk2π)(cosn12π,sinn12π)=(cosnk2πcosn12πsinnk2πsinn12π,sinnk2πcosn12π+cosnk2πsinn12π)=(cosnk+12π,sinnk+12π)=wnk+1

    如果你想问倒数第二个等号怎么等于过去的,请查看

    和角公式百度链接:https://baike.baidu.com/item/和角公式/8782319?fr=aladdin

    update:lyc建议这里补充一下,我就加一下吧
    众所周知 ( a x ) y = a x ∗ y (a^x)^y=a^{x*y} (ax)y=axy
    那么知道这个之后就有 ( w n x ) y = ( ( w n 1 ) x ) y = ( w n 1 ) x ∗ y = w n x ∗ y (w_n^x)^y=((w_n^1)^x)^y=(w_n^1)^{x*y}=w_n^{x*y} (wnx)y=((wn1)x)y=(wn1)xy=wnxy

    性质2:对于任意一个正整数x, w n ∗ x k ∗ x = w n k w_{n*x}^{k*x}=w_n^k wnxkx=wnk

    证明:
    w n ∗ x k ∗ x = ( c o s k ∗ x n ∗ x 2 π , s i n k ∗ x n ∗ x 2 π ) = ( c o s k n 2 π , s i n k n 2 π ) = w n k \begin{aligned} w_{n*x}^{k*x}&amp;=(cos\frac{k*x}{n*x} 2π,sin\frac{k*x}{n*x} 2π)\\ &amp;=(cos\frac{k}{n} 2π,sin\frac{k}{n} 2π)\\ &amp;=w_n^k \end{aligned} wnxkx=(cosnxkx2π,sinnxkx2π)=(cosnk2π,sinnk2π)=wnk
    没错,约分大法好,这个等式说明,这两个数在单位圆上对应的点是同一个,这个性质,使 w 2 n 2 k = w n k w_{2n}^{2k}=w_n^k w2n2k=wnk

    性质3:如果n是偶数,那么 w n k + n 2 = − w n k w_n^{k+\frac n2}=-w_n^k wnk+2n=wnk

    证明:
    w n k + n 2 = ( c o s k + n 2 n 2 π , s i n k + n 2 n 2 π ) = ( − c o s k n 2 π , − s i n k n 2 π ) = − ( c o s k n 2 π , s i n k n 2 π ) = − w n k \begin{aligned} w_n^{k+\frac n2}&amp;=(cos\frac{k+\frac n2}{n} 2π,sin\frac{k+\frac n2}{n} 2π)\\ &amp;=(-cos\frac{k}{n} 2π,-sin\frac{k}{n} 2π)\\ &amp;=-(cos\frac{k}{n} 2π,sin\frac{k}{n} 2π)\\ &amp;=-w_n^k \end{aligned} wnk+2n=(cosnk+2n2π,sinnk+2n2π)=(cosnk2π,sinnk2π)=(cosnk2π,sinnk2π)=wnk
    诱导公式大法好,
    s i n ( π + α ) = - s i n α sin(π+α)= -sinα sinπ+α=sinα
    c o s ( π + α ) = - c o s α cos(π+α)=-cosα cosπ+α=cosα
    理解一下,相当于这两者是单位圆上相对的两个点,值自然是取相反数的啦

    3.2.2代入单位根带来的性质

    你现在已经知道单位根是什么啦
    那么,我们回头看这个离散傅里叶变换,它是求值的时候把x的值分别取 w n 0 , w n 1 ⋅ ⋅ ⋅ w n n − 1 w_n^0,w_n^1···w_n^{n-1} wn0,wn1wnn1这n个数,究竟这么做有什么好处呢?
    答案是——你可以比较方便的实现插值!!!
    哇塞,这真是很牛逼的呢,插值是暴力算法的瓶颈,如果能优化,那就可以优化总复杂度了
    那么如何优化呢?

    现在定义对函数f(x)的 离散傅里叶变换为将 w n 0 , w n 1 ⋅ ⋅ ⋅ w n n − 1 w_n^0,w_n^1···w_n^{n-1} wn0,wn1wnn1这n个数作为 x 0 , x 1 ⋅ ⋅ ⋅ x n − 1 x_0,x_1···x_{n-1} x0,x1xn1代入, 离散傅里叶变换的结果为 ( y 0 , y 1 ⋅ ⋅ ⋅ y n − 1 ) (y_0,y_1···y_{n-1}) (y0,y1yn1),容易发现,这是一个插值的过程
    然后有一个结论:
    一个多项式A(x)在进行离散傅里叶变换后,将离散傅里叶变换的结果的n个y作为系数组成多项式B(x),原来的n个单位根取倒数进行求值,结果的每个数除以n,其结果就是A(x)的各项系数
    文字说的可能不太清晰,用数字来表达就是这样的:


    A ( x ) = a 0 x 0 + a 1 x 1 + ⋅ ⋅ ⋅ + a n − 1 x n − 1 A(x)=a_0x^0+a_1x^1+···+a_{n-1}x^{n-1} A(x)=a0x0+a1x1++an1xn1
    w n 0 , w n 1 ⋅ ⋅ ⋅ w n n − 1 w_n^0,w_n^1···w_n^{n-1} wn0,wn1wnn1作为x分别带入求值
    得到 ( y 0 , y 1 ⋅ ⋅ ⋅ y n − 1 ) (y_0,y_1···y_{n-1}) (y0,y1yn1)
    将这些y作为系数,产生一个新的多项式B(x)
    B ( x ) = y 0 x 0 + y 1 x 1 + ⋅ ⋅ ⋅ + y n − 1 x n − 1 B(x)=y_0x^0+y_1x^1+···+y_{n-1}x^{n-1} B(x)=y0x0+y1x1++yn1xn1
    w n − 0 , w n − 1 ⋅ ⋅ ⋅ w n − ( n − 1 ) w_n^{-0},w_n^{-1}···w_n^{-(n-1)} wn0,wn1wn(n1)作为x分别带入求值
    得到的 ( z 0 , z 1 ⋅ ⋅ ⋅ z n − 1 ) (z_0,z_1···z_{n-1}) (z0,z1zn1)
    对于每个 z k z_k zk,有 z k = a k ∗ n z_k=a_k*n zk=akn

    证明

    z k = ∑ i = 0 n − 1 y i ∗ ( w n − k ) i = ∑ i = 0 n − 1 ( ∑ j = 0 n − 1 a j ∗ ( w n i ) j ) ∗ ( w n − k ) i = ∑ i = 0 n − 1 ∑ j = 0 n − 1 a j ∗ ( w n i ) j ∗ ( w n − k ) i = ∑ j = 0 n − 1 ∑ i = 0 n − 1 a j ∗ w n i ∗ j ∗ w n − k ∗ i = ∑ j = 0 n − 1 a j ( ∑ i = 0 n − 1 ( w n j − k ) i ) \begin{aligned} z_k&amp;=\sum_{i=0}^{n-1}y_i*(w_n^{-k})^i\\ &amp;=\sum_{i=0}^{n-1}(\sum_{j=0}^{n-1}a_j*(w_n^i)^j)*(w_n^{-k})^i\\ &amp;=\sum_{i=0}^{n-1}\sum_{j=0}^{n-1}a_j*(w_n^i)^j*(w_n^{-k})^i\\ &amp;=\sum_{j=0}^{n-1}\sum_{i=0}^{n-1}a_j*w_n^{i*j}*w_n^{-k*i}\\ &amp;=\sum_{j=0}^{n-1}a_j(\sum_{i=0}^{n-1}(w_n^{j-k})^i) \end{aligned} zk=i=0n1yi(wnk)i=i=0n1(j=0n1aj(wni)j)(wnk)i=i=0n1j=0n1aj(wni)j(wnk)i=j=0n1i=0n1ajwnijwnki=j=0n1aj(i=0n1(wnjk)i)
    然后对于 ∑ i = 0 n − 1 ( w n j − k ) i \sum_{i=0}^{n-1}(w_n^{j-k})^i i=0n1(wnjk)i,容易发现
    如果 j − k = 0 j-k= 0 jk=0那么 w n j − k w_n^{j-k} wnjk就是 w n 0 = 1 w_n^0=1 wn0=1,所以 n n n 1 1 1结果为 n n n
    若果 j − k ≠ 0 j-k\neq0 jk̸=0那么通过等比数列求和( x 0 + x 1 + ⋅ ⋅ ⋅ + x n − 1 = x n − 1 x − 1 x^0+x^1+···+x^{n-1}=\frac {x^n-1}{x-1} x0+x1++xn1=x1xn1)可以发现,结果 ( w n j − k ) n − 1 w n j − k − 1 = ( w n n ) j − k − 1 w n j − k − 1 = 1 j − k − 1 w n j − k − 1 = 0 \frac {(w_n^{j-k})^n-1}{w_n^{j-k}-1}=\frac {(w_n^n)^{j-k}-1}{w_n^{j-k}-1}=\frac {1^{j-k}-1}{w_n^{j-k}-1}=0 wnjk1(wnjk)n1=wnjk1(wnn)jk1=wnjk11jk1=0
    所以说,这个系数只有在 j − k = 0 j-k=0 jk=0 j = k j=k j=k时才为n,其它都为0
    所以 z k = a k ∗ n z_k=a_k*n zk=akn
    证毕


    对于这个结论,你会发现,如果你用离散傅里叶变换,你的插值就变成了一次求值,你现在的瓶颈也就变成了只有求值这个操作了,NICE!
    现在暴力带入的求值的复杂度为 O ( n 2 ) O(n^2) O(n2),所以整个算法的复杂度也为 O ( n 2 ) O(n^2) O(n2)

    3.3 快速傅里叶变换(使整个算法复杂度优化到 O ( n l o g n ) O(nlogn) O(nlogn)

    现在的复杂度变成 O ( n 2 ) O(n^2) O(n2)了,你可能会说,这不是和暴力一样的复杂度嘛,学了老半天,还是个大常数 O ( n 2 ) O(n^2) O(n2),真没用
    别着急,现在算法瓶颈在于求值,只要优化它的复杂度,算法就能变优
    然后快速傅里叶变换就来了

    Q:傅里叶就可以为所欲为吗?
    A:没错,傅里叶就是可以为所欲为!
    解释:来一个傅里叶百度百科的链接,他作为一名数学家、物理学家,在计算机发明100+年前就弄出了这个傅里叶变换!!!太巨了orz

    现在要做的是,对于一个多项式 A ( x ) = a 0 x 0 + a 1 x 1 + ⋅ ⋅ ⋅ + a n − 1 x n − 1 A(x)=a_0x^0+a_1x^1+···+a_{n-1}x^{n-1} A(x)=a0x0+a1x1++an1xn1,我们需要快速的获得代入 w n 0 , w n 1 ⋅ ⋅ ⋅ w n n − 1 w_n^0,w_n^1···w_n^{n-1} wn0,wn1wnn1的结果(如果这个ok那么代入 w n − 0 , w n − 1 ⋅ ⋅ ⋅ w n − ( n − 1 ) w_n^{-0},w_n^{-1}···w_n^{-(n-1)} wn0,wn1wn(n1)也行)
    这个快速傅里叶的一个思路就是分治
    首先把这个多项式按次数奇偶分组
    A ( x ) = ( a 0 x 0 + a 2 x 2 + ⋅ ⋅ ⋅ + a n − 2 x n − 2 ) + ( a 1 x 1 + a 3 x 3 + ⋅ ⋅ ⋅ + a n − 1 x n − 1 ) A(x)=(a_0x^0+a_2x^2+···+a_{n-2}x^{n-2})+(a_1x^1+a_3x^3+···+a_{n-1}x^{n-1}) A(x)=(a0x0+a2x2++an2xn2)+(a1x1+a3x3++an1xn1)

    A 1 ( x ) = ( a 0 x 0 + a 2 x 1 + ⋅ ⋅ ⋅ + a n − 2 x n 2 − 1 ) A_1(x)=(a_0x^0+a_2x^1+···+a_{n-2}x^{\frac n2-1}) A1(x)=(a0x0+a2x1++an2x2n1)
    A 2 ( x ) = ( a 1 x 0 + a 3 x 1 + ⋅ ⋅ ⋅ + a n − 1 x n 2 − 1 ) A_2(x)=(a_1x^0+a_3x^1+···+a_{n-1}x^{\frac n2-1}) A2(x)=(a1x0+a3x1++an1x2n1)
    那么有
    A ( x ) = A 1 ( x 2 ) + x ∗ A 2 ( x 2 ) A(x)=A_1(x^2)+x*A_2(x^2) A(x)=A1(x2)+xA2(x2)
    对于所有的值
    如果小于 n 2 \frac n2 2n,那么直接带入这个 k k k,有
    A ( w n k ) = A 1 ( w n 2 k ) + w n k ∗ A 2 ( w n 2 k ) = A 1 ( w n 2 k ) + w n k ∗ A 2 ( w n 2 k ) \begin{aligned} A(w_n^k)&amp;=A_1(w_n^{2k})+w_n^k*A_2(w_n^{2k})\\ &amp;=A_1(w_{\frac n2}^k)+w_n^k*A_2(w_{\frac n2}^k) \end{aligned} A(wnk)=A1(wn2k)+wnkA2(wn2k)=A1(w2nk)+wnkA2(w2nk)
    如果大于等于 n 2 \frac n2 2n,同样带入,用 k + n 2 k+\frac n2 k+2n表示,有
    A ( w n k + n 2 ) = A 1 ( w n 2 k + n ) + w n k + n 2 ∗ A 2 ( w n 2 k + n ) = A 1 ( w n 2 k ∗ w n n ) + w n k + n 2 ∗ A 2 ( w n 2 k ∗ w n n ) = A 1 ( w n 2 k ) − w n k ∗ A 2 ( w n 2 k ) \begin{aligned} A(w_n^{k+\frac n2})&amp;=A_1(w_n^{2k+n})+w_n^{k+\frac n2}*A_2(w_n^{2k+n})\\ &amp;=A_1(w_n^{2k}*w_n^n)+w_n^{k+\frac n2}*A_2(w_n^{2k}*w_n^n)\\ &amp;=A_1(w_n^{2k})-w_n^k*A_2(w_n^{2k}) \end{aligned} A(wnk+2n)=A1(wn2k+n)+wnk+2nA2(wn2k+n)=A1(wn2kwnn)+wnk+2nA2(wn2kwnn)=A1(wn2k)wnkA2(wn2k)
    然后现在如果知道 A 1 ( x ) A_1(x) A1(x) A 2 ( x ) A_2(x) A2(x)代入 w n 2 0 , w n 2 1 ⋅ ⋅ ⋅ w n 2 n 2 − 1 w_{\frac n2}^0,w_{\frac n2}^1···w_{\frac n2}^{\frac n2-1} w2n0,w2n1w2n2n1的值,那 A ( x ) A(x) A(x)代入 w n 0 , w n 1 ⋅ ⋅ ⋅ w n n − 1 w_n^0,w_n^1···w_n^{n-1} wn0,wn1wnn1的值也可以 O ( n ) O(n) O(n)计算出来了,然后递归解决问题
    递归会有终止条件,当n=1的时候带入 w 1 0 w_1^0 w10的值就是那个多项式的 a 0 a_0 a0,就可以直接return了
    考虑时间复杂度的分析 T ( n ) = 2 ∗ T ( n 2 ) + O ( n ) T(n)=2*T(\frac n2)+O(n) T(n)=2T(2n)+O(n)
    总复杂度是 O ( n l o g n ) O(nlogn) O(nlogn)的,完成!
    在FFT说完之际,贴一个经典的揭示暴力多项式乘法和FFT不同的图:在这里插入图片描述
    然后你就会写FFT了,贴一个FFT的递归写法:

    #include<cstdio>
    #include<cctype>
    #include<cmath>
    namespace fast_IO
    {
        const int IN_LEN=10000000,OUT_LEN=10000000;
        char ibuf[IN_LEN],obuf[OUT_LEN],*ih=ibuf+IN_LEN,*oh=obuf,*lastin=ibuf+IN_LEN,*lastout=obuf+OUT_LEN-1;
        inline char getchar_(){return (ih==lastin)&&(lastin=(ih=ibuf)+fread(ibuf,1,IN_LEN,stdin),ih==lastin)?EOF:*ih++;}
        inline void putchar_(const char x){if(oh==lastout)fwrite(obuf,1,oh-obuf,stdout),oh=obuf;*oh++=x;}
        inline void flush(){fwrite(obuf,1,oh-obuf,stdout);}
    }
    using namespace fast_IO;
    #define getchar() getchar_()
    #define putchar(x) putchar_((x))
    typedef long long LL;
    #define rg register
    template <typename T> inline void read(T&x)
    {
        char cu=getchar();x=0;bool fla=0;
        while(!isdigit(cu)){if(cu=='-')fla=1;cu=getchar();}
        while(isdigit(cu))x=x*10+cu-'0',cu=getchar();
        if(fla)x=-x;  
    }
    template <typename T> void printe(const T x)
    {
        if(x>=10)printe(x/10);
        putchar(x%10+'0');
    }
    template <typename T> inline void print(const T x)
    {
        if(x<0)putchar('-'),printe(-x);
        else printe(x);
    }
    const int maxn=2097153;
    const double Pi=acos(-1.0);
    struct complex
    {
        double x,y;
        inline complex operator +(const complex b)const{return (complex){x+b.x,y+b.y};}
        inline complex operator *(const complex b)const{return (complex){x*b.x-y*b.y,x*b.y+y*b.x};}
        inline complex operator -(const complex b)const{return (complex){x-b.x,y-b.y};}
    }a[maxn],b[maxn];
    int n,m,allsum;
    void FFT(int lenth,complex*A,const int fla)
    {
        if(lenth==1)return;
        complex A1[lenth>>1],A2[lenth>>1];
        for(rg int i=0;i<lenth;i+=2)A1[i>>1]=A[i],A2[i>>1]=A[i+1];
        FFT(lenth>>1,A1,fla),FFT(lenth>>1,A2,fla);
        const complex w=(complex){cos(Pi*2.0/lenth),sin(Pi*2.0/lenth)*fla};
        complex k=(complex){1,0};
        lenth>>=1;
        for(rg int i=0;i<lenth;i++,k=k*w)
        {
            A[i]=A1[i]+k*A2[i];
            A[i+lenth]=A1[i]-k*A2[i];
        }
    }
    int main()
    {
        read(n),read(m),allsum=n+m;
        for(rg int i=0;i<=n;i++)read(a[i].x);
        for(rg int i=0;i<=m;i++)read(b[i].x);
        rg int lenth=1;while(lenth<=n+m)lenth<<=1;
        FFT(lenth,a,1),FFT(lenth,b,1);
        for(rg int i=0;i<=lenth;i++)a[i]=a[i]*b[i];
        FFT(lenth,a,-1);
        for(rg int i=0;i<=n+m;i++)print((int)(a[i].x/lenth+0.5)),putchar(' ');
        return flush(),0;
    }
    

    对代码的一些解释:
    那个FFT()函数是用来求值的,前面已经证明过插值就是把单位根取倒数,所以单位根的标号传-1就好了
    另外的部分都是模拟

    4 优化

    4.1 递归转迭代优化

    之前贴的代码在某评测网站上运行最大数据点所花的时间为2493ms,题目数据范围是 n ≤ 1 e 6 n\leq1e6 n1e6,可见这个FFT的速度有一定的常数,这个时候就要考虑一些优化

    首先想到的优化常数的算法自然是把递归转迭代了,而这种优化也是最为常见的

    某一个写FFT的人发现如下性质:对于 a x a_x ax这个数,递归到最后所在的位置刚好是x的二进制位全部翻转的那一位,比如说4的二进制是(100),最后到了1(001),更多的感兴趣可以自己手模
    考虑到你可能连前面的都没看懂,没有能力手模,我还是把它画出来吧

    0,1,2,3,4,5,6,7
    0,2,4,6
    1,3,5,7
    0,4
    2,6
    1,5
    3,7
    0
    4
    2
    6
    1
    5
    3
    7
    000
    100
    010
    110
    001
    101
    011
    111
    000
    001
    010
    011
    100
    101
    110
    111

    你可以自行对比最后两排看看是不是这样
    然后预处理二进制翻转的数组,然后就可以非递归的从底层一层一层往上做了,在我的代码中处理的是Reverse数组
    贴一波代码

    #include<cstdio>
    #include<cctype>
    #include<cmath>
    namespace fast_IO
    {
    	const int IN_LEN=10000000,OUT_LEN=10000000;
    	char ibuf[IN_LEN],obuf[OUT_LEN],*ih=ibuf+IN_LEN,*oh=obuf,*lastin=ibuf+IN_LEN,*lastout=obuf+OUT_LEN-1;
    	inline char getchar_(){return (ih==lastin)&&(lastin=(ih=ibuf)+fread(ibuf,1,IN_LEN,stdin),ih==lastin)?EOF:*ih++;}
    	inline void putchar_(const char x){if(oh==lastout)fwrite(obuf,1,oh-obuf,stdout),oh=obuf;*oh++=x;}
    	inline void flush(){fwrite(obuf,1,oh-obuf,stdout);}
    }
    using namespace fast_IO;
    #define getchar() getchar_()
    #define putchar(x) putchar_((x))
    typedef long long LL;
    #define rg register
    template <typename T> inline void swap(T&a,T&b){T c=a;a=b;b=c;}
    template <typename T> inline void read(T&x)
    {
        char cu=getchar();x=0;bool fla=0;
        while(!isdigit(cu)){if(cu=='-')fla=1;cu=getchar();}
        while(isdigit(cu))x=x*10+cu-'0',cu=getchar();
        if(fla)x=-x;  
    }
    template <typename T> void printe(const T x)
    {
        if(x>=10)printe(x/10);
        putchar(x%10+'0');
    }
    template <typename T> inline void print(const T x)
    {
        if(x<0)putchar('-'),printe(-x);
        else printe(x);
    }
    const int maxn=2097153;const double PI=acos(-1.0);
    int n,m;
    struct complex
    {
    	double x,y;
    	inline complex operator +(const complex b)const{return (complex){x+b.x,y+b.y};}
    	inline complex operator -(const complex b)const{return (complex){x-b.x,y-b.y};} 
    	inline complex operator *(const complex b)const{return (complex){x*b.x-y*b.y,x*b.y+y*b.x};}
    }a[maxn],b[maxn];
    int lenth=1,Reverse[maxn];
    inline void init(const int x)
    {
    	rg int tim=0;
    	while(lenth<=x)lenth<<=1,tim++;
    	for(rg int i=0;i<lenth;i++)Reverse[i]=(Reverse[i>>1]>>1)|((i&1)<<(tim-1));
    }
    inline void FFT(complex*A,const int fla)
    {
    	for(rg int i=0;i<lenth;i++)if(i<Reverse[i])swap(A[i],A[Reverse[i]]);
    	for(rg int i=1;i<lenth;i<<=1)
    	{
    		const complex w=(complex){cos(PI/i),fla*sin(PI/i)};
    		for(rg int j=0;j<lenth;j+=(i<<1))
    		{
    			complex K=(complex){1,0};
    			for(rg int k=0;k<i;k++,K=K*w)
    			{
    				const complex x=A[j+k],y=A[j+k+i]*K;
    				A[j+k]=x+y;
    				A[j+k+i]=x-y;
    			}
    		}
    	}
    }
    int main()
    {
    	read(n),read(m);
    	init(n+m);
    	for(rg int i=0;i<=n;i++)read(a[i].x);
    	for(rg int i=0;i<=m;i++)read(b[i].x);
    	FFT(a,1),FFT(b,1);
    	for(rg int i=0;i<lenth;i++)a[i]=a[i]*b[i];
    	FFT(a,-1);
    	for(rg int i=0;i<=n+m;i++)print((int)(a[i].x/lenth+0.5)),putchar(' ');
    	return flush(),0;
    }
    

    跑的最慢的点是607ms,常数大大变小了了,大约是 1 4 \frac 14 41的关系
    这个版本就是比较常见的了

    4.2 其它优化

    update by 2019.2.24:
    一般的FFT题中可能会用到多次,由于单位根的计算和乘法都比较慢,这个时候我们可以预处理所有的单位根,开一个数组即可,能够提升FFT的精度(在某写情况下会用到)和速度
    另外我们发现,对于DFT和IDFT的判断会占据比较多的判定时间,如果不想写两个的话其实还有一个技巧,由于IDFT是代入反的单位根,所以将数组非 0 0 0部分翻转后做一遍DFT即可

    5 总结

    这一篇FFT的博客奋战完成了,我也彻底的理解了FFT的具体过程,希望你也能从中获益我可是写的很仔细的,毕竟我以后自己要看
    对于FFT其实代码并不长,能解决的问题却有很多,之后我会写相关的博客来进行总结
    此外不得不提到的是,FFT全程都是double运算,所以就有一些精度上的问题需要注意,为了解决这个问题,有一个叫NTT的算法,是特殊模数的模意义下的多项式乘法,大致和FFT差不多,学有余力可以去学习
    如果你已经把NTT也学完了,可以观测一下我的博客多项式运算学习,学会NTT后,你能在那篇博客学到更多的多项式相关操作
    撒花结束!
    (字数10000+的一篇博客,如果发现有错误欢迎指正!)^ _ ^
    update by 2018.11.28: 对部分的公式的对齐进行了改善
    update by 2018.12.10:lyc的帮助下,对文章进行了部分修正

    展开全文
  • 傅里叶变换 傅里叶变换通常用于在数字信号处理...{2m+1}e^{-\frac{2\pi i}{N}(2m+1)k} X[k]=∑m=0N/2−1​x2m​e−N2πi​(2m)k+∑m=0N/2−1​x2m+1​e−N2πi​(2m+1)k 数论变换NTT 多项式乘法 格密码运算中的FFT实现

    傅里叶变换

    傅里叶变换通常用于在数字信号处理中,对数字信号进行频域和时域的转换,分为离散的和连续的,如下图:
    在这里插入图片描述
    格密码是模整数下的操作,所以只考虑离散情况。

    离散傅里叶变换

    序列{ x n x_n xn} n = 0 N − 1 _{n=0}^{N-1} n=0N1,DFT为:
    在这里插入图片描述
    IDFT为:
    在这里插入图片描述
    IDFT除了指数的符号相反结果需要乘以归一化因子1/N外,与DFT是相同的。
    所以计算DFT的程序代码稍作修改也可以用于IDFT。

    快速傅里叶变换FFT

    O( N 2 N^2 N2)太慢了

    原本的DFT中,将{ x n x_n xn}变换为{ X n X_n Xn}。
    N个点求出N个DFT,每个Xi都与所有的{ x n x_n xn}有关,换言之,要求出 X i X_i Xi,需要知道所有的 x n x_n xn

    可知DFT时间复杂度为O( N 2 N^2 N2)

    对于O( N 2 N^2 N2),通常使用二分法、分治法提高效率。
    FFT正是通过分而治之的思想提高DFT的速度的。

    一分为二

    将一个大问题分解为规模更小但问题相同的两个小问题,通常是分治法的思路。

    DFT计算公式:
    X [ k ] = ∑ n = 0 N − 1 x n e − 2 π i n k / N X[k] = \sum_{n=0}^{N-1}x_ne^{-2\pi ink/N} X[k]=n=0N1xne2πink/N

    N通常是2的幂,那么0~N-1可以按照下标的奇偶平等划分。
    则DFT公式变为:

    X [ k ] = ∑ m = 0 N / 2 − 1 x 2 m e − 2 π i N ( 2 m ) k + ∑ m = 0 N / 2 − 1 x 2 m + 1 e − 2 π i N ( 2 m + 1 ) k X[k] = \sum_{m=0}^{N/2-1}x_{2m}e^{-\frac{2\pi i}{N}(2m)k} + \sum_{m=0}^{N/2-1}x_{2m+1}e^{-\frac{2\pi i}{N}(2m+1)k} X[k]=m=0N/21x2meN2πi(2m)k+m=0N/21x2m+1eN2πi(2m+1)k

    可以看成是{ x n x_n xn}的偶数下标和奇数下标组成的两个序列{ x 1 n x1_n x1n}和{ x 2 n x2_n x2n},两个序列的DFT之和。

    即“DFT of even-indexed part of xn”和“DFT of odd-indexed part of xn”。

    X [ k ] = ∑ m = 0 N / 2 − 1 x 2 m e − 2 π i N ( 2 m ) k + ∑ m = 0 N / 2 − 1 x 2 m + 1 e − 2 π i N ( 2 m + 1 ) k X[k] = \sum_{m=0}^{N/2-1}x_{2m}e^{-\frac{2\pi i}{N}(2m)k} + \sum_{m=0}^{N/2-1}x_{2m+1}e^{-\frac{2\pi i}{N}(2m+1)k} X[k]=m=0N/21x2meN2πi(2m)k+m=0N/21x2m+1eN2πi(2m+1)k

    数论变换NTT

    多项式乘法

    格密码运算中的FFT实现

    展开全文
  • 数据结构实验——用链表实现简单的多项式乘法 提示:文章写完后,目录可以自动生成,如何生成可参考右边的帮助文档 文章目录数据结构实验——用链表实现简单的多项式乘法前言一、实验要求二、实验内容和实验步骤:1...
  • C语言 多项式乘法 算法

    千次阅读 2021-01-13 16:46:21
    多项式乘法 什么是多项式? 由若干个单项式相加组成的代数式叫做多项式(若有减法:减一个数等于加上它的相反数)。 多项式中的每个单项式叫做多项式的项,这些单项式中的最高项次数,就是这个多项式的次数。 ...
  • 我们提出了一种通用的流水线多项式乘法体系结构,以在大约((nlg n)= 4 + n = 2)个时钟周期内计算两个n次多项式的乘积。 此外,我们引入了一些通用的优化技术来减少所需的ROM.storage。 在Spartan-6 FPGA上进行...
  • 多项式乘法(java源代码)
  • 多项式算法1:多项式乘法

    千次阅读 2019-06-02 11:27:26
    多项式简介 对于数域F\mathbb FF,若有∀i∈{1,2,3,⋯&ThinSpace;,n}\forall i\in\{1,2,3,\cdots,n \}∀i∈{1,2,3,⋯,n},则 f(x)=a0+a1x+a2x2+⋯+anx=∑i=1naixif(x)=a_0+a_1x+a_2x^2+\cdots+a_nx=\sum_{i=1}^...
  • 在本文中,我们介绍了几种优化技术,以利用数论变换(NTT)加速多项式乘法。 我们建议预先计算常数多项式的N TT,以减少NTT计算的次数。 为了降低位反转操作的成本,引入了一种优化技术来即时执行它。 我们还提出了...
  • 多项式乘法之FFT快速傅里叶变换

    千次阅读 2018-11-15 00:19:15
    快速傅里叶应用于多项式,来加速多项式乘法使其复杂度从o()变成o(n()),好了先一步步来,这个算法确实很难,但也不是完全理解不了,只是方式不对。 首先,正常的多项式相乘是: A[n]=+x+++~~+,B[n]=~~; C=A[n]*B...
  • 利用卷积求多项式乘法

    千次阅读 2019-12-18 20:17:59
    利用卷积求多项式乘法 对于多项式乘法f(x)=(a0+a1x+……+anxn)∗(b0+b1x+……+bmxm)f(x)=(a_0+a_1x+……+a_nx^n)*(b_0+b_1x+……+b_mx^m)f(x)=(a0​+a1​x+……+an​xn)∗(b0​+b1​x+……+bm​xm),其中 f(x)=c0+c1x...
  • 1.多项式乘法实现

    千次阅读 2018-06-14 17:32:37
    对于稀疏多项式,采用数组存储效率低下,因此考虑采用链表结构,节点包括系数,指数,next指针三个域。多项式的运算中,主要要考虑的是同类项合并的问题,这实际是一个数组元素去重的问题(合并幂相同的节点),因此...
  • import java.util.Random;...//环多项式赋值 public static void fuzhi(int q, Random random, int[] h) { for (int i = 0; i < h.length; i++) { h[i] = random.nextInt(q);//从区间[0,q-1]选择一
  • clc; clear; syms s s1=s^2+s+1; %创建两个多项式 s2=2*s+5; sp1=sym2poly(s1);... %计算两个多项式乘法 ss=collect(ss,s) %展开括号,按降幂排列 结果如下 sp = 2 7 7 5 ss = 2*s^3 + 7*s^2 + 7*s + 5
  • C语言实现GF(2^8)下多项式乘法
  • 文章目录一、解析/思路二、代码三、Main1、...Func1:接收多项式f与拆分的多项式gi,返回fgi的结果 Func2:接收分步相乘的结果T,返回异或结果 二、代码 三、Main f=[0 1 0 0 1 1 0 1]; g=[1 0 0 1 1 0 0 1]; m=[0 0 0

空空如也

空空如也

1 2 3 4 5 ... 20
收藏数 22,190
精华内容 8,876
关键字:

多项式乘法