HDU 5780 gcd （欧拉函数）

要做这题你首先要知道一个不常见的公式：点击打开链接

gcd(a^m-b^m,a^n-b^n)=a^(gcd(m,n))-b^(gcd(m,n))  (a>b)

即：

gcd ({x}^{a}-1x​a​​−1,{x}^{b}-1x​b​​−1)={x}^{gcd(a,b)}-1x​gcd(a,b)​​−1成立，相当于计算\sum \sum {x}^{gcd(a,b)}-1∑∑x​gcd(a,b)​​−1 

由此O(n^2)平方可以算出答案。但显然n^2大大的超时了。

比赛的时候我想的是要找一个O(N)的解法，但实际上O(N)也不行，因为题目给的T是300，n是1e6，O(N)也是1e8超时！

先想O(N),再优化。

记s[d]=最大公约数为d的对数，答案\sum s[d]*({x}^{d}-1)∑s[d]∗(x​d​​−1). 

我们的目的就是找到[1,n]的的每个d对应的s[d]，就可以O(n)的做出来了。

题解说：

先用筛法计算出欧拉函数。把正方形分成上三角和下三角计算，最大公约数为d的数量s[d]=2*(phi[1]+phi[2]+...+phi[n/d])-1。

竟然和欧拉函数扯上了关系，我极其吃力的理解为其下：

我们考虑把数这样排列：

（1，1）

（1，2）

（2，2）

（1，3）

（2，3）

（3，3）

（1，4）

（2，4）

（3，4）

（4，4）

（....）

（n，n）

也就是说一共有n中结尾，每种结尾有i组数。

d=1时，对于每一种结尾，显然phi[i]即为这种结尾gcd(i,n)==1的个数。所以s[1]=（phi[1]+phi[2]+...+phi[n])

d=2时，只有结尾为偶数的情况下，两个数的gcd才有可能等于2；所以只需要考虑结尾为2，4，6，8...n(或n-1）的情况，一共是n/2种；那么每种结尾里面到底又有多少组数的gcd正好是2呢。

观察规律发现 正好是 1，1，2，2，4....正好s[2]=（phi[1]+phi[2]+...phi[n/2]；其实就相当于把d/=2后，在前n/2种结尾里找gcd=1的数

d=3时，只有结尾为3的倍数的情况下，两个数的gcd才可能等于3，所有只需要考虑3，6，9....n（或n-1，或n-2）的情况，一共是n/种，相当于把d/=3后，在前n/3种结尾里找gcd=1的数。

写到这里，也总算明白了，s[d]=2*(phi[1]+phi[2]+...+phi[n/d])-1大概是怎么来的了，这里 *2 -1 ，应该不用说什么吧？

这样，O(n)的写法也算是写出来了。只要先预处理一下欧拉函数，并且算一个前缀和，等要算的时候直接用就是了。

但是前面也说了，O（n）还是会超时，这里题解又说了：

观察s[d]计算公式，可以发现对不同的d，若n/d相同，s[d]不发生变化。根据s[d]分段计算，相同的一段的{x}^{d}-1x​d​​−1可以用等比公式求。复杂度为（n+T\sqrt{n}lognn+T√​n​​​logn)

也就是说，会出现大片n/d相等的情况，比如10/4=10/5=2;10/6=10/7=10/8=10/9=10/10=1；所以在求和的时候，没必要一个个的求，而可以一段一段的求。

令 j=n/(n/d)；j即为n/d这个值最边缘的位置了，比如n=10，d=6，则j=10，即找到了10这个10/d=1最后的位置，循环中d=j+1，就可以把中间的数跳过了，大大优化的时间！

而跳过的这些，它们的s[d]都是一样的，∑s[d]∗(x​d​​−1)，只要搞一下等比数列求和，把后半部分的和算出来就行了。

最后，就是代码了！这题基本上是对着题解和ranklist第一名的代码强行理解的，不得不承认，这种难度的题目对我来说太吃力了。尽管我已经知道了那个gcd公式，比比人有优势，但是还是做不出来

$【代码】$

/* ***********************************************Author        :angon************************************************ */#include <stdio.h>#include <string.h>#include <iostream>#include <algorithm>#include <stack>#include <vector>#include <queue>#include <set>#include <map>#include <string>#include <math.h>#include <stdlib.h>#include <time.h>using namespace std;#define REP(i,k,n) for(int i=k;i<n;i++)#define REPP(i,k,n) for(int i=k;i<=n;i++)#define scan(d) scanf("%d",&d)#define scann(n,m) scanf("%d%d",&n,&m)#define mst(a,k)  memset(a,k,sizeof(a));#define LL long long#define N 1000005#define mod 1000000007/*inline int read(){    int s=0;    char ch=getchar();    for(; ch<'0'||ch>'9'; ch=getchar());    for(; ch>='0'&&ch<='9'; ch=getchar())s=s*10+ch-'0';    return s;}inline void print(int x){    if(!x)return;    print(x/10);    putchar(x%10+'0');}*/LL power(LL x,LL n)    //x^n log(n)；{    LL res=1;    while(n>0)    {        if(n & 1)            res=(res*x)%mod;        x=(x*x)%mod;        n >>= 1;    }    return res;}LL inv[N];int prime[N],sphi[N];void init(){    sphi[1]=1;    for(int i=2;i<=N;++i){        if(!sphi[i]){            prime[++prime[0]]=i;            sphi[i]=i-1;        }        for(int j=1;j<=prime[0]&&i*prime[j]<=N;++j)            if(i%prime[j])sphi[i*prime[j]]=sphi[i]*(prime[j]-1);            else sphi[i*prime[j]]=sphi[i]*prime[j];    }    for(int i=1;i<=N;++i) sphi[i]=(sphi[i]+sphi[i-1])%mod;    inv[1] = inv[0] = 1;    for(int i = 2;i < N;i++)    inv[i] = inv[mod%i]*(mod-mod/i)%mod;}LL sn(LL q,LL n){    if(q==1) return n;    return (power(q,n)-1)*inv[q-1]%mod;}int main(){    init();    int t,x,n;    scan(t);    while(t--)    {        scann(x,n);        LL ans=0;        for(int i=1,j;i<=n;i=j+1)        {            j = n/(n/i);            LL sd = 2*sphi[n/i]-1;             ans = (ans + sd * (power(x,i) * sn(x,j-i+1)%mod -(j-i+1)) % mod) % mod;        }        printf("%I64d\n",ans);    }    return 0;}