题目大意:定义 fk(n)=∑1≤i≤n,(i,n)=1ik,求fk(∏wi=1paii) f k ( n ) = ∑ 1 ≤ i ≤ n , ( i , n ) = 1 i k , 求 f k ( ∏ i = 1 w p i a i ) ,其中 k≤1e2, w≤1e3 k ≤ 1 e 2 , w ≤ 1 e 3 。
题解(前半部分在胡扯可以跳过):
首先互质不好限制,考虑从总数中减去。
注意到涉及的状态只会是n的因数,那么对于n<=1e9是随便跑的,但这个题却不行。
首先要解决的第一个问题是如何求 ∑ni=1ik ∑ i = 1 n i k ,其中给出n的因式分解;
把这个东西拉格朗日差值出来,高斯消元一下求出系数表示(反正只有100这么大)。
这样因为要对大质数取模所以直接把n取模即可O(wlgw)的时间对任意n求出答案。
但是没有办法处理后半部分。
向莫比乌斯反演方向思索:
注意到后面的求和是一个k+1次多项式,显然可以高斯消元求出它的系数 {ak+1} { a k + 1 } ,然后可以考虑对g的每一项求和;
后面那个 gi(n) g i ( n ) 显然是一个积性函数啦!
因此枚举次数之后就可以对每一个 pajj p j a j 单独计算答案:
gi(pc)=(pc)i−pk(pc−1)i g i ( p c ) = ( p c ) i − p k ( p c − 1 ) i ,然后乘起来即可。
显然g函数可以log时间求出,这样问题再 O(k3+kwlgc) O ( k 3 + k w l g c ) 级别内做出这个题。
代码:
#include<iostream>
#include<cstring>
#include<cstdio>
#include<algorithm>
#define N 110
#define K N
#define W 1010
#define gc getchar()
#define lint long long
#define inv(x) fast_pow(x,mod-2)
#define mod 1000000007
#define debug(x) cerr<<#x<<"="<<x
#define sp <<" "
#define ln <<endl
using namespace std;
inline int inn()
{
int x,ch;while((ch=gc)<'0'||ch>'9');
x=ch^'0';while((ch=gc)>='0'&&ch<='9')
x=(x<<1)+(x<<3)+(ch^'0');return x;
}
int k,w,a[N][N],b[N],p[W],c[W];
inline int fast_pow(int x,int k,int ans=1)
{ for(;k;k>>=1,x=(lint)x*x%mod) (k&1)?ans=(lint)ans*x%mod:0;return ans; }
inline int f(int n)
{
int ans=0;
for(int i=1;i<=n;i++)
(ans+=fast_pow(i,k))%=mod;
return ans;
}
inline int g(int i,int p,int c)
{
return (fast_pow(fast_pow(p,c),i)-(lint)fast_pow(p,k)*fast_pow(fast_pow(p,c-1),i)%mod+mod)%mod;
}
inline int Gauss(int n)
{
for(int i=1;i<=n;i++)
{
int x=inv(a[i][i]);b[i]=(lint)b[i]*x%mod;
for(int j=i;j<=n;j++) a[i][j]=(lint)a[i][j]*x%mod;
for(int j=1,y;j<=n;j++)
if(i^j)
{
y=a[j][i],b[j]=(b[j]-(lint)y*b[i]%mod+mod)%mod;
for(int k=i;k<=n;k++)
a[j][k]=(a[j][k]-(lint)y*a[i][k]%mod+mod)%mod;
}
}
return 0;
}
int main()
{
k=inn(),w=inn();
for(int i=1;i<=w;i++) p[i]=inn(),c[i]=inn();
for(int i=1;i<=k+2;b[i]=f(i),i++)
for(int j=1;j<=k+2;j++)
a[i][j]=fast_pow(i,j-1);
Gauss(k+2);int ans=0;
for(int i=0;i<=k+1;i++)
{
int res=b[i+1];
for(int j=1;j<=w;j++)
res=(lint)res*g(i,p[j],c[j])%mod;
(ans+=res)%=mod;
}
return !printf("%d\n",ans);
}