bzoj 2137: submultiple

Time Limit: 10 Sec Memory Limit: 259 MB

Submit: 239 Solved: 113
[Submit][Status][Discuss]

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

数据范围着实让人头大，，，前九个是一种算法，后11个是一种算法。
先推一下式子：
f(n)= (d|n)∑g(d)^k = πf(pi^ai)=π(g(1)^k+g(pi)^k+...+g(pi^ai)^k)
=π(∑(i=1 to ai+1)i^k)

对数据分治ai小的暴力算出x^k的前缀和，大的高斯消元算出系数再套公式。
至于高斯消元：
先得出 a0+a1*x+a2*x^2+...a(k+1)*x^(k+1)=1^k+2^k+...+x^k,同理可得：
a0+a1*(x+1)+a2*(x+1)^2+...a(k+1)*(x+1)^(k+1)=1^k+...+x^k+(x+1)^k
下式减上式可得：a0*0+a1*(x+1-1)+a2*((x+1)^2-x^2)+....+a(k+1)*((x+1)^(k+1)-x^(k+1))=(x+1)^k
不难发现a0=0(用手指头想想也知道不会带一个常数的，因为0的多少次方都为0)，所以带k+1个等式消元就可以把a1-ak+1求出来啦。
再之后就是套公式环节，注意是求(ai +1)的k次方前缀和。

#include<bits/stdc++.h>
#define ll long long
#define maxn 100005
#define ha 1000000007
using namespace std;
ll a[maxn],n,k,mx=0;
ll ans=1,ci[maxn];
ll b[20][20];inline ll ksm(ll x,ll y){ll an=1;for(;y;y>>=1,x=x*x%ha) if(y&1) an=an*x%ha;return an;
}inline void work(){mx++; ll now=1;for(int i=1;i<=mx;i++){ci[i]=ci[i-1]+ksm(i,k);if(ci[i]>=ha) ci[i]-=ha;}for(int i=1;i<=n;i++){now=now*ci[a[i]+1]%ha;}printf("%lld\n",now);
}inline void solve(){ll len=k+1,now,pre;for(int i=1;i<=len;i++){now=pre=1;for(int j=1;j<=len;j++){now=now*(i+1)%ha;pre=pre*i%ha;if(j==k) b[i][len+1]=now;b[i][j]=(now-pre+ha)%ha;}}for(int i=1;i<=len;i++){if(!b[i][i]){for(int j=i+1;j<=len;j++) if(b[j][i]){for(int l=1;l<=len+1;l++) swap(b[j][l],b[i][l]);break;}}for(int j=i+1;j<=len;j++) if(b[j][i]){ll A=b[i][i],B=b[j][i],C=A/B;while(B){C=A/B;for(int l=i;l<=len+1;l++) b[i][l]=(b[i][l]-C*b[j][l]+ha)%ha;for(int l=i;l<=len+1;l++) swap(b[i][l],b[j][l]);A=b[i][i],B=b[j][i];}}}for(int i=len;i;i--){ll w=b[i][len+1];for(int j=i+1;j<=len;j++) w=(w-b[j][j]*b[i][j]%ha+ha)%ha;b[i][i]=w*ksm(b[i][i],ha-2)%ha;}for(int i=1;i<=n;i++){ll tot=0,now=a[i]%ha+1;for(int j=1;j<=len;j++,now=now*((a[i])%ha+1)%ha) tot=(tot+b[j][j]*now)%ha;ans=ans*tot%ha;}printf("%lld\n",ans);
}int main(){scanf("%lld%lld",&n,&k);for(int i=1;i<=n;i++){scanf("%lld",a+i);mx=max(mx,a[i]);}if(mx<=100000){work();return 0;}solve();return 0;
}