一开始发现f(xy)=f(x)+f(y),这不是积性函数能求?
思路参考:https://www.cnblogs.com/hua-dong/p/11508790.html
质数p对sum{i=1~n}f(i)的贡献是sum{k=1,p^k<=n}(n/p^k)
就像阶乘算多少个0,
5倍数多1个0
25的倍数多两个0
但算25的倍数时5的倍数已经对25的倍数做出过贡献,所以减掉之后25的贡献为1
这里f(i)也是一样,p^1的贡献是1 p^2贡献是2
但算p的倍数时p的倍数已经对p^2的倍数做出过贡献,所以减掉之后p^2的贡献也是1
也可以这样想,假设n=16 ,p=2
2对至少有1个质因子2的数贡献 =n/2
2对至少为2个质因子2的数贡献 =n/4
所以2对仅有有1个质因子2的数 (f(2),f(6),f(10),f(14)) 的贡献 =(n/2-n/4)*1=4
4的贡献 f(4)+f(12) =(n/4-n/8)*2
8的贡献 f(8) =(n/8-n/16)*3
16的贡献f(16)=n/16*4;
4个加起来 发现:ans= n/2+n/4+n/8+n/16
对于i*f(i)
2的贡献 f(2)+f(6)+f(10)+f(14) = 1*( (2+4+6+8+...16)-(4+8+12+16) )
=1*( 2*(1+2+3.....8) - 4*(1+2+3+4) ) =1*(2sum(8)-4sum(4) )
4的贡献 f(4)+f(12)=2*(4sum(4)-8sum(2))
8的贡献 f(8)=3*(8sum(4)-16sum(1))
16的贡献 f(16)=4*(16sum(1))
4个加起来 发现:ans= 2*sum(n/2)+4*sum(n/4)+8*sum(n/8)+16*sum(n/16);
所以我们答案= n以内每个质数p的贡献和
对于p<=sqrt(n)的质数,我们枚举质数k,质数的倍数t,暴力把p,p^2,p^3的贡献把柿子求出来
对于p>sqrt(n)的质数,n/p<sqrt(n),显然有t=1,那么
对于前半部分,我们设
众所周知,在min25筛中,对于所有可能的n/i=m的值 g(m)是可以套模板维护出来的(整数分块)
因此 [w[m],w[m+1] ]区间函数和也是很快能维护出来的=g(m)-g(m+1)
比如n=10 那么g(10) g(5) g(3) g(2) g(1)能求出来 w[1]=10 w[2]=5 w[3]=3....
比如n=10 那么[6,10]=g[1]-g[2] , [4,5]=g(2)-g(3) 的和都能维护出来
那么我们枚举所有(n/i)=t<sqrt(n)的值 对于每个t,看有多少个n/质数在[t,t+1]内,这区间内每个质数的贡献是t =n/质数
因此ans+=(n+1)*(g[t]-g[t+1])*t
对于后半部分:我们设
同理min25 可求g(m)
因此ans-= (g[t]-g[t+1])*sum(t)=(1+t)*t/2 * (g[t]-g[t+1])
436ms
//min25筛求f前缀和
//要求:积性函数+f(p)是多项式(或者是完全积性函数的和 )+ f(p^c)非常好求
//51nod 1239 欧拉函数前缀和
#include<stdio.h>
#include<math.h>
#include<iostream>
#include<string.h>
using namespace std;
#define ll long long
const ll mod = 998244353;
const ll inv2=(1+mod)/2;
ll qmod(ll a,ll b){ll res=1;while(b){if(b&1)res=res*a%mod;a=a*a%mod;b=b>>1;}return res;}
namespace Min25{
const int N = 1e6+5;
ll Sqr,isp[N],pri[N],sump[N],tot=0;//质数上限
void init(int n){//线筛预处理质数 与质数的前缀和(f(p)的前缀和)
isp[1]=1;
for (ll i=2;i<=n;++i){
if (!isp[i]) pri[++tot]=i,sump[tot]=1ll*(sump[tot-1]+i)%mod;
for (int j=1;i*pri[j]<=n;++j){
isp[i*pri[j]]=1;
if (i%pri[j]==0) break;
}
}
}
ll w[N],id1[N],id2[N],m=0,g[N],h[N];
void getg0(ll n){//对每个n/i 求 g(n/i,0),h(n/i,0);
//定义: g(x,p)=sum{i=1~x}[i的最小质因子>pri[p] or i是质数]g0(i) ,要求g0(p)是完全积性函数,不是的话要分开凑。
//g(n/i,0)代表 所有数当成质数带入 本来g0(p)=? 仅对质数有效
//g[1]=g(n/1,0)....g[m]=g(n/m,0) m是整数分块的上限
memset(h,0,sizeof(h));
m=0;Sqr=sqrt(n);
for (ll i=1,j;i<=n;i=j+1){
j=n/(n/i);w[++m]=n/i;//整数分块 ; 下标映射 ;
if (w[m]<=Sqr) id1[w[m]]=m;
else id2[n/w[m]]=m;
//合数也当成质数求和,忽略h0(1) h0(i)=1 g0(i)=i 需要改
h[m]=(w[m]-1+mod)%mod;
g[m]=((w[m]+2)%mod)*((w[m]-1)%mod)%mod;
if (g[m]&1) g[m]+=mod;g[m]/=2;
}
}
void getgp(ll n){//对每个n/i 求g(n/i,|P|) h(n/i,|P|)
// 定义:g(x,p)=sum{i=1~x}[i的最小质因子>p or i是质数]g0(i) ,其中g0(i)=i; 维护时省略第二维 f(p)=g0(p)-h0(p)凑完全积性
// 定义:h(x,p)=sum{i=1~x}[i的最小质因子>p or i是质数]h0(i) ,其中h0(i)=1; 维护时省略第二维
// 所以我们这个函数求的是 g(x,|P|)=sum{i=1~x}[i是质数]g0(i)
// g[1]=g(n/1,|P|) g[2]=g(n/2,|P|) ...g[m]=g(n/m,|P|) m是整数分块的上限
for (int j=1;j<=tot;++j)
for (int i=1;i<=m&&1ll*pri[j]*pri[j]<=w[i];++i){
int k=(w[i]/pri[j]<=Sqr)?id1[w[i]/pri[j]]:id2[n/(w[i]/pri[j])];//不用改
g[i]=(g[i]-1ll*pri[j]*(g[k]-sump[j-1])%mod)%mod;g[i]=(g[i]+mod)%mod;//需要改 g(n,j)= g(n,j-1) - g0(pri[j])* [ g(k,j-1)-sum_{i=1~j-1}g0(pri[i]) ];
h[i]=(h[i]-h[k]+j-1);h[i]=(h[i]+mod)%mod; //需要改 同上
}
//g[i]-g[i+1] 代表[n/(i+1),n/i]区间内的质数p的g0(p)的和
}
int fpc(int p,int c){//f(p^c) 需要改
return (qmod(p,c)-qmod(p,c-1)+mod)%mod;
}
int S(ll x,int y,ll n){//定义:S(x,y)=sum{i=1~x}[i的最小质因子>=pri[y] ]*f(i); 求积性函数前缀和ans=S(n,1,n)+f(1)
//在此之前,f(n/i,|P|) (n/i)<=m 的各项已经用getgp()维护出来了
if (x<=1||pri[y]>x) return 0;
int k=(x<=Sqr)?id1[x]:id2[n/x];
ll res=(1ll*g[k]-h[k]-sump[y-1]+y-1);(res=res+mod)%mod;// res= g(n,|P|) - sum_{i=1~j}f(Pi) 需要改
//if (y==1) res+=2;//特判 f(2)=f(p)+2 其他题可删
for (int i=y;i<=tot&&1ll*pri[i]*pri[i]<=x;++i){//S的递推式 不改
ll p1=pri[i],p2=1ll*pri[i]*pri[i];
for (int e=1;p2<=x;++e,p1=p2,p2*=pri[i])
(res+=(1ll*S(x/p1,i+1,n)*fpc(pri[i],e)%mod+fpc(pri[i],(e+1))))%=mod;//pri[i]^e 是f(pri[i]^e)的值,需要改
}
return res%mod;
}
}
int main(){
ll n;
using namespace Min25;
scanf("%lld",&n);
Sqr=sqrt(n);init(Sqr);
getg0(n);getgp(n);
/*
for(int i=1;i<=n;i++){
printf("%d %d %d\n",n/i,i,g[i]);
}*/
ll ans=0;
for(int i=1;i<=tot;i++){//p<sqrt(n)的贡献
for(ll j=pri[i];j<=n;j*=pri[i]){
ll t=(n/j)%mod;
(ans+=t*((n+1)%mod)%mod )%=mod;
(ans-=(j%mod)*(mod+(1+t)%mod*t%mod*inv2%mod)%mod)%=mod;
ans=(ans+mod)%mod;
}
}
cout<<ans<<endl;
for(ll k=1;k<Sqr;k++){//p>sqrt(n) 枚举n/p=k的值
(ans+=((n+1)%mod)*k%mod*(h[k]-h[k+1]+mod)%mod)%=mod;
(ans-=mod+(1+k)*k%mod*inv2%mod*(g[k]-g[k+1]+mod)%mod)%=mod;
ans=(ans+mod)%mod;
}
cout<<(ans+mod)%mod<<endl;
//system("pause");
return 0;
}