简单的数学题

题意:

求\(\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^2\sum_{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^2\sum_{d|T}{\mu(\frac{T}{d})d}\)
然后如果\(n<=1e7\)直接线筛就可以做了
但是杜教筛似乎不好筛以狄利克雷卷积的形式定义的函数
但是看后面的东西我们想到了什么?
\(\mu*id=\phi\)
所以后面那个东西就是\(\phi(T)\)
所以我们就只需要去筛\(f(T)=T^2\phi(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)^2\sum_{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\n",(ans % mod + mod) % mod) ;
	return 0 ;
}