BZOJ 3426 && CodeChef/CHANGE

Description：

$~~~~~~$ 经典问题，给你 $N(N<=50)$ 种面值的硬币，每种可以使用无限次，第 $i$ 硬币的面值为 $D_i(D_i<=500)$ ，所有的 $D_i$ 互质，现在你要支付 $C(C<=10^{100})$ 元钱，问有多少种支付方法。答案对 $10^9+7$ 取模。多组数据，数据组数 $T<=5$ 。

$~~~~~~$
$~~~~~~$
$~~~~~~$

Solution：

$~~~~~~$ 首先对于面值为 $D_i$ 的硬币构造多项式:

$$F_i=\sum_{p=0}^{+\infty}x^{D_i*p}$$
$~~~~~~$ 令 $$G(x)=\prod_{i=1}^N F_i$$
$~~~~~~$ 那么显然把 $G(x)$ 中 $x^C$ 项的系数就是答案。
$~~~~~~$ 但是此题 $C$ 的范围异常的大，这让我们不能直接计算乘法之后的多项式，得考虑别的算法。
$~~~~~~$ 根据幂级数求和，易得 $$F_i=\frac{1}{1-x^{D_i}}$$
$~~~~~~$ 所以 $$G(x)=\prod_{i=1}^N \frac{1}{1-x^{D_i}}=\frac{1}{(1-x)^N}\prod_{i=1}^N \frac{1}{1+x+...+x^{D_i-1}}$$
$~~~~~~$ 把 $G(x)$ 分式分解得： $$G(x)=\frac{A(x)}{(1-x)^N}+\sum_{i=1}^N \frac{B_i(x)}{1+x+...+x^{D_i-1}}=\prod_{i=1}^N \frac{1}{1-x^{D_i}}$$
$~~~~~~$ 乱证明一下发现 $A(x)$ 是一个次数小于 $N$ 的多项式， $B_i(x)$ 是一个次数小于 $D_i-1$ 次的多项式。
$~~~~~~$ 考虑多项式 $B_i(x)$ ，那么有
$~~~~~~$ $$\frac{B_i(x)}{1+x+...+x^{D_i-1}}=\prod_{j=1}^N \frac{1}{1-x^{D_j}}-\frac{A(x)}{(1-x)^N}-\sum_{j=1,j!=i}^N \frac{B_j(x)}{1+x+...+x^{D_j-1}}$$
$~~~~~~$ 令 $$R_i(x)=\frac{A(x)}{(1-x)^N}+\sum_{j=1,j!=i}^N \frac{B_j(x)}{1+x+...+x^{D_j-1}}$$
$~~~~~~$ 那么 $$\frac{B_i(x)}{1+x+...+x^{D_i-1}}=\prod_{j=1}^N \frac{1}{1-x^{D_j}}-R_i(x)$$
$~~~~~~$ 将上式两边乘以左边分母有：
$~~~~~~$ $$B_i(x)=\prod_{j=1}^N \frac{1}{1-x^{D_j}}\times(1+x+...+x^{D_j-1})-{R_i(x)\times(1+x+...+x^{D_i-1})}$$
$~~~~~~$ $$=\frac{1}{(1-x)\prod_{j=1,j!=i}^N(1-x^{D_j})}-{R_i(x)\times(1+x+...+x^{D_i-1})}$$
$~~~~~~$ 那么对于 $x=\omega_{D_i}^t~(1<=t<D_i)$ ，根据等比数列求和， $1+x+...+x^{D_i-1}=0$ ，因为所有的 $D_i$ 互质，当 $x=\omega_{D_i}^t$ 时， $R_i(x)\times(1+x+...+x^{D_i-1})=0$
$~~~~~~$ 所以 $$B_i(\omega_{D_i}^t)=\frac{1}{(1-\omega_{D_i}^t)\prod_{j=1,j!=i}^N(1-\omega_{D_i}^{t\times D_j})}$$
$~~~~~~$ 考虑使用等比数列求和可以十分容易的证明下式 $(j!=0)$ ： $$-\frac{1}{D_i}\times\sum_{p=0}^{D_i-1}(p+1)*\omega_{D_i}^{~p\times j}=\frac{1}{{1-\omega_{D_i}^{~j}}}$$
$~~~~~~$ 所以 $$B_i(\omega_{D_i}^t)=(-\frac{1}{D_i})^N\times(\sum_{p=0}^{D_i-1}(p+1)\times\omega_{D_i}^{~p\times t})\times \prod_{j=1,j!=i}^N(\sum_{p=0}^{D_i-1}(p+1)\times\omega_{D_i}^{~p\times t \times D_j})$$
$~~~~~~$ 令 $$B_i(\omega_{D_i}^t)=(-\frac{1}{D_i})^N\times\sum_{j=0}^{D_i-1}b_j\times\omega_{D_i}^j$$
$~~~~~~$ 如果要求出 $b_i$ ，那么需要做 $N$ 次形如下式的卷积运算(以下所有的指数和下标运算都是在模 $D_i$ 意义下进行的) $$(\sum_{p=0}^{D_i-1}a_p\omega_{D_i}^p) \times (\sum_{p=0}^{D_i-1}(p+1)\omega_{D_i}^{p\times q})=\sum_{p=0}^{D_i-1}c_p\omega_{D_i}^p$$
$~~~~~~$ $$c_p=\sum_{k=0}^{D_i-1}(k+1)\times a_{p-q\times k}=\sum_{k=0}^{D_i-1}k\times a_{p-q\times k}+\sum_{k=0}^{D_i-1}a_k=\sum_{k=0}^{D_i-2}(k+1)\times a_{p-q\times (k+1)}+\sum_{k=0}^{D_i-1}a_k$$
$~~~~~~$ $$=\sum_{k=0}^{D_i-2}(k+1)\times a_{(p-q)-q\times k}+\sum_{k=0}^{D_i-1}a_k=\sum_{k=0}^{D_i-1}(k+1)\times a_{(p-q)-q\times k}+\sum_{k=0}^{D_i-1}a_k-D_i\times a_{(p-q)-q\times(D_i-1)}$$
$~~~~~~$ $$=c_{p-q}+\sum_{k=0}^{D_i-1}a_k-D_i\times a_{(p-q)-q\times(D_i-1)}$$
$~~~~~~$ 所以对于一次卷积运算，可以先求出 $c_0$ ,通过不断的将下标加 $q$ 求出下一个 $c$ ，所以这样就能够在 $O(D_i)$ 的时间里计算一次卷积。所以可以在 $O(N\times D_i)$ 的时间里求出 $b_0,b_1,...,b_{D_i-1}$
$~~~~~~$ 观察发现，对于不同的 $t$ $$B_i(\omega_{D_i}^t)=\sum_{p=0}^{D_i-2}(r_p-r_{D_i-1})\times(\omega_{D_i}^t)^p$$
$~~~~~~$ $r_0,r_1,...,r_{D_i-1}$ 都是相同的，对应等于 $t=1$ 时的 $b_0,b_1,...,b_{D_i-1}$ 。
$~~~~~~$ 由于 $B_i(x)$ 是一个 $D_i-2$ 次多项式，且对于 $B_i(x)$ 当 $x=\omega_{D_i}^t (t=1,2..,D_i-1)$ 时都满足 $$B_i(\omega_{D_i}^t)=\sum_{p=0}^{D_i-2}(r_p-r_{D_i-1})\times(\omega_{D_i}^t)^p$$
$~~~~~~$ 那么断定 $$B_i(x)=\sum_{p=0}^{D_i-2}(r_p-r_{D_i-1})\times x^p$$
$~~~~~~$ 考虑 $$\frac{B_i(x)}{1+x+...+x^{D_i-1}}$$ 对答案的贡献。
$~~~~~~$ $$\frac{B_i(x)}{1+x+...+x^{D_i-1}}=\frac{B_i(x)\times (1-x)}{1-x^{D_i}}=B_i(x)\times (1-x) \times \sum_{p=0}^{+\infty}x^{D_i\times p}$$
$~~~~~~$ 由于 $B_i(x)$ 是一个 $D_i-2$ 次多项式 ，显然 $\frac{B_i(x)}{1+x+...+x^{D_i-1}}$ 中 $x^{D_i\times u+v}(0 < = u,0<=v < D_i)$ 的系数为 $(r_v-r_{D_i-1})-(r_{v-1}-r_{D_i-1})=r_v-r_{v-1}$
$~~~~~~$ 所以把 $C$ 对于 $D_i$ 取模后就可以算出 $\frac{B_i(x)}{1+x+...+x^{D_i-1}}$ 对答案的贡献了。
$~~~~~~$ 所以求出所有 $D_i(x)$ 对答案的贡献的时间复杂度是 $O(N^2\times D)$
$~~~~~~$ 之后 $\frac{A(x)}{(1-x)^N}$ 对答案的贡献，乱证明一下发现 $A(x)$ 是一个多项式，并且每一项的系数都是常数。
$~~~~~~$ 并且 $$\frac{1}{(1-x)^N}=\sum_{i=0}^{+\infty}h(i)x^i$$
$~~~~~~$ 其中 $h(x)$ 是一个关于 $x$ 的 $N-1$ 次多项式
$~~~~~~$ 又因为 $A(x)$ 是一个系数为常数的多项式，所以可以断定
$~~~~~~$ $$\frac{A(x)}{(1-x)^N}=\sum_{i=0}^{+\infty}S(i)x^i$$
$~~~~~~$ 其中 $S(x)$ 是一个关于 $x$ 的 $N-1$ 次多项式
$~~~~~~$ 由于我们能够求出 $B_i(x)$ 对于答案的贡献，并且对于 $C$ 比较小的时候的答案我们是可以 $dp$ 求出的，那么就可以求出 $C$ 比较小的时 $\frac{A(x)}{(1-x)^N}$ 对答案的贡献。
$~~~~~~$ 所以先用 $dp$ 求出 $C=0,1,2,...,N-1$ 时的答案，然后求出 $B_i(x)$ 的贡献，那么就可以求出 $S(0),S(1),...,S(N-1)$ ，然后插值或者高斯消元可以解出多项式 $S(x)$ 各项的系数。之后就得出了 $S(x)$ 的表达式 。之后计算对于题目给出的 $C$ 的答案， $S(C)$ 就是 $\frac{A(x)}{(1-x)^N}$ 对答案的贡献。我用的是高斯消元，这部分的时间复杂度是 $O(N^3)$
$~~~~~~$ 至此此题终结。
$~~~~~~$
$~~~~~~$
$~~~~~~$

