技術標籤：資訊競賽解題數學

題目

傳送門

題解

（由於本蒟蒻對富文字編輯器情有獨鍾，下面的推導公式可能都是以無法複製的圖片展示的）

題意很簡單，求，模數p是≥5e8的素數，所以放心求逆元；

下面開始推式子：

發現有個gcd，沒有好的方法直接快速求和，所以把gcd提出來：

後兩個的部分很熟悉，我們發現可以用莫比烏斯反演，設，（用m-大小f表示莫比烏斯反演推導的函式是為了與後面做區分），

感性分析可得，（定義），

所以有，

將帶入原式子：，

然後我們發現迴圈裡有個，已目前能力（是我太菜了）無法直接求和，除非把 i 變為某個數的因子形式，所以改變列舉順序，不列舉 d、i，而列舉 d*i 和 d；

令，有；

觀看尾部的，我們又雙叒發現有點熟悉，因為根據剛學過的狄利克雷卷積的知識，，也就是說，

這就很棒了，因為我們已經可以把帶的一堆換成尤拉函式：，

設，則；

顯然，最多隻有種取值，如果我們能較快求得 f 函式的字首和，就可以用數論分塊做；

所以我們用一用杜教篩：

設，，

；

根據套路，g 函式要根據 f 函式設定為最適合的一個完全積性函式，此處將 g 設定為，則，

由，也就是，可以得到；

由此我們得到了快速求 h 函式字首和的方法，直接套杜教篩：，

此時再用一個數論分塊，就可以較快求得 f 函式的字首和，

最後把兩個過程套起來，根據杜教篩的時間複雜度分析，全部平攤下來是的。

程式碼

#include<cstdio>
#include<cstring>
#include<iostream>
#include<algorithm>
#include<cmath>
#include<vector>
#include<queue>
#include<map>
#include<tr1/unordered_map>
#define ll long long
#define MAXN 10000005
using namespace std;
using namespace tr1;
inline ll read(){
	ll x=0;bool f=1;char s=getchar();
	while((s<'0'||s>'9')&&s>0){if(s=='-')f^=1;s=getchar();}
	while(s>='0'&&s<='9')x=(x<<1)+(x<<3)+s-'0',s=getchar();
	return f?x:-x;
}
ll n,p[MAXN],sp[MAXN],MOD=read();
unordered_map<ll,ll>ph;
bool nop[MAXN];
vector<ll>h;
inline ll ksm(ll a,ll b){
	ll res=1;
	for(;b;b>>=1){
		if(b&1)res=res*a%MOD;
		a=a*a%MOD;
	}
	return res;
}
ll NI=ksm(6,MOD-2);//提前求逆元
inline void build(){
	nop[0]=nop[1]=1;
	for(int i=1;i<MAXN-4;i++)p[i]=i;
	for(int a=2;a<MAXN-4;a++){
		if(!nop[a])h.push_back(a),p[a]=a-1;
		for(int i=0,u;i<h.size()&&h[i]*a<MAXN-4;i++){
			u=a*h[i],nop[u]=1;
			if(a%h[i]==0)p[u]=p[a]*h[i],i=MAXN;
			else p[u]=p[a]*p[h[i]];
		}
	}
	sp[1]=p[1];
	for(int i=2;i<MAXN-4;i++)sp[i]=(sp[i-1]+p[i]*i%MOD*i%MOD)%MOD;
}
inline ll sumx(ll n){n%=MOD;//n的範圍比模數大就離譜
	return (n*(n+1)>>1)%MOD;
}
inline ll sumy(ll n){n%=MOD;
	return n*(n+1)%MOD*(n*2+1)%MOD*NI%MOD;
}
inline ll sumphi(ll n){
	if(n<MAXN-4)return sp[n];
	if(ph.find(n)!=ph.end())return ph[n];
	ll res=sumx(n)*sumx(n)%MOD;
	for(ll i=2,ls;i<=n;i=ls+1)
		ls=n/(n/i),res=(res-(sumy(ls)-sumy(i-1)+MOD)%MOD*sumphi(n/i)%MOD+MOD)%MOD;
	return ph[n]=res;
}
inline ll getans(ll n){
	ll res=0;
	for(ll i=1,ls;i<=n;i=ls+1)
		ls=n/(n/i),res=(res+(sumphi(ls)-sumphi(i-1)+MOD)%MOD*sumx(n/i)%MOD*sumx(n/i)%MOD)%MOD;
	return res;
}
int main()
{
	build();n=read();
	printf("%lld\n",getans(n));
	return 0;
}