简单的数学题

题意:

求(sum_{i=1}^{n}sum_{j=1}^{n}{ijgcd(i,j)},n<=1e10)

---

题解

先把式子化简一下
(sum_{i=1}^{n}sum_{j=1}^{n}{ijgcd(i,j)}\sum_{d=1}^{n}sum_{i=1}^{n}sum_{j=1}^{n}{ijgcd(i,j)}\sum_{d=1}^{n}{d^3}sum_{i=1}^{n/d}sum_{j=1}^{n/d}{ij[gcd(i,j)=1]})
然后还是设(F(t)=sum_{i=1}^{n/d}sum_{j=1}^{n/d}{ij[t|gcd(i,j)]})
(F(t)=t^2sum_{i=1}^{n/dt}sum_{j=1}^{n/dt}{ij})
(B(n)=sum_{i=1}^{n}i)
所以原式就是(sum_{d=1}^{n}{d^3}{f(1)})
(sum_{d=1}^{n}{d^3}sum_{i=1}^{n/d}{mu(i)i^2B(frac{n}{id})^2})
还是老套路设(T=id)
(sum_{T=1}^{n}{B(frac{n}{T})^2}sum_{d|T}{mu(frac{T}{d})(frac{T}{d})^2d^3})
(sum_{T=1}^{n}{B(frac{n}{T})^2}T^2sum_{d|T}{mu(frac{T}{d})d})
然后如果(n<=1e7)直接线筛就可以做了
但是杜教筛似乎不好筛以狄利克雷卷积的形式定义的函数
但是看后面的东西我们想到了什么?
(mu*id=phi)
所以后面那个东西就是(phi(T))
所以我们就只需要去筛(f(T)=T^2phi(T)了)
那么我们考虑用什么作为(g)函数
(sum_{d|T}{phi(d)d^2g(frac{T}{d})})
似乎(g=id^2)就很好
这样就成了(sum_{d|T}{phi(d)id(T)^2})
然后把(id(T)^2)提前
就成了(id(T)^2sum_{d|T}{phi(d)}\=id(T)^2id(T)\=id(T)^3)
这样就可以筛了

代码

#include<map>
#include<cstdio>
#include<cstring>
#include<iostream>
#include<algorithm>
# define LL long long
const int M = 4000005 ;
using namespace std ;

bool notp[M] ;
int pnum , p[M] , upp = 4000000 ;
LL mod , n , sum[M] , ans , inv2 , inv6 ;
map < LL , LL > psum ;

inline LL dc(LL n) { 
	n %= mod ;
	return n % mod * (n + 1) % mod * inv2 % mod; 
}
inline LL dc2(LL n) { 
	n %= mod ;
	return n % mod * (n + 1) % mod * (n * 2 + 1) % mod * inv6 % mod ; 
}
inline LL Fpw(LL Base , LL k) {
	LL temp = 1 ; Base %= mod ; k %= (mod - 1) ;
	while(k) {
		if(k & 1) temp = temp * Base % mod ;
		Base = Base * Base % mod ; k >>= 1 ;
	}
	return temp ;
}
inline void Presolve(int n) {
	sum[1] = 1 ; 
	for(int i = 2 ; i <= n ; i ++) {
		if(!notp[i]) {
			p[++pnum] = i ;
			sum[i] = (i - 1 + mod) % mod ;
		}
		for(int j = 1 ; j <= pnum && 1LL * i * p[j] <= n ; j ++) {
			notp[i * p[j]] = true ;
			if(i % p[j] == 0) sum[i * p[j]] = sum[i] * p[j] % mod ;
			else sum[i * p[j]] = sum[i] * sum[p[j]] % mod ;
		}
	}
	for(int i = 1 ; i <= n ; i ++)
		sum[i] = (sum[i - 1] + sum[i] * i % mod * i % mod) % mod ;
}

inline LL Sum(LL n) {
	if(n <= upp) return sum[n] ;
	if(psum[n]) return psum[n] ;
	LL ret = dc(n) * dc(n) % mod ;
	for(LL l = 2 , r ; l <= n ; l = r + 1) {
		r = n / (n / l) ;
		LL tp = ((dc2(r) - dc2(l - 1)) % mod + mod) % mod ;
		ret = ((ret - tp * Sum(n / l) % mod) % mod + mod) % mod ;
	}
	ret %= mod ;
	psum[n] = ret ; return ret ;
}
int main() {
	scanf("%lld%lld",&mod,&n) ;
	inv2 = Fpw(2 , mod - 2) ;
	inv6 = Fpw(6 , mod - 2) ;
	Presolve(min(n , 4000000LL)) ;
	for(LL l = 1 , r ; l <= n ; l = r + 1) {
		r = n / (n / l) ;
		LL tp = ((Sum(r) - Sum(l - 1)) % mod + mod) % mod ;
		ans = (ans + (tp * dc(n / l) % mod * dc(n / l) % mod + mod) % mod + mod) % mod ;
	}
	printf("%lld
",(ans % mod + mod) % mod) ;
	return 0 ;
}