HDU 5780 gcd

Description

Little White learned the greatest common divisor, so she plan to solve a problem: given x,n,query $\sum_{i=1}^n\sum_{j=1}^n gcd(x^i−1,x^j−1)$

Solution

你妹这场BC跪成狗。。
T3没判1挂了。。
这道T5也不难，怎么会有点恐惧？？

首先的，可以$x^i-1$可以因式分解为$(x-1)\sum_{j=0}^{i-1}x^j$
令$f(x)=\sum_{i=0}^n-1 x^i$
看后面那个式子，显然对于某$(i,j)$,$gcd(x^i-1,x^j-1)=(x-1)*\sum_{k=0}^{gcd(i,j)}x^k$
再把外面的$(x-1)$乘进来，得到$gcd(x^i-1,x^j-1)=x^{gcd(i,j)}-1$
化简：
原式

$$=\sum_{i=1}^n \sum_{j=1}^n x^{gcd(i,j)}-1$$
$$=\sum_d \sum_{i=1}^{\lfloor \frac n d \rfloor} \sum_{j=1}^{\lfloor \frac n d \rfloor} (x^d-1)*[gcd(i,j)=1]$$
然后我就不知道自己在干啥了，抛了个莫比乌斯反演进去，得到了一个 $O(nlnn)$的做法，被卡掉了。。
事实上因为两个 $\sum$上限相同，可以直接用欧拉函数(当然莫比乌斯函数抛进去也可以通过卷积化掉变成欧拉函数)。。
继续：
$$=\sum_{d} (x^d-1)* (\sum_{i=1}^{\lfloor \frac n d \rfloor }2*\varphi[i] -1)$$
然后前缀和处理后面的 $\sum$，分段处理即可
当然一段 $x^d-1$可以用等比数列公式求和，所以还要快速幂和逆元

果然还是没见过世面，要推大半个小时才可以推出来

Code

#include<iostream>
#include<stdio.h>
#include<string.h>
#include<time.h>
#include<stdlib.h>
#include<math.h>
#include<algorithm>
#include<vector>
#include<string>
#include<queue>
#include<stack>
#include<set>
#include<map>
using namespace std;
typedef double db;
typedef unsigned ud;
typedef long long ll;
typedef unsigned long long ull;
typedef vector<int> vec;
typedef pair<int,int> pii;
typedef vector<pii> vecp;
#define pb push_back
#define ph push
#define fi first
#define se second
const int INF=(ud)-1>>1;
const ll inf=(ull)-1ll>>1;
template<class T>void rd(T &a){
    a=0;char c;
    while(c=getchar(),!isdigit(c));
    do a=a*10+(c^48);
        while(c=getchar(),isdigit(c));
}
template<class T>void nt(T x){
    if(!x)return;
    nt(x/10);
    putchar(48+x%10);
}
template<class T>void pt(T x){
    if(!x)putchar('0');
    else nt(x);
}
template<class T>void Max(T &a,T b){
    if(a<b)a=b;
}
template<class T>void Min(T &a,T b){
    if(a==-1||a>b)a=b;
}
const int M=1e6+5;
int prime[M],t;
bool mark[M];
inline bool chk(int x){
    int cnt=0;
    for(int i=0;prime[i]*prime[i]<=x;++i){
        if(x%prime[i]==0){
            x/=prime[i];
            if(x%prime[i]==0)return 0;
        }
    }
    return 1;
}
const int P=1e9+7;
int pw[M],phi[M],s[M];
inline void Mod_add(int &a,int b){
    if((a+=b)>=P)a-=P;
}
inline int Mod(int a){
    if(a>=P)return a-P;
    if(a<0)return a+P;
    return a;
}
inline int Mod_pow(int x,int a){
    int res=1;
    for(int i=0;1ll<<i<=a;++i){
        if(a&1<<i)res=1ll*res*x%P;
        x=1ll*x*x%P;
    }
    return res;
}
int x,n;
inline int calc(int a,int b){
    int w=Mod(Mod_pow(x,b-a+1)-1);
    int s=Mod_pow(x-1,P-2);
    int k=Mod_pow(x,a);
    return (1ll*k*w%P*s-(b-a+1))%P;
}
inline void gao(){
    cin>>x>>n;
    if(x==1){puts("0");return;}
    int res=0;
    for(int d=1;d<=n;++d){
        int last=(n/(n/d));
        Mod_add(res,1ll*calc(d,last)*s[n/d]%P);
        d=last;
    }
    cout<<res<<endl;
}
inline void Euler(){
    phi[1]=1;
    for(int i=2;i<M;++i){
        if(!mark[i])prime[t++]=i,phi[i]=i-1;
        for(int j=0,p;j<t&&(p=prime[j])*i<M;++j){
            mark[p*i]=1;
            if(i%p==0){phi[p*i]=p*phi[i];break;}
            phi[p*i]=phi[i]*(p-1);
        }
    }
    s[0]=-1;
    for(int i=1;i<M;++i)
        s[i]=Mod(s[i-1]+(phi[i]<<1));
}
int main(){
    Euler();
    int _;
    for(cin>>_;_--;)gao();
    return 0;
}