【bzoj4804】欧拉心算 解题报告

【bzoj4804】欧拉心算

Description

给出一个数字(N),计算

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

Input

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

Output

按读入顺序输出答案


很多方法

可以推式子到

[sum_{T=1}^nlfloorfrac{n}{T} floor^2sum_{k|T}varphi(k)mu(frac{T}{k}) ]

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

也可以得到

[sum_{d=1}^n(varphi(d)(2sum_{i=1}^{lfloorfrac{n}{d} floor}varphi(i)-1)) ]

然后搞就行了

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

[sum_{i=1}^ksum_{j=1}^i[gcd(i,j)=1]=sum_{i=1}^kvarphi(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

原文地址:https://www.cnblogs.com/butterflydew/p/10142253.html