hdu 5780 gcd（线性筛+快速幂+数论）

题目描述

传送门

题目大意：
∑1<=a,b<=ngcd(xa−1,xb−1)

题解

首先对于gcd(xa−1,xb−1)进行化简
更相减损gcd(a,b)=gcd(a−b,b)
那么gcd(xa−1,xb−1)=gcd(xa−xb,xb−1)=gcd(xb(xa−b−1),xb−1)
又因为xb,xb−1 互质，所以
gcd(xa−1,xb−1)=gcd(xa−b−1,xb−1)
gcd(xa−1,xb−1)=xgcd(a,b)−1

那么原始的式子就变成了∑1<=a,b<=nxgcd(a,b)−1

∑k=1n(xk−1)∑i=1n∑j=1n[gcd(i,j)=k]

∑k=1n(xk−1)∑i=1⌊nk⌋∑j=1⌊nk⌋[gcd(i,j)=1]

因为[n=1]=∑d|nμ(d),所以

∑k=1n(xk−1)∑i=1⌊nk⌋∑j=1⌊nk⌋∑d|gcd(i,j)μ(d)

∑k=1n(xk−1)∑d=1⌊nk⌋∑i=1⌊nk⌋[d|i]∑j=1⌊nk⌋[d|j]μ(d)

∑k=1n(xk−1)∑d=1⌊nk⌋⌊⌊nk⌋d⌋2μ(d)

这样可以O(T∗n34) 利用分块求解，但是如何继续优化呢？
考虑对于F[m]=∑mi=1∑mj=1[gcd(i,j)=1] 进行更直观的理解

F[m]=2∗(∑a=1m∑b=1a[gcd(a,b)=1])−1

F[m]=2∗(∑a=1mϕ(a))−1

对于(xk−1)这一部分，我们对于每一块讲将所有的-1提出来，然后利用等比数列求和公式计算即可。

代码

#include<iostream>#include<cstdio>#include<algorithm>#include<cstring>#include<cmath>#define LL long long #define N 1000003#define p 1000000007using namespace std;int n,T,prime[N],pd[N],mu[N],SUM[N];LL a[N],x,ans,smu[N],phi[N];LL gcd(LL x,LL y){    LL r;    while (y) {        r=x%y;        x=y; y=r;    }     return x;}void init(){    mu[1]=1; phi[1]=1;    for (int i=2;i<=1000000;i++) {        if (!pd[i]) {            prime[++prime[0]]=i;            mu[i]=-1; phi[i]=i-1;        }        for (int j=1;j<=prime[0];j++) {            if (i*prime[j]>1000000) break;            pd[i*prime[j]]=1;            if (i%prime[j]==0) {                mu[i*prime[j]]=0;                phi[i*prime[j]]=phi[i]*prime[j];                break;            }            else {                mu[i*prime[j]]=-mu[i];                phi[i*prime[j]]=phi[i]*(prime[j]-1);            }        }    }    for (int i=1;i<=1000000;i++) phi[i]=(phi[i]+phi[i-1])%p;}LL quickpow(LL num,int x){    LL ans=1; LL base=num%p;    while (x) {        if (x&1) ans=ans*base%p;        x>>=1;        base=base*base%p;    }    return ans;}int main(){    scanf("%d",&T);    init();    while (T--) {        scanf("%I64d%d",&x,&n);        ans=0; int j=0;        for (int i=1;i<=n;i=j+1) {            j=min(n/(n/i),n);            LL t=2*phi[n/i]-1; t%=p;              LL a1=quickpow(x,i); LL q=quickpow(x,(j-i+1));            LL t1=a1*(q-1)%p*quickpow(x-1,p-2)%p-(j-i+1);            //cout<<t<<" "<<t1<<endl;            ans=(ans+t*t1%p)%p;        }        if (x==1) printf("0\n");        else  printf("%I64d\n",(ans%p+p)%p);    }}