【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