Code：

#include <cstdio>
#include <cstdlib>
#include <cstring>
#include <cmath>
#include <algorithm>
#include <iostream>

using namespace std;

const long long Mod=1e9+7;

int N;
int T;
int D[60]={0};
long long B[60][510]={{0}};
long long help[510]={0};
long long A[110]={0};
long long F[110]={0};
long long fc[110][110]={{0}};

struct Bignum
{
    int cnmbdctr[200];
}C={0};

void read_(struct Bignum &R)
{
    int len=0;
    char ch=getchar();
    for(;ch<'0' || ch>'9';ch=getchar());
    for(;ch>='0' && ch<='9';ch=getchar())
        R.cnmbdctr[++len]=ch-'0';
    for(int i=(len>>1);i>=1;i--)
        swap(R.cnmbdctr[i],R.cnmbdctr[len-i+1]);
    R.cnmbdctr[0]=len;
    return;
}

int Module(Bignum R,long long MM)
{
    long long remains=0;
    for(int i=R.cnmbdctr[0];i>=1;i--)
        remains=(remains*10+R.cnmbdctr[i])%MM;
    return (int)remains;
}

long long power(long long a,long long k)
{
    long long o=1;
    for(;k>0;k>>=1)
    {
        if(k&1) o=o*a%Mod;
        a=a*a%Mod;
    }
    return o;
}

