【bzoj4804】欧拉心算 解题报告

【bzoj4804】欧拉心算

Description

给出一个数字\(N\)，计算

\[\sum_{i=1}^n\sum_{j=1}^n \varphi(\gcd(i,j)) \]

Input

第一行为一个正整数\(T\)，表示数据组数。
接下来\(T\)行为询问，每行包含一个正整数\(N\)。
\(T\le 5000,N\le 10^7\)

Output

按读入顺序输出答案

---

很多方法

可以推式子到

\[\sum_{T=1}^n\lfloor\frac{n}{T}\rfloor^2\sum_{k|T}\varphi(k)\mu(\frac{T}{k}) \]

然后把后面的筛出来就行了

也可以得到

\[\sum_{d=1}^n(\varphi(d)(2\sum_{i=1}^{\lfloor\frac{n}{d}\rfloor}\varphi(i)-1)) \]

然后搞就行了

后面的式子化解用到了定义

\[\sum_{i=1}^k\sum_{j=1}^i[\gcd(i,j)=1]=\sum_{i=1}^k\varphi(i) \]

---

#include <cstdio>
#include <cctype>
#define ll long long
const int BufferSize=1<<16;
namespace Fast{
    const int LEN=10000000;
    char inp[LEN],outp[LEN];
    int tmp[20];
    int inpos,outpos;
    void init(){
        fread(inp,1,LEN,stdin);
        inpos=0; outpos=0;
    }
    char GetChar(){return inp[inpos++];}
    int read(){
        int ret=0; char ch=GetChar();
        while (ch<'0'||ch>'9') ch=GetChar();
        while ('0'<=ch&&ch<='9') ret=ret*10+ch-'0',ch=GetChar();
        return ret;
    }
    void PutChar(char ch){outp[outpos++]=ch;}
    void print(ll x){
        int pos=0;
        if (!x) tmp[++pos]=0;
        else
            while (x) tmp[++pos]=x%10,x/=10;
        for (int i=pos;i>=1;--i) PutChar(tmp[i]+'0');
    }
    void Print(){fwrite(outp,1,outpos,stdout);}
}
#define ll long long
const int N=1e7+1;
int pri[N],ispri[N],a[N],b[N],cnt;
ll f[N];
void init()
{
    b[1]=f[1]=1;
    for(register int i=2;i<N;i++)
    {
        if(!ispri[i])
        {
            pri[++cnt]=i;
            f[i]=i-2;
            a[i]=1;
            b[i]=i;
        }
        for(register int d,j=1,x;j<=cnt&&i*pri[j]<N;j++)
        {
            x=i*pri[j];
            ispri[x]=1;
            if(i%pri[j])
            {
                a[x]=1;
                b[x]=pri[j];
                f[x]=f[i]*f[pri[j]];
            }
            else
            {
                a[x]=a[i]+1;
                b[x]=b[i]*pri[j];
                d=i/b[i];
                if(d==1)
                    f[x]=b[i]/pri[j]*(pri[j]-1)*(pri[j]-1);
                else
                    f[x]=f[d]*f[b[i]*pri[j]];
                break;
            }
        }
    }
    for(int i=2;i<N;i++) f[i]+=f[i-1];
}
int main()
{
    Fast::init();
    init();
    int T,n;T=Fast::read();
    while(T--)
    {
        n=Fast::read();
        ll ans=0;
        for(register int d,l=1,r;l<=n;l=r+1)
        {
            d=n/l,r=n/d;
            ans+=1ll*d*d*(f[r]-f[l-1]);
        }
        Fast::print(ans),Fast::PutChar(' ');
    }
    Fast::Print();
    return 0;
}

---

2018.12.19