【BZOJ2137】submultiple
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
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
【样例说明】
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
题解：数据明显分为两部分，一部分pi很小，一部分K很小，需要分别处理。
不难发现，对于$n=prodlimits_ip_i^{c_i}$，$ans=prodlimits_i(1^k+2^k+...{(c_i+1)}^k)$。这就是伯努利数的形式。
当pi很小时，我们可以预处理出$i^k$的前缀和，然后暴力计算。当k很小时，我们知道伯努利数可以表示成一个k+1次的多项式形式，可以暴力算出前k+1个值得到k+1个方程，然后进行模意义下的高斯消元求出多项式的系数，最后将p带入多项式即可。
#include <cstdio>
#include <cstring>
#include <iostream>
using namespace std;
typedef long long ll;
const ll P=1000000007;
int n,m;
ll ans;
ll v[80000],s[100010];
ll A[20][20];
inline ll rd()
{
ll ret=0,f=1; char gc=getchar();
while(gc<'0'||gc>'9') {if(gc=='-') f=-f; gc=getchar();}
while(gc>='0'&&gc<='9') ret=ret*10+(gc^'0'),gc=getchar();
return ret*f;
}
inline ll pm(ll x,ll y)
{
x%=P;
ll z=1;
while(y)
{
if(y&1) z=z*x%P;
x=x*x%P,y>>=1;
}
return z;
}
int main()
{
n=rd(),m=rd(),ans=1;
int i,j,k;
for(i=1;i<=n;i++) v[i]=rd();
if(m<=12)
{
for(i=1;i<=m+1;i++)
{
for(j=1;j<=m+1;j++) A[i][j]=pm(i,j);
for(j=1;j<=i;j++) A[i][m+2]=(A[i][m+2]+pm(j,m))%P;
}
for(i=1;i<=m+1;i++)
{
for(j=i;j<=m+1;j++) if(A[j][i]) break;
if(i!=j) for(k=i;k<=m+2;k++) swap(A[i][k],A[j][k]);
ll tmp=pm(A[i][i],P-2);
for(k=i;k<=m+2;k++) A[i][k]=A[i][k]*tmp%P;
for(j=1;j<=m+1;j++) if(i!=j)
{
tmp=A[j][i];
for(k=i;k<=m+2;k++) A[j][k]=(A[j][k]-A[i][k]*tmp%P+P)%P;
}
}
for(i=1;i<=n;i++)
{
ll tmp=0;
for(j=1;j<=m+1;j++) tmp=(tmp+A[j][m+2]*pm(v[i]+1,j))%P;
ans=ans*tmp%P;
}
printf("%lld",ans);
return 0;
}
for(i=1;i<=100000;i++) s[i]=(s[i-1]+pm(i,m))%P;
for(i=1;i<=n;i++) ans=ans*s[v[i]+1]%P;
printf("%lld",ans);
return 0;
}