Description

设函数g(N)表示N的约数个数。现在给出一个数M，求出所有M的约数x的g(x)的K次方和。

Input

第一行输入N，K。N表示M由前N小的素数组成。接下来N行，第i+1行有一个正整数Pi，表示第Ai小的素数有 Pi次。等式：

Output

输出一个数，表示答案。只需输出最后答案除以1000000007的余数。

Sample Input

2 3
1
3

Sample Output

900
【样例说明】
M=2^1*3^3=54
M的约数有1,2,3,6,9,18,27,54.约数个数分别为1,2,2,4,3,6,4,8.
Answer=1^3+2^3+2^3+4^3+3^3+6^3+4^3+8^3=900

编号 N K Pi<=
1 50 3 10000
2 50 100 10000
3 50 20101125 10000
4 999 17651851 100000
5 5000 836954247 100000
6 4687 1073741823 100000
7 4321 123456789 100000
8 5216 368756432 100000
9 8080 2^31-1 100000
10 10086 3 2^63-1
11 64970 3 2^63-1
12 71321 3 2^63-1
13 350 5 2^31-1
14 250 6 2^31-1
15 110 7 2^31-1
16 99 8 2^31-1
17 80 9 2^31-1
18 70 10 2^31-1
19 60 11 2^31-1
20 50 12 2^31-1

解题思路：

拿到题感觉一脸不可做，反正不是反演的样子。

先来考虑一下样例：

$2^1*3^3=54$

考虑如何将答案分类。将贡献重新累加起来。

首先$g(d)=\sum\limits_{i=1}^{n}{(1+c_i)}$（）其中$c_i$为第$i$项的次数）

现在考虑如何将$M$的所有因数表示出来：

以54为例

其中每个因数可以按照2的次数分类，可以是0~1次共有2种选法。

对于每种选法，对应的数中3的次数有0~3共4种选法。

而对应3的每一种选法其答案都是$[（2的次数+1）*（3的次数+1）]^k$

那么这种结果可以看做两个多项式乘积的形式,这两个多项式都是(1+2+3+...)的形式，而由于那个指数上的k是可以带入括号的。

那么答案就变成了$f(54)=(1^3+2^3)*(1^3+2^3+3^3+4^3)=9*100=900$

现在答案就可以解出来了，答案就变成了：

$f(M)=\prod\limits_{i=1}^{n}\sum\limits_{j=1}^{p_i}{j^k}$

n是可以枚举的，对于前45分的点，可以暴力快速幂求解。

剩下的怎么办。

后面的那个高次求和$\sum\limits_{i=1}^{n}{i^k}$是非常公式化的，怎么可能没有公式。

高次求和公式会形如：

$\sum\limits_{i=1}^{n}{i^k}=\sum\limits_{i=1}^{k+1}{a_i*n^i}$

所以可以这样求解(orz zwz)：

$a_1x^1+a_2x^2+...+a_{k+1}x^{k+1}+(x+1)^k=a_1(x+1)+a_2(x+1)^2+...+a_{k+1}(x+1)^{k+1}$

二项式定理展开移项：

$\sum\limits_{i=1}^{k+1}{C_i^0x^0}+\sum\limits_{i=2}^{k+1}{C_i^1x^1}+...+\sum\limits_{i=k+1}^{k+1}{C_i^kx^k}=\sum\limits_{i=0}^{k}{C_k^ix^i}$

这个杨辉三角打出矩阵消元就可以解出系数了。

代码：

 #include<cstdio>

 #include<cstring>

 #include<algorithm>

 typedef long long lnt;

 const lnt mod=(lnt)(1e9+);

 lnt n,k;

 lnt maxs;

 lnt C[][];

 lnt F[];

 lnt p[];

 lnt a[];

 lnt b[][];

 lnt c[];

 lnt ksm(lnt x,lnt y)

 {

     lnt ans=;

     while(y)

     {

         if(y&)ans=ans*x%mod;

         x=x*x%mod;

         y=y/;

     }

     return ans;

 }

 lnt squ(lnt x)

 {

     return x%mod*x%mod;

 }

 int main()

 {

     scanf("%lld%lld",&n,&k);

     for(int i=;i<=n;i++)

         scanf("%lld",&p[i]),

         maxs=std::max(p[i],maxs);

     if(maxs<=)

     {

         lnt ans=;

         for(int i=;i<=maxs+;i++)F[i]=(F[i-]+ksm(i,k))%mod;

         for(int i=;i<=n;i++)ans=(ans*F[p[i]+])%mod;

         printf("%lld\n",ans);

         return ;

     }

     if(k==)

     {

         lnt ans=;

         lnt Inv=ksm(,mod-);

         for(int i=;i<=n;i++)

             ans=ans*squ(p[i]%mod+)%mod*squ(p[i]%mod+)%mod*Inv%mod;

         printf("%lld\n",ans);

         return ;

     }

     C[][]=;

     for(int i=;i<=k+;i++)

     {

         C[i][]=;

         for(int j=;j<=i;j++)

             C[i][j]=(C[i-][j-]+C[i-][j])%mod;

     }

     for(int i=;i<=k+;i++)

     {

         for(int j=i;j<=k+;j++)

         {

             b[i][j]=C[j][i-];

         }

     }

     for(int i=;i<=k+;i++)a[i]=C[k][i-];

     for(int i=k+;i>=;i--)

     {

         a[i]=(a[i]*ksm(b[i][i],mod-))%mod;

         for(int j=i-;j;j--)

         {

             (a[j]-=b[j][i]*a[i]%mod)%=mod;

         }

     }

     lnt ans=;

     for(int i=;i<=n;i++)

     {

         lnt tmp=;

         for(int j=;j<=k+;j++)

             tmp=(tmp+(a[j]*ksm(p[i]%mod+,j))%mod)%mod;

         ans=(ans*tmp)%mod;

     }

     printf("%lld\n",(ans%mod+mod)%mod);

     return ;

 }

巴特西

BZOJ2137: submultiple(生成函数，二项式定理)