SP34112题解
曾经尝试min26筛爆草，今天直接写了个PN发现本地开O2都跑了60s
显然有 \(\sigma^*(p^k)=p^k+1\)。
数据范围显然是不允许正常的筛法通过的（也许 zzq 的 \(O(\frac{n^{\frac{2}{3}}}{\log^2n})\) 能飞过去），于是考虑和构造有关的杜教筛和 PN 筛。
杜教筛由于其需要递归所以直接弃疗，考虑 PN 筛怎么做。
PN 筛最后的卷积是 \(O(\sqrt n)\) 的，只需要快速计算 \(g\) 即可。我们会快速计算的函数似乎只有 \(id^k\) 和 \(\sigma^k\)。
显然前者在这里是行不通的，所以考虑令 \(g=\sigma\)。
那么有：
\[\sum_{i=0}^{k}h(p^i)g(p^{k-i})=p^k+1
\]
计算 \(f(p^k)-f(p^{k-1})\)：
\[h(p^k)+\sum_{i=0}^{k-1}h(p^i)(g(p^{k-i})-g(p^{k-i-1}))=f(p^k)-f(p^{k-1})
\]
\[h(p^k)+\sum_{i=0}^{k-1}h(p^i)(\frac{p^{k-i+1}-p^{k-i}}{p-1})=f(p^k)-f(p^{k-1})
\]
\[h(p^k)+\sum_{i=0}^{k-1}h(p^i)p^{k-i}=f(p^k)-f(p^{k-1})
\]
再来一个差分，设有 \(F(p^k)=h(p^k)+\sum_{i=0}^{k-1}h(p^i)p^{k-i}=f(p^k)-f(p^{k-1})\)，计算 \(F(p^k)-p\times F(p^{k-1})\)：
\[h(p^k)+h(p^{k-1})\times p-p\times h(p^{k-1})=f(p^k)-f(p^{k-1})-p\times f(p^{k-1})+p\times f(p^{k-2})
\]
\[h(p^k)=f(p^k)-f(p^{k-1})-p\times f(p^{k-1})+p\times f(p^{k-2})
\]
对于特殊情况有：
\[h(1)=1
\]
\[h(p)=0
\]
\[h(p^2)=-p
\]
别的情况都是 \(h(p^k)=0\) 啦。
对于 \(g\) 设一个范围线性筛，剩下的范围 \(O(\sqrt n)\) 计算。容易发现 \(h\) 只在平方数处可能有值，也就是说复杂度是：
\[\sum_{B\leq\lfloor\frac{n}{i^2}\rfloor}\sqrt{\frac{n}{i^2}}
\]
\[\sqrt{n}\int_0^{\sqrt{\lfloor\frac{n}{B}\rfloor}}x^{-1}{\rm d}x
\]
复杂度为：
\[O(\sqrt{n}\ln\sqrt{n}-\sqrt{n}\ln\sqrt{B}+B)
\]
这里头顶上有个 \(\sqrt{n}\ln\sqrt{n}\) 特别大，于是只能看情况调整 \(B\)。。。
不过其实令 \(B=\sqrt{n}\)，就能够省下很大一部分常数啦。
而且其实瓶颈是线性筛。。。
upd:草，怎么spojrk3了
#include<cstdio>
#include<cmath>
typedef unsigned long long ull;
typedef __uint128_t LL;
typedef long long ll;
const int M=7072005,mod=1e9+7;
struct Barrett{
ull B,m;
Barrett(const ull&m=2):m(m),B((LL(1)<<63)/m){}
friend inline ull operator/(const ull&a,const Barrett&mod){
ull r=LL(mod.B)*a>>63;return(r+1)*mod.m<=a?r+1:r;
}
inline bool check(const ull&a){
return a==m*(a/ *this);
}
}B[M];
ull m,n[50005],N[M];int T,Block,top,pri[M],tmp[M],sigma[M];ull Sum[M];
inline ull Solve(const ull&n){
if(n<=m)return Sum[n];ull x(1),y,ans(0);while(x*x<=n)y=n/B[x],ans+=(y&1?y*(y+1>>1):(y>>1)*(y+1)),ans+=x*y,++x;
--x;return ans-(x&1?x*x*(x+1>>1):(x>>1)*x*(x+1));
}
inline void sieve(const int&n){
sigma[1]=1;for(int i=1;i<=n;++i)B[i]=Barrett(i);
for(int i=2;i<=n;++i){
if(!sigma[i])pri[++top]=i,tmp[i]=1,sigma[i]=i+1;
for(int x,j=1;j<=top&&(x=i*pri[j])<=n;++j){
if(B[pri[j]].check(i)){
tmp[x]=tmp[i];sigma[x]=sigma[i]*pri[j]+tmp[i];break;
}
sigma[x]=sigma[i]*sigma[pri[j]];tmp[x]=sigma[i];
}
}
for(int i=1;i<=n;++i)Sum[i]=Sum[i-1]+sigma[i];
}
inline ull DFS(const ull&n,const int&k){
const ull&N=n/B[pri[k]]/B[pri[k]];return !N||k==top+1?0:DFS(n,k+1)-pri[k]*(Solve(N)+DFS(N,k+1));
}
signed main(){
scanf("%d",&T);for(int i=1;i<=T;++i)scanf("%lld",n+i),n[i]>m&&(m=n[i]);sieve(m=sqrt(m));
for(int i=1;i<=T;++i)printf("%llu\n",DFS(n[i],1)+Solve(n[i]));
}
