GDOI2016 之 GDSOI 第一题 互补约数

互补约数

题目大意

已知函数

求它的前缀和函数

输入格式

输入包含一行，一个正整数 n。

输出格式

输出只有一行， F(n)。

样例输入

Sample Input 1
10

Sample Input 2
1000

输出格式

Sample Output 1
32

Sample Output 2
12776

数据范围

20分算法

虽然侥幸进了决赛，但还是只能想到20分的。
不难看出，题目让我们求的是

∑ni=1∑⌊ni⌋j=1gcd(i,j)。

这样子枚举一下i，j，再打个gcd，求一下和就可以了。
这样20分就弄到手了。

100分算法

100分算法要用到莫比乌斯反演。
想学的可以自已查一下。
这里简单介绍它的内容。
mu(i)表示莫比乌斯函数。
如果满足

Fn=∑d|nGd

则满足

Gn=∑d|nFd∗mu(nd)

还有变式。
如果满足

Fi=∑⌊ni⌋d=1Gi∗d

则有

Gi=∑⌊ni⌋d=1Fi∗d∗mu(d)

如果想学的可以看我写的博客，链接如下:
http://blog.csdn.net/xianhaoming/article/details/51519100

好，介绍完了，现在进入正题。
我们设fk表示gcd(i,j)=k的对数，gk表示k|gcd(i,j)的对数。其中一定满足题目gcd(i,j)≤n√。
则现在满足

gk=∑⌊n√k⌋i=1fk∗i

根据莫比乌斯反演，得

fk=∑⌊n√k⌋i=1gk∗i∗mu(i)

那这样我们要求的答案就变成了

Ans=∑n√k=1k*fk

=∑n√k=1 k*∑⌊n√k⌋l=1gk∗l∗mu(l)

我们设

s=l*k

则

Ans=∑n√s=1 gs*(∑d|s mu(d)∗sd)

我们发现∑d|s mu(d)∗sd这一段只跟s有关，所以我们把此式子换成Ws，则

Ans=∑n√s=1 gs*Ws

我们先说一下W数组怎么预处理。
我们发现很多mu(i)的值为0，对W没有贡献。
于是我们可以这么处理：

    E:=trunc(sqrt(N));    for i:=1 to E do    if mu[i]<>0 then    for l:=1 to E div i do    M[i*l]:=M[i*l]+mu[i]*l;

现在我们看一下gs怎么求。
易得

gs=∑ns2i=1⌊ni∗s2⌋

现在我们尝试优化这个式子的时间复杂度。
我们发现对于多对（i，j），存在

⌊ni∗s2⌋=⌊nj∗s2⌋

这是我们可以分块，把一段⌊ni∗s2⌋的值相等的一段一次求完，这样就可以大大降低时间复杂度。详细做法见程序。

Code (pascal)

var    bz:array[0..400000] of boolean;    mu,jh,ss:array[0..400000] of int64;    j,k,l,i:longint;    cqy:extended;    n,m,lj,P,O,bj,le,r,ans,pp:int64;function min(a,b:INT64):INT64;    begin        if a<b then exit(a)        else exit(b);    end;begin    assign(input,'gcd.in'); reset(input);    assign(output,'gcd.out'); rewrite(output);    readln(n);    m:=trunc(sqrt(n));    mu[1]:=1;    for i:=2 to m+1 do    begin        if bz[i]=false then        begin            inc(o);            ss[o]:=i;            mu[i]:=-1;        end;        for l:=1 to o do        begin            if ss[l]*i>m then break;            bz[ss[l]*i]:=true;            mu[ss[l]*i]:=-mu[i];            if i mod ss[l]=0 then            begin                mu[ss[l]*i]:=0;                break;            end;        end;    end;    for i:=1 to m+1 do    if mu[i]<>0 then    for l:=1 to m div i do    jh[i*l]:=jh[i*l]+mu[i]*l;    for i:=1 to m do    begin        le:=1;        lj:=0;        cqy:=n/i/i;        bj:=n div i+1;        while le<=bj do        begin            pp:=trunc(cqy/le);            if pp<=0 then break;            r:=min(bj,trunc(cqy/pp));            lj:=lj+(r-le+1)*pp;            le:=r+1;        end;        ans:=ans+lj*jh[i];    end;    writeln(ans);    close(input);    close(output);end.