51nod 1348 乘积之和（NTT+分治）

题目描述

传送门

题目大意：给出由N个正整数组成的数组A，有Q次查询，每个查询包含一个整数K，从数组A中任选K个（K <= N）把他们乘在一起得到一个乘积。求所有不同的方案得到的乘积之和，由于结果巨大，输出Mod 100003的结果即可

题解

一次卷积的上限大概是1016左右，选择两个109左右的NTT模数，每次做到时候两个分别做，再用中国剩余定理合并即可。
分治，f[i]表示选择i个数的乘积和，用NTT将区间(l,mid)，(mid+1,r)合并起来得到(l,r)的答案。

代码

#include<iostream>
#include<cstdio>
#include<algorithm>
#include<cstring>
#include<cmath>
#define N 200003
#define LL long long
#define mod 100003
using namespace std;
LL f[20][N],M[3],ans[N];
int n,R[N],m,val[N];
LL quickpow(LL num,LL x,LL p)
{
    LL base=num%p; LL ans=1;
    while (x) {
        if (x&1) ans=ans*base%p;
        x>>=1;
        base=base*base%p;
    }
    return ans;
}
void NTT(LL a[N],int n,int opt,LL p)
{
    for (int i=0;i<n;i++) 
     if (i>R[i]) swap(a[i],a[R[i]]);
    for (int i=1;i<n;i<<=1) {
        LL wn=quickpow(3,(p-1)/(i<<1),p);
        for (int p1=i<<1,j=0;j<n;j+=p1) {
            LL w=1;
            for (int k=0;k<i;k++,w=w*wn%p) {
                LL x=a[j+k],y=w*a[j+k+i]%p;
                a[j+k]=(x+y)%p; a[j+k+i]=(x-y+p)%p;
            }
        } 
    }
    if (opt==-1) reverse(a+1,a+n);
}
LL mul(LL num,LL x,LL p)
{
    LL ans=0; LL base=num;
    while (x) {
        if (x&1) ans=(ans+base)%p;
        x>>=1;
        base=(base+base)%p;
    }
    return ans;
}
LL china(LL a1,LL a2)
{
    LL MM=M[1]*M[2];
    LL x=quickpow(M[2],M[1]-2,M[1]),x1=quickpow(M[1],M[2]-2,M[2]);
    LL A=(mul(a1*M[2]%MM,x%MM,MM)+mul(a2*M[1]%MM,x1%MM,MM))%MM;
    return A%mod;
}
void solve(int l,int r,int x)
{
    if (l==r) {
        f[x][0]=1; f[x][1]=val[l]%mod;
        return;
    }
    int mid=(l+r)/2;
    int len=r-l+1; int L=0,n1=1;
    for (n1=1;n1<=len;n1<<=1) L++;
    solve(l,mid,x+1); 
    LL a[n1+10],b[n1+10],a1[n1+10],b1[n1+10];
    for (int i=0;i<=mid-l+1;i++) a[i]=a1[i]=f[x+1][i];
    for (int i=mid-l+2;i<=n1;i++) a[i]=0,a1[i]=0;
    solve(mid+1,r,x+1);
    for (int i=0;i<=r-mid;i++) b[i]=b1[i]=f[x+1][i];
    for (int i=r-mid+1;i<=n1;i++) b[i]=0,b1[i]=0;;
    for (int i=0;i<=n1;i++) R[i]=(R[i>>1]>>1)|((i&1)<<(L-1));
    NTT(a,n1,1,M[1]); NTT(b,n1,1,M[1]);
    for (int i=0;i<n1;i++) a[i]=a[i]*b[i]%M[1];
    NTT(a,n1,-1,M[1]);
    LL inv=quickpow(n1,M[1]-2,M[1]);
    for (int i=0;i<=len;i++) a[i]=a[i]*inv%M[1];
    NTT(a1,n1,1,M[2]); NTT(b1,n1,1,M[2]);
    for (int i=0;i<=n1;i++) a1[i]=a1[i]*b1[i]%M[2];
    NTT(a1,n1,-1,M[2]);
    inv=quickpow(n1,M[2]-2,M[2]);
    for (int i=0;i<=len;i++) a1[i]=a1[i]*inv%M[2];
    for (int i=0;i<=len;i++) f[x][i]=china(a[i],a1[i]);
}
int main()
{
    scanf("%d%d",&n,&m);
    M[1]=998244353,M[2]=1004535809;
    for (int i=1;i<=n;i++) scanf("%d",&val[i]);
    solve(1,n,1);
    for (int i=1;i<=m;i++) {
        int x; scanf("%d",&x);
        printf("%lld\n",f[1][x]);
    }
}