bzoj5300: [Cqoi2018]九连环(fft优化高精+快速幂)

17 篇文章 0 订阅
10 篇文章 0 订阅

传送门
题意不好描述(自己看样例解释)


首先可以推出一个递推式: f n = f n − 1 + 2 f n − 2 + 1 f_n=f_{n-1}+2f_{n-2}+1 fn=fn1+2fn2+1
然后可以构造两个等式:

  1. ( f n + f n − 1 + 1 ) = 2 ( f n − 1 + f n − 2 + 1 ) (f_n+f_{n-1}+1)=2(f_{n-1}+f_{n-2}+1) (fn+fn1+1)=2(fn1+fn2+1)
  2. ( f n − 2 f n − 1 − 1 2 ) = − ( f n − 1 − f n − 1 2 ) (f_n-2f_{n-1}-\frac12)=-(f_{n-1}-f_n-\frac12) (fn2fn121)=(fn1fn21)

然后变形一波可以得到 f n = ⌊ 2 n + 1 − 1 3 ⌋ f_n=\left\lfloor\frac{2^{n+1}-1}3\right\rfloor fn=32n+11
接着我们惊奇的发现没有模数。
只能 f f t fft fft加快速幂了。
代码:

#include<bits/stdc++.h>
#define ri register int
using namespace std;
typedef long double ldb;
struct Cp{
    ldb x,y;
    friend inline Cp operator+(const Cp&a,const Cp&b){return (Cp){a.x+b.x,a.y+b.y};}
    friend inline Cp operator-(const Cp&a,const Cp&b){return (Cp){a.x-b.x,a.y-b.y};}
    friend inline Cp operator*(const Cp&a,const Cp&b){return (Cp){a.x*b.x-a.y*b.y,a.x*b.y+a.y*b.x};}
    friend inline Cp operator*(const Cp&a,const ldb&b){return (Cp){a.x*b,a.y*b};}
    friend inline Cp operator/(const Cp&a,const ldb&b){return (Cp){a.x/b,a.y/b};}
};
const ldb pi=acos(-1.0);
int lim,tim;
vector<Cp>A,B;
vector<int>pos;
inline void init(const int&up){ 
    lim=1,tim=0;
    while(lim<=up)lim<<=1,++tim;
    pos.resize(lim),A.resize(lim),B.resize(lim),pos[0]=0;
    for(ri i=0;i<lim;++i)pos[i]=(pos[i>>1]>>1)|((i&1)<<(tim-1));
}
inline void fft(vector<Cp>&a,const int&type){
    for(ri i=0;i<lim;++i)if(i<pos[i])swap(a[i],a[pos[i]]);
    for(ri mid=1;mid<lim;mid<<=1){
        Cp wn=(Cp){cos(pi/mid),type*sin(pi/mid)};
        for(ri j=0,len=mid<<1;j<lim;j+=len){
            Cp w=(Cp){1,0},a0,a1;
            for(ri k=0;k<mid;++k,w=w*wn){
                a0=a[j+k],a1=a[j+k+mid]*w;
                a[j+k]=a0+a1,a[j+k+mid]=a0-a1;
            }
        }
    }
    if(type==-1)for(ri i=0;i<lim;++i)a[i]=a[i]/(ldb)lim;
}
struct poly{
    vector<Cp>a;
    poly(int k=0,Cp x=(Cp){0,0}){a.resize(k+1),a[k]=x;}
    inline int deg()const{return a.size()-1;}
    inline poly extend(const int&k){poly ret=*this;return ret.a.resize(k+1),ret;}
    inline Cp&operator[](const int&k){return a[k];}
    inline const Cp&operator[](const int&k)const{return a[k];}
    friend inline poly operator*(const poly&a,const poly&b){
        int n=a.deg(),m=b.deg();
        init(n+m);
        poly ret(lim);
        for(ri i=0;i<=n;++i)A[i]=a[i];
        for(ri i=0;i<=m;++i)B[i]=b[i];
        for(ri i=n+1;i<lim;++i)A[i]=(Cp){0,0};
        for(ri i=m+1;i<lim;++i)B[i]=(Cp){0,0};
        fft(A,1),fft(B,1);
        for(ri i=0;i<lim;++i)A[i]=A[i]*B[i];
        return fft(A,-1),ret.a=A,ret;
    }
    friend inline poly operator*(const poly&a,const ldb&b){
        poly ret(a.deg());
        for(ri i=0;i<=a.deg();++i)ret[i]=a[i]*b;
        return ret;
    }
    friend inline poly operator/(const poly&a,const ldb&b){
        poly ret(a.deg());
        for(ri i=0;i<=a.deg();++i)ret[i]=a[i]/b;
        return ret;
    }
};
inline poly update(const poly&a){
    vector<int>res;
    int sum=0;
    for(ri val,i=0,up=a.deg();i<=up;++i)val=((int)(a[i].x+0.5))+sum,sum=val/10,res.push_back(val%10);
    while(sum)res.push_back(sum%10),sum/=10;
    while(!res[res.size()-1]&&res.size()!=1)res.pop_back();
    poly ret(res.size()-1);
    for(ri i=0;i<=ret.deg();++i)ret[i]=(Cp){res[i],0};
    return ret;
}
inline void ksm(int p){
    poly ret,a;
    ret[0]=(Cp){1,0},a[0]=(Cp){2,0};
    while(p){
        a=update(a);
        if(p&1)ret=update(ret),ret=ret*a;
        a=a*a,p>>=1;
    }
    ret=update(ret),reverse(ret.a.begin(),ret.a.end());
    vector<int>ans;
    for(ri val,divv=0,i=0;i<=ret.deg();++i)val=(int)ret[i].x+divv,ans.push_back(val/3),divv=val%3*10;
    int pos=0;
    while(!ans[pos]&&pos!=ans.size()-1)++pos;
    for(ri i=pos;i<ans.size();++i)cout<<ans[i];
    puts("");
}
char s[10005];
int main(){
    int n,m;
    scanf("%d",&m);
    while(m--)scanf("%d",&n),ksm(n+1);
    return 0;
}
  • 1
    点赞
  • 1
    收藏
    觉得还不错? 一键收藏
  • 0
    评论

“相关推荐”对你有帮助么?

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

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

抵扣说明:

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

余额充值