【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

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

题解：数据明显分为两部分，一部分pi很小，一部分K很小，需要分别处理。

不难发现，对于$n=\prod\limits_ip_i^{c_i}$，$ans=\prod\limits_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;

}

巴特西

【BZOJ2137】submultiple 高斯消元求伯努利数