RT，复杂度应该是对的。

#include<iostream>
#include<cstring>
#include<cstdio>
#include<cmath>
#include<algorithm>
#include<map>
#pragma comment(linker,"/stack:200000000")
#pragma GCC target("sse,sse2,sse3,ssse3,sse4.1,sse4.2,avx,avx2,popcnt,tune=native")
using namespace std;
typedef __int128 ll;
ll MOD=1;
ll n[10005],k[10005],t,MX;
ll e[2000005];
map<ll,ll> szs;
int b[12000005],p[35][5],GG,cur,pi[2000005],d[2000005],sz[5],mrk_cnt,r,K,block=2000000;
bool a[2000005];
void shai(int n){
	a[0]=a[1]=true;
	e[1]=1;
	for(int i=2;i<=block;i++){
		if(!a[i]){
			b[++r]=i;
			d[i]=1;
			e[i]=4;
		}
		pi[i]=r;
		for(int j=1;j<=r&&i*b[j]<=block;j++){
			a[i*b[j]]=true;
			d[i*b[j]]=1;
			e[i*b[j]]=e[i]*e[b[j]];
			if(i%b[j]==0){
				d[i*b[j]]=d[i]+1;
				e[i*b[j]]=e[i]/(d[i]*3+1)*(d[i*b[j]]*3+1);
				break;
			}
		}
	}
	for(int i=1;i<=block;i++){
		e[i]+=e[i-1];
	}
}
inline void write(__int128 x){
    if(x>9)
        write(x/10);
    putchar(x%10+'0');
}
ll getphi(ll x,int s){
	if(!s){
		return x;
	}
	if(s<=2){
		return p[x%sz[s]][s]+(x/sz[s])*p[sz[s]][s];
	}
	if(x<=b[s]*b[s]){
		return pi[x]-s+1;
	}
	if(x<=b[s]*b[s]*b[s]&&x<9000){
		int sx=pi[int(pow(x,1.0/2.0))];
        ll ans=pi[x]-(sx+s-2)*(sx-s+1)/2;
        for (register int i=s+1;i<=sx;i++) {
        	ans+=pi[x/b[i]];
		}
        return ans;
	}
	return getphi(x,s-1)-getphi(x/b[s],s-1);
}
ll getpi(ll x){
	if(x<=block){
		return pi[x];
	}
	ll ans=getphi(x,pi[ll(pow(x,1.0/3.0))])+pi[ll(pow(x,1.0/3.0))]-1;
	for(register ll i=pi[ll(pow(x,1.0/3.0))]+1,ed=pi[ll(pow(x,1.0/2.0))];i<=ed;i++){
		ans-=getpi(x/b[i])-i+1;
	}
	return ans;
}
inline __int128 read(){
    __int128 x(0),f(1);
    char ch=getchar();
    while(ch<'0'||ch>'9'){
        if(ch=='-')
            f=-1;
        ch=getchar();
    }
    while(ch>='0'&&ch<='9'){
        x=x*10+ch-'0';
        ch=getchar();
    }
    return x*f;
}
__int128 pre(__int128 x,long long s,__int128 k){
	__int128 ans=1;
	while(s){
		if(s&1){
			ans=(ans*x)%k;
		}
		x=(x*x)%k;
		s>>=1;
	}
	return ans;
}
ll S(ll ns,ll kk){
	if(ns<b[kk]){
		return 0;
	}
	ll ans=(4*getpi(ns)%MOD-kk*4%MOD+MOD)%MOD;
	for(ll j=kk+1;j<=r&&b[j]*b[j]<=ns;j++){
		for(ll kp=1,l=b[j];l<=ns;kp++,l*=b[j]){
			ans+=((S(ns/l,j)+(kp!=1))%MOD*(kp*3+1))%MOD;
			ans%=MOD;
		} 
	}
	return max(ans,(ll)0);
}
int main(){
	t=read();
	for(int i=1;i<=t;i++){
		n[i]=read();
		MX=max(MX,n[i]);
	}
	shai(1);
	sz[0]=1;
	for(int i=0;i<=30;i++){
		p[i][0]=i;
	}
	for(int i=1;i<=64;i++){
		MOD*=2;
	}
	for(int i=1;i<=2;i++){
		sz[i]=sz[i-1]*b[i];
		for(int j=1;j<=30;j++){
			p[j][i]=p[j][i-1]-p[j/b[i]][i-1];
		}
	}
	for(int i=1;i<=t;i++){
		if(n[i]<=100000){
			write(e[n[i]]);
			printf("\n");
			continue;
		}
		write((S(n[i],(ll)0)+(ll)1)%MOD);
		printf("\n");
	}
	return 0;
}

一个min25和meissel-lehmer的合体，不知道为什么TLE。