G tsy's number

推公式看这里
https://www.cnblogs.com/acjiumeng/p/10743919.html

推出来公式为
$\sum _{t=1}(t{)}^3\ast \frac{{\lfloor \frac{n}{t}\rfloor }^3\ast (\lfloor \frac{n}{t}\rfloor +1{)}^2\ast (\lfloor \frac{n}{t}\rfloor \ast 2+1)}{12}\sum _{d\mid t}\varphi (d)\ast \mu (\frac{t}{d})$t=1∑​(t)3∗12⌊tn​⌋3∗(⌊tn​⌋+1)2∗(⌊tn​⌋∗2+1)​d∣t∑​ϕ(d)∗μ(dt​)

但是这位博主的代码实现感觉不太好，我综合了群里大佬的各种写法，写出了一下代码

1. 求和公式最后出现了12，而模数是1<<30 ,这样不好求，我们发现直接求 $sum\left(i^2\right),sum\left(i^3\right)$sum(i2),sum(i3)是可以接受的，这样就避免了除法
2. 关于杜教筛求解积性函数
   bzoj 4804
   $f\left(n\right)=\sum _{d\mid t}\varphi \left(d\right)\ast \mu \left(\frac{t}{d}\right)$f(n)=d∣t∑​ϕ(d)∗μ(dt​)
   $f\left(p^k\right)=\varphi \left(p^k\right)\ast u\left(1\right)+\varphi \left(p^{(k-1)}\right)\ast u\left(p\right)=\varphi \left(p^k\right)-\varphi \left(p^{(k-1)}\right)$f(pk)=ϕ(pk)∗u(1)+ϕ(p(k−1))∗u(p)=ϕ(pk)−ϕ(p(k−1))
   当 $k==1$k==1
   $f\left(p\right)=p-2$f(p)=p−2
   当 $k\&gt;=2$k>=2
   $f\left(p^k\right)=p^{(k-2)}\ast (p-1{)}^2$f(pk)=p(k−2)∗(p−1)2
   假设 $n=n^{\prime }\ast p^k$n=n′∗pk p 为k的最小的素因子k = 1时
   $f\left[n\right]=f\left[n\right]\ast f\left[p\right]$f[n]=f[n]∗f[p]
   当 $k==2$k==2
   $f\left[n\right]=f\left[n/p/p\right]\ast (p-1)\ast (p-1)$f[n]=f[n/p/p]∗(p−1)∗(p−1)
   当 $k\&gt;=3$k>=3
   $f\left[n\right]=f\left[n/p\right]\ast p$f[n]=f[n/p]∗p
   根据线性筛的性质，一个数一定是被它最小的素因子筛到，于是我们便可得到代码

#include <bits/stdc++.h>
using namespace std;
typedef long long LL;
const LL     mod = 1<<30;
const int maxn = 1e7+10;
int  pow1[maxn],pow2[maxn];//,pow3[maxn];//,pow4[maxn];
// \sum{d|n} phi*u
LL prime[maxn],cnt,f[maxn];
bool vis[maxn];
void init(int n){
	// pow2[i] = 
	for(LL i = 1;i <= n; ++i){
		pow1[i] = pow1[i-1]+i,pow1[i] %= mod;
		pow2[i] = pow2[i-1]+i*i%mod,pow2[i] %= mod;
	}
	f[1] = 1;
	for(int i = 2;i <= n; ++i){
		if(!vis[i]) prime[cnt++] = i,f[i] = i-2;
		for(int j = 0;j < cnt; ++j){
			if(prime[j] * i > n)break;
			vis[prime[j]*i] = 1;
			if(i%prime[j] == 0){

				if(i/prime[j]%prime[j] == 0)
					f[i*prime[j]] = f[i]*prime[j]%mod;
				else
					f[i*prime[j]] = f[i/prime[j]]*(prime[j]-1)%mod*(prime[j]-1)%mod;
				// if(i*prime[j]== 4)
				// cout<<f[i*prime[j]]<<" "<<prime[j]<<endl;
				break;
			}
			else
				f[i*prime[j]] = f[i]*f[prime[j]]%mod;

		}
	}
	// cout<<f[4]<<endl;
	for(int i = 1;i <= n; ++i){
		f[i] = 1ll*i*i%mod*i%mod*f[i]%mod;
		f[i] = (f[i-1]+f[i])%mod;
	}
}
LL solve(LL n){
	LL ans = 0;
	for(int i = 1;i <= n; ){
		LL t = n/i;
		LL t2= min(n/t,n);
		// n/t*pow2[n/t]*pow3[n/t]
		ans = (ans + (f[t2]-f[i-1])%mod*t%mod*pow1[t]%mod*pow2[t]%mod)%mod;
		i = t2+1;
	}
	return ans;
}
int main(void){
	init(maxn-1);
	// cout<<f[4]<<endl;
	int t;cin>>t;
	while(t--){
		LL n;scanf("%lld",&n);
		printf("%lld\n",(solve(n)%mod+mod)%mod);
	}
	return 0;
}