求
∏i=1n(∑j=1igcd(i,j))mod109+7(n≤5∗107)
本来想出一道莫比乌斯反演的题目的,但是
∑太简单了,就换成了
∏试试看,结果发现一些奇怪的东西…
定义
f(n)=∑i=1ngcd(i,n),则有:
{f(pk)=1(pk−pk−1)+p(pk−1−pk−2)...+pk=k(pk−pk−1)+pk(枚举gcd大小并求出次数)f(n)=f(p)f(q)(p⊥q)
第二个式子是证明它是积性函数的:
令
n=pcqd(p⊥q),则需要证明
∑i=1ngcd(i,n)=∑j=1pcgcd(pc,j)∑k=1qdgcd(qd,k).
此时我们计算一下
paqb(a≤c,b≤d)作为gcd的出现次数.
从左边看则为
paqbn−⌊pa+1qbn⌋−⌊paqb+1n⌋+⌊pa+1qb+1n⌋.
从右边看则为
(pc−a−⌊pc−a−1⌋)(qd−b−⌊qd−b−1⌋).
很明显,拆开后是等价的.
那么我们只要求出
f(pk)及其出现次数,用一下快速幂即可.
本算法的瓶颈在于线性筛,因为快速幂是高度压缩了运算次数.
#include<cstdio>
#include<cstring>
#include<algorithm>
using namespace std;
typedef long long ll;
const int mod=1e9+7,N=5e7+10;
ll ans;
ll power_mod(ll a,ll b) {
ll c=1;a%=mod;
while(b) {
if(b&1)c=c*a%mod;
a=a*a%mod;b=b>>1;
}
return c;
}
void calc(int x,ll y,ll p) {
int cnt=0;
do {
x/=p;cnt++;
ans=ans*power_mod(y*p+cnt*y*(p-1),x-x/p)%mod;
y*=p;
}while(x>=p);
}
int n,prime[3001144],tot;bool v[N];
void get_prime() {
for(int i=2;i<=n;i++) {
if(!v[i])prime[++tot]=i,calc(n,1,i);
for(int j=1;i*prime[j]<=n;j++) {
v[i*prime[j]]=1;
if(i%prime[j]==0)break;
}
}
}
int main() {
scanf("%d",&n);ans=1;
get_prime();printf("%lld\n",ans);
return 0;
}