CF2081G2 Hard Formula (Hard Version) 题解

目录

验题人来报个到。

题目描述

给定 \(N\) ，求 \((\sum_{i=1}^Ni\bmod\varphi(i))\bmod 2^{32}\) 。

数据范围

• \(1\le N\le 10^{12}\) ，但可以加强到 \(10^{13}\) 。

时间限制 \(\texttt{10s}\) ，空间限制 \(\texttt{1024MB}\) 。

分析

一、\(N\le 10^8\)

线性筛即可做到 \(\mathcal O(N)\) 。

二、\(N\le 10^{10}\)

记 \(g(n)=\frac n{\varphi(n)}=\prod_{p|n}\frac p{p-1},g'(n)=\lfloor g(n)\rfloor\) ，则：

\[res=\frac{N(N+1)}2-\sum_{i=1}^Ng'(i)\varphi(i)\\ \]

问题转化为计算 \(\sum_{i=1}^Ng'(i)\varphi(i)\) 。

\(\texttt{Key observation}\) ：若 \(g'(np)\neq g'(n)\) ，则 \(p\le n\) ，唯二反例是 \(n=1,p=2\) 和 \(n=2,p=3\) 。

证明：

\[\because g(n)\cdot\frac p{p-1}\ge g'(n)+1\\ \therefore p\le\frac{g'(n)+1}{g'(n)+1-g(n)}\\ \therefore p\le\frac{\lfloor\frac n{\varphi(n)}\rfloor\varphi(n)+\varphi(n)}{\lfloor\frac n{\varphi(n)}\rfloor\varphi(n)+\varphi(n)-n}\\ \texttt{let}\ x=\lfloor\frac n{\varphi(n)}\rfloor\varphi(n)+\varphi(n)-n\ge 1\\ \therefore p\le\frac{x+n}x\le n+1\\ \]

若 \(p=n+1\) ，由取等条件 \(x=1\) 知 \(\varphi(n)\mid n+1\) 。

若 \(n\gt 2\) ，则 \(\varphi(n)\ge 2\) ，与 \(p\) 为素数矛盾； \(n\le 2\) 手动验证即可。

因此大于 \(\sqrt N\) 的质因子 \(p\) 不会改变 \(g'(n)\) 的值。

筛出 \(\sqrt N\) 以内的质因子，仿照 \(\min25\) 筛合数部分的套路，令：

\[S(n,j,u)=\sum_{i=2}^n[\text{minp}(i)\gt p_j]\lfloor g(i)\cdot u\rfloor\varphi(i)\\ \]

则：

\[S(n,j,u)=\sum_{id_p\gt j,p^e\le n}\varphi(p^e)(\lfloor u\cdot\frac p{p-1}\rfloor+S(\lfloor\frac n{p^e}\rfloor,id_p,u\cdot\frac p{p-1}))\\ \]

显然我们不可能枚举所有 \(p\) ，但是当 \(p\gt\sqrt n\) 且 \(u\cdot\frac p{p-1}\) 不再进位时内层贡献为 \(\varphi(p)\) ，可以通过 \(\min25\) 筛计算 \(\varphi\) 在素数位置的前缀和。

时间复杂度 \(\mathcal O(\frac{n^{0.75}}{\log n})\) ，实测跑 \(10^{11}\) 轻轻松松。

#include<bits/stdc++.h>
#define ui unsigned
#define int long long
using namespace std;
const int maxn=2e6+5;
int B,m,n,cnt;
int p[maxn>>3],phi[maxn];
bitset<maxn> b;
ui g1[maxn],g2[maxn],s1[maxn];
int id1[maxn],id2[maxn],val[maxn];
void init(int N)
{
    phi[1]=1;
    for(int i=2;i<=N;i++)
    {
        if(!b[i]) p[++cnt]=i,phi[i]=i-1;
        for(int j=1;j<=cnt&&i*p[j]<=N;j++)
        {
            b[i*p[j]]=1;
            if(i%p[j]==0) {phi[i*p[j]]=phi[i]*p[j];break;}
            phi[i*p[j]]=phi[i]*(p[j]-1);
        }
    }
    if(n<=N)
    {
        ui res=0;
        for(int i=1;i<=n;i++) res+=i%phi[i];
        printf("%u\n",res),exit(0);
    }
    for(int i=1;i<=cnt;i++) s1[i]=s1[i-1]+p[i];
}
int id(int x)
{
    return x<=B?id1[x]:id2[n/x];
}
ui dfs(int n,int j,double u)
{
    if(p[j]>=n) return 0;
    int i=j+1,_u=floor(u);
    ui res=0;
    if(j==cnt) res=_u*(g1[id(n)]-g2[id(n)]-(s1[j]-j));
    else
    {
        for(;;i++)
        {
            double v=u*p[i]/(p[i]-1);
            ui _v=floor(v),y=p[i]-1;
            for(int x=p[i],e=1;x<=n;x*=p[i],y*=p[i],e++)
                res+=y*(_v+dfs(n/x,i,v));
            if(i==cnt||(_v==_u&&p[i]*p[i]>n)) break;
        }
        if(n>p[i]) res+=_u*(g1[id(n)]-g2[id(n)]-(s1[i]-i));
    }
    return res;
}
signed main()
{
    ///freopen("data.in","r",stdin);
    ///freopen("data.out","w",stdout);
    auto st=clock();
    scanf("%lld",&n),init(max(B=sqrt(n),1000000ll));
    for(int l=1,r=0;l<=n;l=r+1)
    {
        int x=n/l;
        r=n/x,val[++m]=x,x<=B?id1[x]=m:id2[n/x]=m;
        g1[m]=x*(x+1ull)/2-1,g2[m]=x-1;
    }
    for(int j=1;j<=cnt;j++)
        for(int i=1;p[j]*p[j]<=val[i];i++)
        {
            int k=id(val[i]/p[j]);
            g1[i]-=p[j]*(g1[k]-s1[j-1]),g2[i]-=g2[k]-(j-1);
        }
    printf("%u\n",n*(n+1ull)/2-1-dfs(n,0,1));
    auto ed=clock();
    fprintf(stderr,"time = %.3lf ms\n",1.0*(ed-st)/CLOCKS_PER_SEC);
    return 0;
}

三、\(N\le 10^{12}\)

建出 \(\min25\) 筛的搜索树：共 \(N\) 个节点， \(1\) 号点为根， \(d\) 的父节点为 \(\frac d{\text{maxp}(d)}\) 。

在 \(g'\) 变化的位置统计贡献：

\[res=\sum_{i=1}^N\varphi(i)+\sum_{g'(d)\neq g'(fa_d)}G(d,p=\text{maxp}(d))\\ \]

其中 \(G(d,p)\) 为 \(d\) 子树中所有节点的 \(\varphi\) 之和。

注意到 \(d\) 是非叶节点 \(\iff dp\le N\) ，非叶节点数量很少，满足 \(g'(d)\neq g'(fa_d)\) 的非叶结点更少。

\(N\)	\(g'(d)\neq g'(fa_d)\) 的非叶节点数	搜索时间
\(10^{10}\)	78370	0.015s
\(10^{11}\)	285027	0.058s
\(10^{12}\)	1074403	0.182s
\(10^{13}\)	4010656	0.650s

搜索使用的代码片段如下：

void dfs(ll d,int j,double u)
{
    auto check=[&]()->bool
    {///剪枝,如果子树中不存在g'(u)\neq g'(fa_u)的点,直接返回
        double v=u;
        ll x=d;
        for(int i=j+1;i<=cnt;i++)
        {
            if((__int128)x*p[i]>n) break;
            x*=p[i];
            if(floor(v=v*p[i]/(p[i]-1))>floor(u)) return 1;
        }
        return 0;
    };
    if(!check()) return ;
    for(int i=max(j,1);i<=cnt&&(__int128)d*p[i]*p[i]<=n;i++)
    {
        double v=i!=j?u*p[i]/(p[i]-1):u;
        sum+=floor(v)!=floor(u),dfs(d*p[i],i,v);
    }
}

---

接下来对叶节点和非叶节点分别计算贡献。

对于前者，对每个（搜到的） \(d\) ，它符合要求的叶儿子 \(dq\) 是一段区间，贡献为：

\[\sum_d\sum_{q\gt p,dq^2\gt n\ge dq}[g'(dq)\neq g'(d)]\varphi(dq)\\ \]

将 \(\varphi(dq)\) 拆成 \(\varphi(d)\varphi(q)\) ，在 \(g'\) 序列中二分找到 \(q\) 的左右端点，预处理 \(\varphi\) 在素数位置的前缀和即可。

对于后者，我们需要对大约 \(10^6\) 个满足 \(g'(d)\neq g'(fa_d)\) 的 \(d\) 计算：

\[G(d,p=\text{maxp}(d))=\varphi(\frac dp)\sum_{i=1}^{\lfloor N/\frac dp\rfloor}[\text{minp}(i)=p]\varphi(i)\\ \]

令 \(f(n,p)=\sum_{i=1}^n[\text{minp}(i)\ge p]\varphi(i)\) ，这里 \(\text{minp}(1)=\infty\) ，初始值 \(f(n,0)\) 可以用杜教筛计算。

杜教筛需要轻微卡常：

• map+dfs 会导致复杂度多一只 \(\log\) ，改成递推求所有 \(S(\lfloor\frac ni\rfloor)\) 的值。
• 整数除法比浮点数除法慢很多，统一换成浮点数除法。
• 常规写法要枚举 \(\lfloor\frac ni\rfloor\) 相同的区间，改成枚举 \(i\) 和 \(\lfloor\frac ni\rfloor\) 。

具体参见代码。

则 \(G(d)=\varphi(\frac dp)(f(\lfloor N/\frac dp\rfloor,p)-f(\lfloor N/\frac dp\rfloor,p+1))\) 。

将第一维压入狄利克雷生成函数：

\[F_p(x)=\sum_{i=1}^\infty[\text{minp}(i)\ge p]\frac{\varphi(i)}{i^x}=\prod_{q\ge p}\sum_{k\ge 0}\frac{\varphi(q^k)}{q^{kx}} \]

则 \(f(n,p)=\sum_{i=1}^n[\frac 1{i^x}]F_p(x)\) 。

从 \(F_p(x)\) 向 \(F_{p+1}(x)\) 转移的过程相当于除以 \(\sum_{k\ge 0}\frac{\varphi(q^k)}{q^{kx}}\) 。由于 \(id=\varphi*1\) ，所以除以 \(\varphi\) 相当于除以 \(id\) 再乘上 \(1\) 。

除以 \(id\) ：

\[H_p(x)=\frac{F_p(x)}{\sum_{k\ge 0}\frac{p^k}{p^{kx}}}=F_p(x)(1-\frac p{p^x})\\ h(n,p)=f(n,p)-p\cdot f(\lfloor\frac np\rfloor,p)\\ \]

乘上 \(1\) ：

\[F_{p+1}(x)=H_p(x)\sum_{k\ge 0}\frac 1{p^{kx}}\\ f(n,p+1)=\sum_{k\ge 0}h(\lfloor\frac n{p^k}\rfloor,p)=h(n,p)+f(\lfloor\frac np\rfloor,p+1)\\ \]

单次递推代价为 \(\mathcal O(\sqrt N)\) ，但是总共有 \(\mathcal O(\frac{\sqrt N}{\log\sqrt N})\) 个素数 \(p\) ，所以直接递推时间复杂度仍然过高。

对 \(p\le N^\frac 14\) 的部分递推，这里只有 \(\mathcal O(\frac{N^\frac 14}{\log N^\frac 14})\) 个素数。

若 \(p=\text{maxp}(d)\le N^\frac 14\) ，在递推到 \(p\) 时计算 \(G(d)\) 。

若 \(p=\text{maxp}(d)\gt N^\frac 14\) ，由 \(g'(d)\neq g'(fa_d)\) 知 \(d\gt p^2\gt\sqrt N\) 。

换一种方式计算 \(G(d)\) ：

\[G(d)=\sum_{k\ge 0}\varphi(dp^k)f(\lfloor\frac N{dp^k}\rfloor,p+1)\\ \]

其中 \(\lfloor\frac N{dp^k}\rfloor\le\sqrt N\) ，因此我们只需考虑前 \(\sqrt N\) 个数。

对 \(1\le i\le\sqrt N\) ，若 \(\varphi(i)\) 对 \(f(x,y)\) 产生贡献，则 \(i\le x,\text{minp}(i)\ge y\) 。这是一个二维偏序。

时间复杂度 \(\mathcal O(\frac{N^\frac 34}{\log N}+\sqrt N\log N)\) 。

#include<bits/stdc++.h>
#define ui unsigned
#define int long long
#define div _div
using namespace std;
const int N1=1e8+5,N2=2e6+5,N3=3e5+5;
int B,k1,k2,m,n,cnt;
ui res;
ui p[N1>>4],phi[N1];
bitset<N1> b;
double g[N3];
ui c[N2],f[N2],h[N2],l[N2],s[N3];
int id1[N2],id2[N2],val[N2];
struct node
{
    ui u,v,w;
}e1[N2],e2[N2];
void init(int N)
{
    phi[1]=1;
    for(int i=2;i<=N;i++)
    {
        if(!b[i]) p[++cnt]=i,phi[i]=i-1;
        for(int j=1;j<=cnt&&i*p[j]<=N;j++)
        {
            b[i*p[j]]=1;
            if(i%p[j]==0) {phi[i*p[j]]=phi[i]*p[j];break;}
            phi[i*p[j]]=phi[i]*(p[j]-1);
        }
    }
    if(n<=N)
    {
        for(int i=1;i<=n;i++) res+=i%phi[i];
        printf("%u\n",res),exit(0);
    }
    for(B=sqrt(n);p[cnt]>B;cnt--) ;
    for(int i=1;i<=cnt;i++) g[i]=p[i]/(p[i]-1.0),s[i]=s[i-1]+phi[p[i]];
    for(int i=1;i<=N;i++) phi[i]+=phi[i-1];
}
int div(int x,int j)
{
    return (double)x/j;
}
ui id(int x)
{
    return x<=B?id1[x]:id2[n/x];
}
void get_phi(int n)
{
    static int tmp[N2];
    for(int l=1,r=0;l<=n;l=r+1)
    {
        int x=n/l;
        r=n/x,val[++m]=x,x<=B?id1[x]=m:id2[n/x]=m;
    }
    for(int i=m;i>=1;i--)
    {
        int x=val[i],t=sqrt(x);
        if(x<=N1-5) f[i]=phi[x];
        else
        {
            f[i]=x*(x+1ull)/2;
            /*ver 1:for(int l=2,r=0;l<=x;l=r+1) r=x/(x/l),f[i]-=(r-l+1)*f[id(x/l)];*/
            /*ver 2:for(int j=2;j<=t;j++) f[i]-=f[id(x/j)];
                    for(int j=x/t-1;j>=1;j--) f[i]-=(x/j-x/(j+1))*f[id(j)];*/
            /*ver 3:for(int j=2;j<=t;j++) f[i]-=f[id(div(x,j))];
                    for(int j=x/t-1;j>=1;j--) f[i]-=(div(x,j)-div(x,j+1))*f[id(j)];*/
            tmp[x/t]=div(x,x/t);
            for(int j=x/t-1;j>=1;j--) tmp[j]=div(x,j),f[i]-=(tmp[j]-tmp[j+1])*f[id(j)];
            for(int j=2;j<=t;j++) f[i]-=f[id(tmp[j])];
            ///ver 5&6 做法常数更小,在1e13部分会提到
            ///f[i]=x*(x+1ull)/2+t*phi[t];
            /*ver 5:for(int j=1;j<=t;j++) f[i]-=div(x,j)*(phi[j]-phi[j-1]);
                    for(int j=2;j<=t;j++) f[i]-=f[id(div(x,j))];*/
            /*ver 6:for(int j=1,u;j<=t;j++) f[i]-=(u=div(x,j))*(phi[j]-phi[j-1])+(j!=1?f[id(u)]:0);*/
        }
    }
}
void dfs(int d,int j,double u,ui _phi)
{
    auto check=[&]()->bool
    {
        double v=u;
        for(int i=j+1,x=d;i<=cnt;x*=p[i++])
        {
            if((__int128)x*p[i]>n) break;
            if(floor(v*=g[i])>floor(u)) return 1;
        }
        return 0;
    };
    if(!check()) return ;
    int i=max(j,1ll),_u=floor(u);
    for(;i<=cnt&&(__int128)d*p[i]*p[i]<=n;i++)
    {
        double v=i!=j?u*g[i]:u;
        if(floor(v)!=_u)
        {
            if(1ll*p[i]*p[i]<=B) e1[++k1]=(node){id(n/d),i,_phi};
            else for(int x=d*p[i],y=_phi*(p[i]-1);x<=n;x*=p[i],y*=p[i]) e2[++k2]=(node){n/x,i+1,y};
        }
        dfs(d*p[i],i,v,_phi*(p[i]-(i!=j)));
    }
    ///part 1:计算叶儿子的贡献
    i+=i==j;
    if(i<=cnt&&d*p[i]<=n&&floor(u*g[i])!=_u)
    {
        int l=i,r=cnt+1;
        while(r-l>1)
        {
            int mid=(l+r)>>1;
            d*p[mid]<=n&&floor(u*g[mid])!=_u?l=mid:r=mid;
        }
        res+=_phi*(s[l]-s[i-1]);
    }
}
void add(int x,int v)
{
    while(x) c[x]+=v,x-=x&-x;
}
ui query(int x)
{
    ui res=0;
    while(x<=cnt) res+=c[x],x+=x&-x;
    return res;
}
signed main()
{
    ///freopen("data.in","r",stdin);
    ///freopen("data.out","w",stdout);
    auto st=clock();
    scanf("%lld",&n);
    init(N1-5),get_phi(n),res=f[id(n)],dfs(1,0,1,1);
    ///part 2:计算 minp(d)<=n^0.25 的关键点贡献
    sort(e1+1,e1+k1+1,[&](node x,node y){return x.v<y.v;});
    for(ui j=1,l=1,r=0;j<=cnt&&p[j]*p[j]<=B;j++)
    {
        for(r=l;r<=k1&&e1[r].v==j;r++) ;
        for(ui i=l;i<r;i++) res+=f[e1[i].u]*e1[i].w;
        for(ui i=1;i<=m;i++) h[i]=f[i]-p[j]*f[id(div(val[i],p[j]))];
        for(ui i=m;i>=1;i--) f[i]=h[i]+f[id(div(val[i],p[j]))];
        for(ui i=l;i<r;i++) res-=f[e1[i].u]*e1[i].w;
        l=r;
    }
    ///part 3:计算 minp(d)>n^0.25 的关键点贡献
    sort(e2+1,e2+k2+1,[&](node x,node y){return x.u<y.u;});
    l[1]=cnt;
    for(ui j=1;j<=cnt;j++) for(ui i=p[j];i<=e2[k2].u;i+=p[j]) if(!l[i]) l[i]=j;
    for(ui i=1,j=1;i<=k2;i++)
    {
        while(j<=e2[i].u) add(l[j],phi[j]-phi[j-1]),j++;
        res+=query(e2[i].v)*e2[i].w;
    }
    printf("%u\n",n*(n+1ull)/2-res);
    auto ed=clock();
    fprintf(stderr,"time = %.3lf ms\n",1.0*(ed-st)/CLOCKS_PER_SEC);
    return 0;
}

四、\(N\le 10^{13}\)

考虑一个杜教筛的变形。记 \(t=\lfloor\sqrt n\rfloor\) ，由 \(xy\le n\) 知 \(\min(x,y)\le t\) 。枚举 \(x,y\) 中 \(\le t\) 的部分，有：

\[S_{f*g}(n)=\sum_{i=1}^tf(i)S_g(\lfloor\frac ni\rfloor)+\sum_{i=1}^tg(i)S_f(\lfloor\frac ni\rfloor)-S_f(t)S_g(t)\\ \]

代入 \(f(x)=\varphi(x),g(x)=1\) ，有：

\[S_{\varphi}(n)=\frac{n(n+1)}2+t\cdot S_{\varphi}(t)-\sum_{i=1}^t\lfloor\frac ni\rfloor\varphi(i)-\sum_{i=2}^tS_{\varphi}(\lfloor\frac ni\rfloor)\\ \]

从而减小杜教筛常数，这也是 \(10^{12}\) 部分杜教筛中 ver 4 和 ver 5 的内容。

---

我们需要对杜教筛进一步优化。

引理： \([1,n]\) 中 \(\text{minp}(i)\ge n^\frac 13\) 的数只有 \(\mathcal O(\frac n{\log n})\) 个。

证明：

显然 \(\text{minp}(i)\ge n^\frac 13\) 的数只有 \(p\) 和 \(pq\) 两种形式，前者个数为 \(\mathcal O(\frac n{\log n})\) 。

对于后者：

\[\begin{aligned} &\sum_{p\in\text{prime},p\ge n^\frac 13}\mathcal O(\frac{\frac np}{\log\frac np})\\ =&\mathcal O(\frac n{\log n})\sum_{p\in\text{prime},p\ge n^\frac 13}\mathcal O(\frac 1p)\\ =&\mathcal O(\frac n{\log n})\cdot\mathcal O(\log\log n-\log\log n^\frac 13)\\ =&\mathcal O(\frac n{\log n})\\ \end{aligned} \]

取 \(C=N^\frac 16\) ，考虑先求出 \(F(\lfloor\frac Ni\rfloor,C)\) ，通过 \(\texttt{Dgf}\) 逆推（ \(10^{12}\) 部分讲了如何正推，转移方程类似），每次乘上 \(\sum_{k\ge 0}\frac{\varphi(q^k)}{q^{kx}}\) 即可。

令 \(f(x)=[\text{minp}(i)\gt C]\varphi(x),g(x)=[\text{minp}(x)\gt C]\) ，则 \((f*g)(x)=[\text{minp}(x)\gt C]x\) 。

代入杜教筛变形的式子：

\[S_f(n)=S_{f*g}(n)+S_f(t)\cdot S_g(t)-\sum_{i=1}^tS_g(\lfloor\frac ni\rfloor)f(i)-\sum_{i=2}^tS_f(\lfloor\frac ni\rfloor)g(i)\\ \]

目标求出 \(S_f(n)=F(n,C)\) 。

其中 \(S_{f*g}(n)\) 和 \(S_g(n)\) 可以用 \(\min25\) 筛素数部分的 \(\texttt{dp}\) 完成，代价为 \(\mathcal O(\frac C{\log C}\cdot\sqrt N)\) 。

根据引理，求和号枚举代价为 \(\mathcal O(\frac t{\log t})\) 。

预处理 \([1,B]\) 中 \(S_g(n)\) 的值，剩余部分用上式递推，代价为：

\[\begin{aligned} &\mathcal O(B+\sum_{i=1}^{\frac NB}\frac{\sqrt \frac Ni}{\log\frac Ni})\\ =&\mathcal O(B+\frac{\sqrt N}{\log N}\sqrt x\mid_1^{\frac NB})\\ =&\mathcal O(B+\frac N{\log N}\cdot\frac1{\sqrt B})\\ =&\mathcal O(\frac{N^\frac 23}{\log^\frac 23N})\\ \end{aligned} \]

至此，我们可以在 \(\mathcal O(\frac{N^\frac 23}{\log^\frac 23N})\) 的时间内完成杜教筛。

---

如果你真的写了 \(10^{12}\) 的部分分，就会发现 part 2 和 part 3 的时间分配极不均衡（\(\texttt{dp}\) 跑得很慢，树状数组跑得飞快），这启示我们重新设置分治阈值。

对 \(\text{minp}(i)\le C\) 的部分 \(\texttt{dp}\) ，时间复杂度 \(\mathcal O(\frac C{\log C}\cdot\sqrt N)\) 。

对 \(\text{minp}(i)\gt C\) 的部分跑树状数组，注意只需考虑 \([1,\frac N{C^2}]\) 中 \(\text{minp}\gt C\) 的数，时间复杂度 \(\mathcal O(\frac{\frac N{C^2}}{\log\frac N{C^2}}\log N)\) 。

取 \(C=N^\frac 16\log^\frac 13N\approx 515\) ，可以得到理论最优复杂度 \(\mathcal O(\frac{N^\frac 23}{\log^\frac 23N})\) ，为平衡常数取 \(C=300\) 最优。

代码实现时，可以先求出 \(F(\lfloor\frac Ni\rfloor,C)\) 的值，一边逆推一边做 part 2 的 \(\texttt{dp}\) ，最后加上 \(\varphi(N)\) 的贡献。

时间复杂度 \(\mathcal O(\frac{N^\frac 23}{\log^\frac 23N})\) 。

#include<bits/stdc++.h>
#define ui unsigned
#define int long long
#define div _div
using namespace std;
const int k=62,N1=5e7+5,N2=6.4e6+5,N3=3e5+5;///对[1,C=p[k]]跑dp,其余跑树状数组
int B,m,n,k1,k2,cnt;
ui res;
ui p[N1>>4],phi[N1],l[N1],s[N3];
bitset<N1> b;
double g[N3];
ui c[N2],f[N2],h[N2];
int id1[N2],id2[N2],val[N2];
vector<ui> vec;
struct node
{
    ui u,v,w;
}e1[N2],e2[N2];
void init(int N)
{
    phi[1]=1;
    for(int i=2;i<=N;i++)
    {
        if(!b[i]) p[++cnt]=i,phi[i]=i-1,l[i]=cnt;
        for(int j=1;j<=cnt&&i*p[j]<=N;j++)
        {
            b[i*p[j]]=1,l[i*p[j]]=j;
            if(i%p[j]==0) {phi[i*p[j]]=phi[i]*p[j];break;}
            phi[i*p[j]]=phi[i]*(p[j]-1);
        }
    }
    if(n<=N)
    {
        for(int i=1;i<=n;i++) res+=i%phi[i];
        printf("%u\n",res),exit(0);
    }
    for(B=sqrt(n);p[cnt]>B;cnt--) ;
    for(int i=1;i<=max(cnt,k);i++) g[i]=p[i]/(p[i]-1.0),s[i]=s[i-1]+p[i];
    vec.reserve(N2);
    for(int i=1;i<=N;i++) if(i==1||l[i]>k) vec.push_back(i);
    vec.push_back(N+1);///防止溢出
}
int div(int x,int j)
{
    return (double)x/j;
}
int id(int x)
{
    return x<=B?id1[x]:id2[n/x];
}
void get_phi()
{///求F(N/i,C)的值
    static ui g[N2],h[N2];
    for(int l=1,r=0;l<=n;l=r+1)
    {
        int x=n/l;
        r=n/x,val[++m]=x,x<=B?id1[x]=m:id2[n/x]=m;
        g[m]=x-1,h[m]=x*(x+1ull)/2-1;
    }
    for(int j=1;j<=k;j++)
        for(int i=1;p[j]*p[j]<=val[i];i++)
        {
            int u=id(div(val[i],p[j]));
            g[i]-=g[u]-(j-1),h[i]-=p[j]*(h[u]-s[j-1]);
        }
    for(int i=1;i<=m;i++) g[i]=(val[i]>p[k]?1+g[i]-k:1),h[i]=(val[i]>p[k]?1+h[i]-s[k]:1);///补上1的贡献,删去素数的贡献
    for(int i=m,j=0,cur=0;i>=1;i--)
    {
        int x=val[i],t=sqrt(x);
        if(x<=N1-5)
        {
            while(vec[j]<=x) cur+=phi[vec[j++]];
            f[i]=cur;
        }
        else
        {
            f[i]=h[i]+f[id(t)]*g[id(t)];
            for(int j=0,u=0;vec[j]<=t;j++) f[i]-=g[u=id(div(x,vec[j]))]*phi[vec[j]]+(j?f[u]:0);
        }
    }
    ///加上下面一段即可得到 \varphi 的块筛
    ///for(int j=k;j>=1;j--)
    ///{
        /*乘以id:h(n,p)=\sum_{k\ge 0}f(\llfoor\frac n{p^k}\rfloor,p)=p\cdot h(\lfloor\frac np\rfloor,p)+f(n,p)
          除以1: F_{p-1}(x)=H_p(x)\cdot(1-\frac 1{p^x}),f(n,p-1)=h(n,p)-h(\lfloor\frac np\rfloor,p)*/
        /*ver 1:for(int i=m;i>=1;i--) h[i]=f[i]+p[j]*h[id(div(val[i],p[j]))];
                for(int i=1;i<=m;i++) f[i]=h[i]-h[id(div(val[i],p[j]))];*/
        ///for(int i=m,u=0;i>=1;i--) h[i]=f[i]+p[j]*h[u=id(div(val[i],p[j]))],f[i]=h[i]-h[u];
    ///}
    ///cerr<<f[id(n)]<<endl;
}
void dfs(int d,int j,double u,ui _phi)
{
    auto check=[&]()->bool
    {
        double v=u;
        for(int i=j+1,x=d;i<=cnt;x*=p[i++])
        {
            if((__int128)x*p[i]>n) break;
            if(floor(v*=g[i])>floor(u)) return 1;
        }
        return 0;
    };
    if(!check()) return ;
    int i=max(j,1ll),_u=floor(u);
    for(;i<=cnt&&(__int128)d*p[i]*p[i]<=n;i++)
    {
        double v=i!=j?u*g[i]:u;
        if(floor(v)!=_u)
        {
            if(i<=k) e1[++k1]=(node){id(n/d),i,_phi};
            else for(int x=d*p[i],y=_phi*(p[i]-1);x<=n;x*=p[i],y*=p[i]) e2[++k2]=(node){n/x,i+1,y};
        }
        dfs(d*p[i],i,v,_phi*(p[i]-(i!=j)));
    }
    ///part 1:计算叶儿子的贡献
    i+=i==j;
    if(i<=cnt&&d*p[i]<=n&&floor(u*g[i])!=_u)
    {
        int l=i,r=cnt+1;
        while(r-l>1)
        {
            int mid=(l+r)>>1;
            d*p[mid]<=n&&floor(u*g[mid])!=_u?l=mid:r=mid;
        }
        res+=_phi*(s[l]-s[i-1]-l+i-1);
    }
}
void add(int x,int v)
{
    if(x>k) while(x) c[x]+=v,x-=x&-x;///仅考虑minp>C的数,但bit不是常数瓶颈,所以优化有限
}
ui query(int x)
{
    ui res=0;
    while(x<=cnt) res+=c[x],x+=x&-x;
    return res;
}
signed main()
{
    ///freopen("data.in","r",stdin);
    ///freopen("data.out","w",stdout);
    auto st=clock();
    scanf("%lld",&n),init(N1-5),get_phi(),dfs(1,0,1,1);
    ///part 2:计算 minp(d)<=C 的关键点贡献
    sort(e1+1,e1+k1+1,[&](node x,node y){return x.v>y.v;});
    for(ui j=k,l=1,r=0;j>=1;j--)
    {
        for(r=l;r<=k1&&e1[r].v==j;r++) ;
        for(ui i=l;i<r;i++) res-=f[e1[i].u]*e1[i].w;
        for(ui i=m,u=0;i>=1;i--) h[i]=f[i]+p[j]*h[u=id(div(val[i],p[j]))],f[i]=h[i]-h[u];
        for(ui i=l;i<r;i++) res+=f[e1[i].u]*e1[i].w;
        l=r;
    }
    ///part 3:计算 minp(d)>C 的关键点贡献
    sort(e2+1,e2+k2+1,[&](node x,node y){return x.u<y.u;});
    l[1]=cnt=3001134;///[1,N1-5]素数个数
    for(ui i=1,j=1;i<=k2;i++)
    {
        while(j<=e2[i].u) add(l[j],phi[j]),j++;
        res+=query(e2[i].v)*e2[i].w;
    }
    res+=f[id(n)],printf("%u\n",n*(n+1ull)/2-res);
    auto ed=clock();
    fprintf(stderr,"time = %.3lf ms\n",1.0*(ed-st)/CLOCKS_PER_SEC);
    return 0;
}

花絮

CF 审核员说不想开太多档部分分所以删掉了 \(10^{13}\) 这一档。

验题人最初只发现了 \(n=2,p=3\) 的反例，但没发现 \(n=1,p=2\) 的反例，以至于出题人最初版本的暴力能过 \(N=10^{10}\) 但过不了 \(N\le 10\) 的数据（什么离谱的 corner case），直到开赛前一周才发现。

CF 上 c++17 语言是 32 位机，并且用不了 __int128 （\(10^{12}\) 不需要，但 \(10^{13}\) 会用到）。

浮点数除法优化效果非常明显，将 \(\mathcal O(\frac{N^{0.75}}{\log N})\) 次整数除法换成浮点数除法，可以从 \(\texttt{TLE}\) 变成 \(\texttt{6s-}\) 。

出题人的 \(10^{10}\) 暴力拿来跑 \(10^{12}\) 的数据范围比我 \(10^{12}\) 做法还快，所以开时限的时候特别尴尬，最后决定为了不卡常而放过了部分常数极小的暴力做法。

比赛因为一些离谱的技术原因被推迟并 unr 了（哭）。主站被 DDos 然后宣布推迟，结果镜像站没推迟并且部分选手看到并保存了题面（原贴已删除）。