【題目地址】

題意簡述

給定你$n,m$求：

$\sum_{i_1=1}^m\sum_{i_2=1}^m\cdots\sum_{i_n=1}^mgcd(i_1,i_2,\cdots,i_n)$

對$2^{64}$取模，$n,m\leq 10^{11}$

---

這個題和CQOI2015選數很像，但是$10^{11}$顯然要杜教篩，所以我們按照套路反演一波，列舉$gcd$，原式轉化為：

$\sum_{d=1}^md\sum_{i_1=1}^m\sum_{i_2=1}^m\cdots\sum_{i_n=1}^m[gcd(i_1,i_2,\cdots,i_n)=d] \\ \sum_{d=1}^md\sum_{i_1=1}^{\lfloor\frac{m}{d}\rfloor}\sum_{i_2=1}^{\lfloor\frac{m}{d}\rfloor}\cdots\sum_{i_n=1}^{\lfloor\frac{m}{d}\rfloor}[gcd(i_1,i_2,\cdots,i_n)=1]$

將$\mu$放入：

$\sum_{d=1}^md\sum_{i_1=1}^{\lfloor\frac{m}{d}\rfloor}\sum_{i_2=1}^{\lfloor\frac{m}{d}\rfloor}\cdots\sum_{i_n=1}^{\lfloor\frac{m}{d}\rfloor}\sum_{w|i_1,w|i_2,\cdots,w|i_n}\mu(w) \\ \sum_{d=1}^md\sum_{w=1}^{\lfloor\frac{m}{d}\rfloor}\mu(w)\left\lfloor\frac{m}{dw}\right\rfloor^n$

根據套路，我們列舉$dw$，令$T=dw$得到：

$\sum_{T=1}^m\left(\left\lfloor\frac{m}{T}\right\rfloor\right)^n\sum_{d|T}\mu(d)\frac{T}{d}$

後面部分就是$\mu * id=\varphi$，所以原式變為：

$\sum_{T=1}^m\left(\left\lfloor\frac{m}{T}\right\rfloor\right)^n\varphi(T)$

所以杜教篩加快速冪即可，而模數十分特殊，所以直接用unsigned long long，但是對於$\frac{n\times(n+1)}{2}$我們不再用逆元，而是直接判奇偶除過去。

$\sum_{i=1}^n\varphi(i)=\sum_{i=1}^ni-\sum_{i=2}^n\sum_{j=1}^{\lfloor\frac{n}{i}\rfloor}\varphi(j)$

上程式碼：

#include<map>
#include<cstdio>
#include<cstring>
#include<iostream>
#include<algorithm>
#define ll unsigned long long
using namespace std;
const int M=1e7+1;
ll n,m;
ll fpow(ll a,ll b){
	ll res=1llu;
	for(;b;b>>=1,a*=a)if(b&1)res*=a;
	return res;
}
ll prime[M],cnt,phi[M];
bool vis[M];
void init(){
	phi[1]=1llu;
	for(ll i=2;i<M;i++){
		if(!vis[i]){
			prime[++cnt]=i;
			phi[i]=i-1llu;
		}
		for(ll j=1,v;j<=cnt&&i*prime[j]<M;j++){
			v=i*prime[j];
			vis[v]=1;
			if(!(i%prime[j])){
				phi[v]=phi[i]*prime[j];
				break;
			}
			phi[v]=phi[i]*phi[prime[j]];
		}
	}
	for(ll i=1;i<M;i++)phi[i]+=phi[i-1];
}
ll S(ll x){
	if(x&1) return (x+1)/2*x;
	else return x/2*(x+1);
}
map <ll,ll> mp;
ll calc(ll x){
	if(x<M) return phi[x];
	if(mp.count(x)) return mp[x];
	ll ans=S(x);
	for(ll i=2,j;i<=x;i=j+1){
		j=(x/(x/i));
		ans-=(j-i+1)*calc(x/i);
	}
	return mp[x]=ans;
}
ll solve(ll x,ll y){
	ll ans=0;
	for(ll i=1,j;i<=x;i=j+1){
		j=(x/(x/i));
		ans+=(calc(j)-calc(i-1))*fpow(x/i,y);
	}
	return ans;
}
int main(){
	cin>>n>>m;
	init();
	cout<<solve(m,n)<<endl;
	return 0;
}

---

End