hdu 5780 gcd

首先是一个公式：

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

由此可得，本题要求x^gcd(a,b)-1;

再考虑每一个gcd值(d)可能对应多个数对，而只要求出对应多少个数对(sd)，就能直接由乘法得出结论。

那么怎么求N以内最大公约数是d的对数，那就是欧拉函数的事了。对每一个d，其值就是:

2*(phi[1]+...+phi[n/d])-1。这个不做过多解释。（减一是因为phi[1]在计算的时候会算上两遍，也就是扩大d倍之后a==b算了一遍，而b==a时又算了一遍）；

但是就是维护了前缀和，一个一个点算仍然超时，因为T太大。

所以考虑把相同sd值的数全都合并，跳着加一段，剩下的等比数列逆元搞。（比如对n=10来说，i=6一直到i=10都相等，因为n/i相等，这样每次延伸到最右端，也就是logn，而预处理的nlogn，是可以接受的）。

#include <cstdio>
#include <cstring>
#include <vector>
#include <cmath>
#include <stack>
#include <cstdlib>
#include <queue>
#include <map>
#include <iostream>
#include <algorithm>
#define mod 1000000007
using namespace std;
typedef long long LL;
int N = 1000100;


LL prime[1000100];
LL phi[1000100];
LL inv[1000100];
LL sum_phi[1000100];
bool is_prime[1000100];
LL quick_mod(LL a,LL b)
{
    LL ans=1;
    a%=mod;
    while(b)
    {
        if(b&1)
        {
            ans=ans*a%mod;
            b--;
        }
        b>>=1;
        a=a*a%mod;
    }
    return ans;
}
void init()
{
    phi[1]=1;
    for(int i=2; i<=N; ++i)
    {
        if(!phi[i])
        {
            prime[++prime[0]]=i;
            phi[i]=i-1;
        }
        for(int j=1; j<=prime[0]&&i*prime[j]<=N; ++j)
            if(i%prime[j])phi[i*prime[j]]=phi[i]*(prime[j]-1);
            else phi[i*prime[j]]=phi[i]*prime[j];
    }

    for (int i=1; i<N; i++)
    {
        phi[i]=(phi[i-1]+phi[i])%mod;
    }
    inv[1] = inv[0] = 1;
    for(int i = 2; i < N; i++)
        inv[i] = inv[mod%i]*(mod-mod/i)%mod;
//    for (int i=0;i<N;i++)
//    {
//        printf ("%lld ",phi[i]);
//        system("pause");
//    }
}
LL he(LL q,LL n)
{
    if (q==1) return n;
    return (quick_mod(q,n)-1)*inv[q-1]%mod;
}
LL n,x;
int main()
{
    init();
    int T;
    scanf ("%d",&T);
    while (T--)
    {
        scanf ("%I64d%I64d",&x,&n);
        int j;
        long long ans=0;
        for (int i=1; i<=n; i=j+1)
        {
            j=n/(n/i);
            LL sd=2*phi[n/i]-1;
            ans=(ans+sd*(quick_mod(x,i)*he(x,j-i+1)%mod-(j-i+1))%mod)%mod;
        }
        printf ("%I64d\n",ans);
    }
    return 0;
}