void Guass()
{
    for(int i=0;i<=N;i++)
    {
        int gg=i;
        for(;gg<=N;gg++)
            if(fc[gg][i]!=0) break;
        if(gg==N+1) continue;
        swap(fc[i],fc[gg]);
        long long Inv=power(fc[gg][i],Mod-2);
        for(int j=i;j<=N+1;j++)
            fc[gg][j]=fc[gg][j]*Inv%Mod;
        for(int j=0;j<=N;j++)
        {
            if(j==i) continue;
            if(fc[j][i]!=0)
            {
                long long mult=fc[j][i];
                for(int p=i;p<=N+1;p++)
                    fc[j][p]=(fc[j][p]-fc[gg][p]*mult%Mod+Mod)%Mod;
            }
        }
    }
    return;
}

long long getansb(int i,int PP)
{
    if(PP==0)
        return (B[i][0]-B[i][D[i]-1]+Mod)%Mod;
    else return (B[i][PP]-B[i][PP-1]+Mod)%Mod;
}

int main()
{
    cin>>T;
    for(;T>0;T--)
    {
        cin>>N;
        memset(C.cnmbdctr,0,sizeof(C.cnmbdctr));
        read_(C);
        for(int i=1;i<=N;i++)
            scanf("%d",&D[i]);
        memset(B,0,sizeof(B));
        memset(help,0,sizeof(help));
        for(int i=1;i<=N;i++)
        {
            long long Inv=power(Mod-D[i],Mod-2);
            Inv=power(Inv,N);
            for(int j=0;j<D[i];j++)
                B[i][j]=j+1;
            for(int p=1;p<=N;p++)
            {
                if(p==i) continue;
                long long Sum=0;
                for(int j=0;j<D[i];j++)
                {
                    help[j]=B[i][j];
                    B[i][j]=0;
                    Sum=(Sum+help[j])%Mod;
                }
                for(int j=0;j<D[i];j++)
                    B[i][0]=(B[i][0]+(j+1)*help[(D[i]-j*D[p]%D[i])%D[i]])%Mod;
                int Cnt=0;
                for(int j=1;j<D[i];j++)
                {
                    int Next=(Cnt+D[p])%D[i];
                    B[i][Next]=(B[i][Cnt]+Sum-(D[i]*help[(Cnt-D[p]*(D[i]-1)%D[i]+D[i])%D[i]])%Mod+Mod)%Mod;
                    Cnt=Next;
                }
            }
            for(int j=0;j<D[i];j++)
                B[i][j]=B[i][j]*Inv%Mod;
        }
        memset(F,0,sizeof(F));
        F[0]=1;
        for(int i=1;i<=N;i++)
        {
            for(int j=D[i];j<=N+20;j++)
                F[j]=(F[j]+F[j-D[i]])%Mod;
        }
        memset(fc,0,sizeof(fc));

        for(int i=0;i<=N;i++)
        {
            fc[i][0]=1;
            long long Cnt=i;
            for(int j=1;j<=N;j++)
            {
                fc[i][j]=Cnt;
                Cnt=Cnt*i%Mod;
            }
            fc[i][N+1]=F[i];
            for(int j=1;j<=N;j++)
                fc[i][N+1]=(fc[i][N+1]-getansb(j,i%D[j])+Mod)%Mod;
        }
        Guass();
        long long Ans=0;
        for(int i=1;i<=N;i++)
        {
            int PP=Module(C,(long long)D[i]);
            Ans=(Ans+getansb(i,PP))%Mod;
        }
        long long Cnt=1;
        for(int j=0;j<=N;j++)
        {
            Ans=(Ans+fc[j][N+1]*Cnt)%Mod;
            Cnt=Cnt*Module(C,Mod)%Mod;
        }
        cout<<Ans<<endl;
    }
    return 0;
}