zoukankan      html  css  js  c++  java
  • [loj150]挑战多项式

    以NTT为运算基础,即默认支持在$o(nlog n)$的时间内多项式乘法

    二次剩余:称$n$为模$p$意义下的二次剩余,当且仅当存在$x$使得$x^{2}equiv n(mod p)$

    当$p$为素数时,$n$为二次剩余当且仅当$pmid n$或$n^{frac{p-1}{2}}equiv 1(mod p)$,且$[1,p)$中恰有$frac{p-1}{2}$个数为二次剩余

    模意义下开根号:当$p$为素数且$n$为模$p$意义下的二次剩余,求$x$使得$x^{2}equiv n(mod p)$

    当$pmid n$,令$x=0$即可,否则即$n^{frac{p-1}{2}}equiv 1(mod p)$时

    随机$a$使得$a^{2}-n$为非二次剩余,根据前面的结论,单次判定为$o(log p)$,且随机次数为期望$o(1)$

    令$w^{2}equiv a^{2}-n(mod p)$(其中$w$类似于虚数单位$i$,并不存在),即有$w^{p-1}equiv -1(mod p)$

    将之化简,即$nequiv (a-w)(a+w)(mod p)$,将$aequiv a^{p}(mod p)$且$wequiv -w^{p}(mod p)$,代入$a-w$中即可得到$nequiv (a^{p}+w^{p})(a+w)(mod p)$

    进一步的,不难得到$(a+w)^{p}equiv a^{p}+w^{p}(mod p)$,那么$nequiv (a+w)^{p+1}(mod p)$,将其开根号也即有$xequiv (a+w)^{frac{p+1}{2}}(mod p)$

    在快速幂过程中,将每一个数用形如$a+bw$的形式描述(利用$w^{2}equiv a^{2}-n$来消除高次项),且可以证明最终$bequiv 0(mod p)$

    牛顿迭代:给定函数$G$,若$G(F_{0})equiv 0(mod x^{frac{n}{2}})$,求$G(F)equiv 0(mod x^{n})$,做法如下——

    根据泰勒展开,有
    $$
    G(F)=sum_{i=0}^{infty}frac{1}{i!}(F-F_{0})^{i}G^{(i)}(F_{0})equiv 0(mod x^{n})
    $$
    (其中$G^{(i)}$指$G$的$i$阶导数)

    若$G(F)$中若修改$F$的$i$次项系数,不影响$G(F)$中比$i$次数低的项系数,即有$Fequiv F_{0}(mod x^{frac{n}{2}})$

    进而,即$(F-F_{0})^{2}equiv 0(mod x^{n})$,因此上式等价于$G(F_{0})+(F-F_{0})G'(F_{0})equiv 0(mod x^{n})$

    化简后,即$Fequiv F_{0}-frac{G(F_{0})}{G'(F_{0})}(mod x^{n})$,再根据具体的$G$计算即可

    多项式求逆:给定多项式$A$,求多项式$F$使得$AFequiv 1(mod x^{n})$

    将之简单化简,也即解$G(F)=F^{-1}-Aequiv 0(mod x^{n})$

    (这里$G$需要保证$A$为单独一项,否则会导致嵌套)

    $G(F)equiv 0(mod x)$直接对常数项求逆元即可,接下来根据$Fequiv F_{0}-frac{G(F_{0})}{G'(F_{0})}(mod x^{n})$,代入并化简即可得到$Fequiv (2-AF_{0})F_{0}(mod x^{n})$,直接NTT计算即可

    由于每一层计算复杂度为$o(nlog n)$,总复杂度也为$o(nlog n)$

    多项式除法:给定多项式$A$和$B$,求多项式$F$使得$AFequiv B(mod x^{n})$

    先利用多项式求逆找到$AF'equiv 1(mod x^{n})$,令$Fequiv F'B(mod x^{n})$显然成立

    多项式开方:给定多项式$A$,求多项式$F$使得$F^{2}equiv A(mod x^{n})$

    将之简单化简,也即解$G(F)=F^{2}-Aequiv 0(mod x^{n})$

    $G(F)equiv 0(mod x)$直接对常数项开根即可,接下来根据$Fequiv F_{0}-frac{G(F_{0})}{G'(F_{0})}(mod x^{n})$,代入并化简即可得到$Fequiv frac{F_{0}^{2}-A}{2F_{0}}(mod x^{n})$,利用多项式除法计算即可

    由于每一层计算复杂度为$o(nlog n)$,总复杂度也为$o(nlog n)$

    多项式求对数:给定多项式$A$,求多项式$F$使得$e^{F}equiv A(mod x^{n})$

    (关于这个$e^{F}$,将其对$F_{0}=0$这个多项式泰勒展开,即有$e^{F}=sum_{i=0}^{infty}frac{1}{i!}F^{i}$)

    如果使用与之前做法相同的思路,会导致其需要在过程中求多项式的exp,而求exp时又需要求对数,会导致嵌套而无法实现

    事实上,对于多项式求对数,有更为简单(不依赖于exp)的做法——

    对原式两边求导,即$F'e^{F}equiv A(mod x^{n})$

    代入$e^{F}equiv A(mod x^{n})$,即$F'equiv frac{A}{A'}(mod x^{n})$

    再对两边同时积分,即$Fequiv int frac{A}{A'}(mod x^{n})$

    利用多项式求逆以及多项式积分法则,即可做到$o(nlog n)$的复杂度

    多项式求exp:给定多项式$A$,求$Fequiv e^{A}(mod x^{n})$

    将之简单化简,也即解$G(F)=ln F-Aequiv 0(mod x^{n})$

    由于$ln F$的常数项一定为0,因此$A$的常数项也必然为0

    $G(F)equiv 0(mod x)$直接取$F=1$即可($A$常数项为0),接下来根据$Fequiv F_{0}-frac{G(F_{0})}{G'(F_{0})}(mod x^{n})$,代入并化简即可得到$Fequiv (A-ln F_{0}+1)F_{0}(mod x^{n})$,利用多项式求对数计算即可

    由于每一层计算复杂度为$o(nlog n)$,总复杂度也为$o(nlog n)$

    多项式求幂:给定多项式$A$,求$Fequiv A^{k}(mod x^{n})$

    直接快速幂复杂度为$o(nlog^{2}n)$,并不足够优秀

    对原式的两边同时取对数,即$ln Fequiv kln A(mod x^{n})$,再同时exp即$F=e^{kln A}$,利用多项式求对数和exp即可$o(nlog n)$计算

      1 #include<bits/stdc++.h>
      2 using namespace std;
      3 #define mod 998244353
      4 int n,k,tmp;
      5 struct num{
      6     int x,y;
      7     num(int a,int b){
      8         x=a,y=b;
      9     }
     10     num operator + (const num &k)const{
     11         return num((x+k.x)%mod,(y+k.y)%mod);
     12     }
     13     num operator * (const num &k)const{
     14         return num((1LL*x*k.x+1LL*y*k.y%mod*tmp)%mod,(1LL*x*k.y+1LL*k.x*y)%mod);
     15     }
     16 };
     17 struct poly{
     18     vector<int>a;
     19     poly(){
     20         a.clear();
     21     }
     22     poly(int x){
     23         a.clear();
     24         a.push_back(x);
     25     }
     26 }a,ans;
     27 int qpow(int n,int m){
     28     int s=n,ans=1;
     29     while (m){
     30         if (m&1)ans=1LL*ans*s%mod;
     31         s=1LL*s*s%mod;
     32         m>>=1;
     33     }
     34     return ans;
     35 }
     36 num qpow(num n,int m){
     37     num s=n,ans=num(1,0);
     38     while (m){
     39         if (m&1)ans=ans*s;
     40         s=s*s;
     41         m>>=1;
     42     }
     43     return ans;
     44 }
     45 int Sqrt(int k){
     46     int a;
     47     while (1){
     48         int x=1LL*rand()*rand()%mod,p=(1LL*x*x+mod-k)%mod;
     49         if ((p)&&(qpow(p,(mod-1)/2)!=1)){
     50             a=x;
     51             break;
     52         }
     53     }
     54     tmp=(1LL*a*a+mod-k)%mod;
     55     int ans=qpow(num(a,1),(mod+1)/2).x;
     56     return min(ans,mod-ans);
     57 }
     58 poly Add(poly a,int l){
     59     while (a.a.size()<(1<<l))a.a.push_back(0);
     60     return a;
     61 }
     62 poly Int(poly a,int n){
     63     for(int i=(1<<n)-1;i;i--)a.a[i]=1LL*a.a[i-1]*qpow(i,mod-2)%mod;
     64     a.a[0]=0;
     65     return a;
     66 }
     67 poly derive(poly a,int n){
     68     for(int i=1;i<(1<<n);i++)a.a[i-1]=1LL*a.a[i]*i%mod;
     69     a.a[(1<<n)-1]=0;
     70     return a;
     71 }
     72 void ntt(poly &a,int n,int p){
     73     for(int i=0;i<(1<<n);i++){
     74         int s=0;
     75         for(int j=0;j<n;j++)
     76             if (i&(1<<j))s+=(1<<n-j-1);
     77         if (i<s)swap(a.a[i],a.a[s]);
     78     }
     79     for(int i=2;i<=(1<<n);i*=2){
     80         int s=qpow(3,(mod-1)/i);
     81         if (p)s=qpow(s,mod-2);
     82         for(int j=0;j<(1<<n);j+=i){
     83             for(int k=0,ss=1;k<(i>>1);k++,ss=1LL*ss*s%mod){
     84                 int x=a.a[j+k],y=1LL*a.a[j+k+(i>>1)]*ss%mod;
     85                 a.a[j+k]=(x+y)%mod;
     86                 a.a[j+k+(i>>1)]=(x+mod-y)%mod;
     87             }
     88         }
     89     }
     90     if (p){
     91         int s=qpow((1<<n),mod-2);
     92         for(int i=0;i<(1<<n);i++)a.a[i]=1LL*a.a[i]*s%mod;
     93     }
     94 }
     95 poly add(poly a,poly b,int n){
     96     for(int i=0;i<(1<<n);i++)a.a[i]=(a.a[i]+b.a[i])%mod;
     97     return a;
     98 }
     99 poly mul(poly a,int k,int n){
    100     for(int i=0;i<(1<<n);i++)a.a[i]=1LL*a.a[i]*k%mod;
    101     return a;
    102 }
    103 poly mul(poly a,poly b,int n){
    104     a=Add(a,n+1),b=Add(b,n+1);
    105     for(int i=(1<<n);i<(1<<n+1);i++)a.a[i]=b.a[i]=0;
    106     ntt(a,n+1,0);
    107     ntt(b,n+1,0);
    108     poly ans;
    109     for(int i=0;i<(1<<n+1);i++)ans.a.push_back(1LL*a.a[i]*b.a[i]%mod);
    110     ntt(ans,n+1,1);
    111     for(int i=0;i<(1<<n);i++)ans.a.pop_back();
    112     return ans;
    113 }
    114 poly inv(poly a,int n){
    115     if (!n)return poly(qpow(a.a[0],mod-2));
    116     poly s=Add(inv(a,n-1),n),ans=mul(a,s,n);
    117     ans=add(mul(ans,mod-1,n),Add(poly(2),n),n);
    118     return mul(ans,s,n);
    119 }
    120 poly sqrt(poly a,int n){
    121     if (!n)return poly(Sqrt(a.a[0]));
    122     poly s=Add(sqrt(a,n-1),n),ans=mul(s,s,n);
    123     ans=mul(add(ans,a,n),(mod+1)/2,n);
    124     return mul(ans,inv(s,n),n);
    125 }
    126 poly ln(poly a,int n){
    127     return Int(mul(derive(a,n),inv(a,n),n),n);
    128 }
    129 poly exp(poly a,int n){
    130     if (!n)return poly(1);
    131     poly s=Add(exp(a,n-1),n),ans=ln(s,n);
    132     ans=add(add(a,mul(ans,mod-1,n),n),Add(poly(1),n),n);
    133     return mul(ans,s,n);
    134 }
    135 poly qpow(poly a,int n,int m){
    136     return exp(mul(ln(a,n),m,n),n);
    137 }
    138 int main(){
    139     srand(time(0));
    140     scanf("%d%d",&n,&k);
    141     a=Add(a,17);
    142     for(int i=0;i<=n;i++)scanf("%d",&a.a[i]);
    143     ans=exp(Int(inv(sqrt(a,17),17),17),17);
    144     ans=add(mul(ans,mod-1,17),add(a,Add(poly(mod+2-a.a[0]),17),17),17);
    145     ans=derive(qpow(add(ln(ans,17),Add(poly(1),17),17),17,k),17);
    146     for(int i=0;i<n;i++)printf("%d ",ans.a[i]);
    147 }
    View Code
  • 相关阅读:
    React入门实例
    【C语言】一些重要的知识点
    【C语言】字符串模块
    【C语言】指针模块
    贝尔曼福特算法
    dijkstra算法
    拓扑序列
    树和图的广度优先遍历
    树和图的深度优先遍历
    回溯剪枝,dfs,bfs
  • 原文地址:https://www.cnblogs.com/PYWBKTDA/p/14736495.html
Copyright © 2011-2022 走看看