51nod 1172 Partial Sums V2

3 篇文章 0 订阅
1 篇文章 0 订阅

题目

给出一个数组A,经过一次处理,生成一个数组S,数组S中的每个值相当于数组A的累加,比如:A = {1 3 5 6} => S = {1 4 9 15}。如果对生成的数组S再进行一次累加操作,{1 4 9 15} => {1 5 14 29},现在给出数组A,问进行K次操作后的结果。(输出结果 Mod 10^9 + 7)

分析

发现,每次处理相当于将A卷上一个 I(ai=1) I ( ∀ a i = 1 )
于是机智的我在wiki又发现狄利克雷卷积满足交换律(我居然才知道)
于是快速幂咯,时间复杂度 O(nlog2n) O ( n l o g 2 n ) ,常数巨大,在51nod的老爷机上根本过不了。
然后就TLE得一塌糊涂。
于是找规律,发现 ansk=i+j=kAiCjj+n1 a n s k = ∑ i + j = k A i C j + n − 1 j
然后NTT
但是, 109+7 10 9 + 7 并没有原根,所以祭出三模数NTT(如果有人想用高精度,我也没办法)。
因为 (109+7)2n1023 ( 10 9 + 7 ) 2 ∗ n ≈ 10 23 ,所以找三个 109 10 9 左右的模数。
假设

ansa0(mod m0) a n s ≡ a 0 ( m o d   m 0 )

ansa1(mod m1) a n s ≡ a 1 ( m o d   m 1 )

ansa2(mod m2) a n s ≡ a 2 ( m o d   m 2 )

因为 m0m1m2 m 0 m 1 m 2 会爆long long
所有,通过CRT(中国剩余定理)合并前两项,于是 M=m0m1 M = m 0 ∗ m 1
ansA(mod M) a n s ≡ A ( m o d   M )

ansa2(mod m2) a n s ≡ a 2 ( m o d   m 2 )

ans=kM+A a n s = k M + A
因为
kM+Aansa2(mod m2) k M + A ≡ a n s ≡ a 2 ( m o d   m 2 )

所以
k(a2A)M1(mod m2) k ≡ ( a 2 − A ) M − 1 ( m o d   m 2 )

那么根据上面的式子求出k,通过 kM+A=ans k M + A = a n s 就可以求出ans了。

#include <cmath>
#include <iostream>
#include <cstdio>
#include <cstdlib>
#include <cstring>
#include <algorithm>
#include <queue>
#include <map>
#include <bitset>
#include <set>
const int maxlongint=2147483647;
const long long mo=1e9+7;
const int N=200005;
using namespace std;
long long mod[3]={469762049,998244353,1004535809},M=mod[0]*mod[1],W[3][N];
long long f[3][N],g[3][N],v[3][N],v_[N],jc[3][N],ny[3][N],inv[N];
int n,m,fn;
long long poww(long long x,long long y,int z)
{
    long long s=1;
    x%=mod[z];
    for(;y;x=x*x%mod[z],y>>=1)
        if(y&1) s=s*x%mod[z];
    return s;
}
long long mul(long long x,long long y,long long mo)
{
    long long s=0;
    x%=mo;
    for(;y;x=(x+x)%mo,y>>=1)
        if(y&1) s=(s+x)%mo;
    return s;
}
void NTT(long long *f,int type,int z)
{
    for(int i=0,p=0;i<fn;i++)
    {
        if(i<p) swap(f[i],f[p]);
        for(int j=fn>>1;(p^=j)<j;j>>=1);
    }
    for(int i=2;i<=fn;i<<=1)
    {
        int half=i>>1,pe=fn/i;
        for(int j=0;j<half;j++)
        {
            long long w=!type?W[z][j*pe]:W[z][fn-j*pe];
            for(int k=j;k<fn;k+=i)
            {
                long long x=f[k],y=f[k+half]*w%mod[z];
                f[k]=(x+y)%mod[z],f[k+half]=(x-y+mod[z])%mod[z];
            }
        }
    }
    if(type)
    {
        long long ni=poww(fn,mod[z]-2,z);
        for(int i=0;i<fn;i++) f[i]=f[i]*ni%mod[z];
    }
}
long long CRT(long long a0,long long a1,long long a2)
{
    long long n0=poww(mod[1],mod[0]-2,0),n1=poww(mod[0],mod[1]-2,1);
    long long A=(mul(a0*mod[1]%M,n0%M,M)+mul(a1*mod[0]%M,n1%M,M))%M,n2=poww(M%mod[2],mod[2]-2,2);
    long long k=(a2-A)%mod[2]*n2%mod[2];
    k=(k%mod[2]+mod[2])%mod[2];
    return (k%mo*(M%mo)%mo+A)%mo;
}
int main()
{
    scanf("%d%d",&n,&m);
    for(fn=1;fn<=n*2+2;fn<<=1);
    inv[0]=inv[1]=1;
    for(int i=2;i<=fn;i++) inv[i]=(-(mo/i)*inv[mo%i]%mo+mo)%mo;
    for(int i=0;i<n;i++) scanf("%lld",&f[0][i]),f[1][i]=f[2][i]=f[0][i];
    for(int j=0;j<3;j++)
    {
        W[j][0]=1,W[j][1]=poww(3,(mod[j]-1)/fn,j);
        for(int i=2;i<=fn;i++) W[j][i]=W[j][i-1]*W[j][1]%mod[j];
        for(int i=0;i<n;i++) v[j][i]=1;
        g[j][0]=1;
        for(int i=1;i<n;i++) g[j][i]=g[j][i-1]*(i+m-1)%mo*inv[i]%mo;
    }
    for(int j=0;j<3;j++)
    {
        NTT(f[j],0,j),NTT(g[j],0,j);
        for(int i=0;i<fn;i++) f[j][i]=f[j][i]*g[j][i]%mod[j];
        NTT(f[j],1,j);
    }
    for(int i=0;i<n;i++) printf("%lld\n",CRT(f[0][i],f[1][i],f[2][i]));
}
  • 1
    点赞
  • 0
    收藏
    觉得还不错? 一键收藏
  • 0
    评论

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

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

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

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

抵扣说明:

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

余额充值