不会……
两篇加在一起都看不懂……
https://www.cnblogs.com/cellular-automaton/p/8241128.html
https://www.luogu.org/blog/cjyyb/solution-p3768
1 //minamoto 2 #include<iostream> 3 #include<cstdio> 4 #include<map> 5 #define ll long long 6 using namespace std; 7 const int N=4e6+5; 8 map<ll,ll> mp; 9 ll sum[N],phi[N],pri[N],cnt,vis[N],mod6,mod; 10 //mod是2的逆元，mod6是6的逆元 11 void init(ll p){ 12 vis[1]=sum[1]=1; 13 for(int i=2;i<N;++i){ 14 if(!vis[i]) pri[++cnt]=i,phi[i]=i-1,sum[i]=1ll*i*i%p*phi[i]%p; 15 for(int j=1;j<=cnt,i*pri[j]<N;++j){ 16 vis[i*pri[j]]=1; 17 int now=i*pri[j]; 18 if(i%pri[j]==0){ 19 phi[now]=phi[i]*pri[j]%p; 20 sum[now]=phi[now]*now%p*now%p; 21 break; 22 } 23 phi[now]=phi[i]*(pri[j]-1)%p; 24 sum[now]=phi[now]*now%p*now%p; 25 } 26 } 27 for(int i=2;i<N;++i) (sum[i]+=sum[i-1])%=p; 28 } 29 ll ksm(ll x,ll y,ll p){ 30 ll res=1; 31 while(y){ 32 if(y&1) (res*=x)%=p; 33 (x*=x)%=p,y>>=1; 34 } 35 return res; 36 } 37 ll calc(ll n,ll p){ 38 n%=p; 39 ll res=((1+n)*n%p)*mod%p; 40 (res*=res)%=p; 41 return res; 42 } 43 ll calcs(ll n,ll p){ 44 n%=p; 45 ll res=(n*(n+1)%p)*(2*n+1)%p; 46 (res*=mod6)%=p; 47 return res; 48 } 49 ll calcsum(ll n,ll p){ 50 if(n<N) return sum[n]; 51 if(mp.count(n)) return mp[n]; 52 ll x=2,res=calc(n,p); 53 while(x<=n){ 54 ll y=n/(n/x); 55 res=((res-(calcs(y,p)-calcs(x-1,p)+p)%p*calcsum(n/x,p)%p)+p)%p; 56 x=y+1; 57 } 58 return mp[n]=res; 59 } 60 int main(){ 61 // freopen("testdata.in","r",stdin); 62 ll p,n; 63 scanf("%lld%lld",&p,&n); 64 mod=ksm(2,p-2,p),mod6=ksm(6,p-2,p); 65 init(p); 66 ll x=1,ans=0; 67 while(x<=n){ 68 ll y=n/(n/x); 69 (ans+=((calcsum(y,p)-calcsum(x-1,p)+p)%p*calc(n/x,p)%p))%=p; 70 x=y+1; 71 } 72 printf("%lld ",ans); 73 return 0; 74 }