题目描述
定义积性函数f(x)f(x),且f(p^k)=p^k(p^k-1) (p是一个质数),求
f(i) 的前缀和
对10^9+7109+7取模。
输入格式
一行一个整数nn。
输出格式
一个整数表示答案。
输入输出样例
输入 #1复制
10
输出 #1复制
263
输入 #2复制
1000000000
输出 #2复制
710164413
说明/提示
f(1)=1,f(2)=2,f(3)=6,f(4)=12,f(5)=20f(1)=1,f(2)=2,f(3)=6,f(4)=12,f(5)=20
f(6)=12,f(7)=42,f(8)=56,f(9)=72,f(10)=40f(6)=12,f(7)=42,f(8)=56,f(9)=72,f(10)=40
SOLUTION:
洛谷的题解和网上的题解中求S()的时候是不同的
我用的是网上的方法
https://www.cnblogs.com/zhoushuyu/p/9187319.html
https://www.luogu.org/problemnew/solution/P5325
CODE:
#include<cstdio>
#include<algorithm>
#include<cstring>
#include<cmath>
#define ll long long
#define int long long
using namespace std;
const ll MOD=1000000007,inv2=500000004,inv3=333333336;
ll prime[1000005],num,sp1[1000005],sp2[1000005];
ll n,Sqr,tot,g1[1000005],g2[1000005],w[1000005],ind1[1000005],ind2[1000005];
bool flag[1000005];
void pre(int n)//预处理,线性筛
{
flag[1]=1;
for(int i=1;i<=n;i++)
{
if(!flag[i])
{
prime[++num]=i;
sp1[num]=(sp1[num-1]+i)%MOD;
sp2[num]=(sp2[num-1]+1ll*i*i)%MOD;
}
for(int j=1;j<=num&&prime[j]*i<=n;j++)
{
flag[i*prime[j]]=1;
if(i%prime[j]==0)break;
}
}
}
ll S(ll x,int y)//第二部分
{
if(x<=1||prime[y]>x)return 0;
ll k=x<=Sqr?ind1[x]:ind2[n/x];
ll ans=(g2[k]-g1[k]+MOD-(sp2[y-1]-sp1[y-1])+MOD)%MOD;
for(int i=y;i<=num&&prime[i]*prime[i]<=x;i++)
{
ll p1=prime[i],p2=1ll*prime[i]*prime[i];
for (int e=1;p2<=x;++e,p1=p2,p2*=prime[i])
ans+=((1ll*S(x/p1,i+1)*p1%MOD*(p1-1)%MOD)%MOD+1ll*p2%MOD*((p2-1)%MOD)%MOD),
ans%=MOD;
/*
ll pe=prime[i];
for(int e=1;pe<=x;e++,pe=pe*prime[i])
{
ll xx=pe%MOD;
ans=(ans+xx*(xx-1)%MOD*(S(x/pe,i)+(e!=1)))%MOD;
}*/
}
return ans%MOD;
}
signed main()
{
scanf("%lld",&n);
Sqr=sqrt(n);
pre(Sqr);
for(ll i=1;i<=n;)
{
ll j=n/(n/i);
w[++tot]=n/i;
g1[tot]=w[tot]%MOD;
g2[tot]=g1[tot]*(g1[tot]+1)/2%MOD*(2*g1[tot]+1)%MOD*inv3%MOD;
g2[tot]--;
g1[tot]=g1[tot]*(g1[tot]+1)/2%MOD-1;
if(n/i<=Sqr)ind1[n/i]=tot;
else ind2[n/(n/i)]=tot;
i=j+1;
}//g1,g2分别表示一次项和二次项,ind1和ind2用来记录这个数在数组中的位置
for(int i=1;i<=num;i++)//由于g数组可以滚动,所以就只开了一维
{
for(int j=1;j<=tot&&prime[i]*prime[i]<=w[j];j++)
{
ll k=w[j]/prime[i]<=Sqr?ind1[w[j]/prime[i]]:ind2[n/(w[j]/prime[i])];
g1[j]-=prime[i]*(g1[k]-sp1[i-1]+MOD)%MOD;
g2[j]-=prime[i]*prime[i]%MOD*(g2[k]-sp2[i-1]+MOD)%MOD;
g1[j]%=MOD,g2[j]%=MOD;
if(g1[j]<0)g1[j]+=MOD;
if(g2[j]<0)g2[j]+=MOD;
}
}
printf("%lld\n",(S(n,1)+1)%MOD);
return 0;
}
对于30\%30%的数据,保证1\le n\le 10^61≤n≤106。
对于100\%100%的数据,保证1\le n\le 10^{10}1≤n≤1010