huntian oy

问题类别	杜教筛	算法类别	数论
其他知识点	无	来源	HDU-6706

Problem Description
One day, Master oy created a new function to celebrate his becoming a 'huntian' in majsoul.

\(f(n,a,b)=\sum_{i=1}^n \sum_{j=1}^i gcd(i^a−j^a,i^b−j^b)[gcd(i,j)=1]\%(10^9+7)\)

Given n, a and b, Master oy wanted Newbie jj who was still a 'chuxin' to answer the value of f(n,a,b).

Input
There are multiple test cases.
The first line contains an integer T, indicating the number of test cases.
For each test case, there are three positive integers n, a and b which are separated by spaces. It's guaranteed that a and b are coprime.

\(1≤n,a,b≤10^9\)

\(T=10^4\), but there are only 10 test cases that n is over \(10^6\).

Output
For each test case, an integer in one line representing your answer.

Sample Input
2
1 2 3
100 2 3

Sample Output
0
101542

题意：求解\(f(n,a,b)=\sum_{i=1}^n \sum_{j=1}^i gcd(i^a−j^a,i^b−j^b)[gcd(i,j)=1]\)

思路：打表可得 \(gcd(i^a−j^a,i^b−j^b)=i-j\)
原式化简为\[\sum_{i=1}^n \sum_{j=1}^i (i-j)[gcd(i,j)=1]\]
然后分成两部分：

\[\sum_{i=1}^n \sum_{j=1}^i i[gcd(i,j)=1]-\sum_{i=1}^n \sum_{j=1}^i j[gcd(i,j)=1]\]

左边部分显然就是：\(\sum_{i=1}^n i\varphi(i)\)

主要简化右边部分：\(\sum_{i=1}^n \sum_{j=1}^i j[gcd(i,j)=1]\)

推导过程：

\[\sum_{i=1}^n\sum_{j=1}^i j\sum_{d|gcd(i,j)}u(d)\]

\[\sum_{i=1}^n\sum_{d|i}u(d)\sum_{j=1}^{i}j[d|j]\]

\[\sum_{i=1}^n\sum_{d|i}u(d)\sum_{j=1}^{\lfloor \frac id \rfloor}jd\]

\[\sum_{i=1}^n\sum_{d|i}u(d)d \frac {(1+{\lfloor \frac id \rfloor}){\lfloor \frac id \rfloor}}{2}\]

\[\sum_{i=1}^n\frac i2 \sum_{d|i}u(d)(1+\lfloor \frac id\rfloor)\]

\[\sum_{i=1}^n(\frac i2 \sum_{d|i}u(d)+\frac i2\varphi(i))\]

\[\sum_{i=1}^n \frac{i\ \varphi(i)}{2}+\sum_{i=1}^n\frac {i[i=1]}{2}\]

题目条件告诉我们n是大于等于1的，因此右边部分就是 \(\frac 12\)

所以总的式子合起来就是：

\[\sum_{i=1}^n i\varphi(i)-\sum_{i=1}^n \frac{i\ \varphi(i)}{2}-\frac 12\]
最终所求即：
\[\sum_{i=1}^n \frac{i\ \varphi(i)}{2}-\frac 12\]

求积性函数前缀和我们可以使用杜教筛：\(h=f*g\)
我们设 \(f(d)=d\varphi(d)\)，显然我们卷上id函数就可以凑成h
\[h=\sum_{d|n}d\varphi(d) \frac nd=\sum_{d|n}\varphi(d)n=n^2\]

利用杜教筛的公式（公式其实不用记，手推一遍也就1分钟的事情 ）：
\[g(1)S(n)=\sum_{i=1}^n h(i)-\sum_{d=2}^n g(d)S(\lfloor \frac nd \rfloor)\]
并根据推得的g、h函数化简
\[S(n)=\sum_{i=1}^n i^2-\sum_{d=2}^n d\times S(\lfloor \frac nd \rfloor)\]

这样就做完了

#include <iostream>
#include <algorithm>
#include <cstdio>
#include <queue>
#include <vector>
#include <set>
#include <map>
#include <unordered_map>
#include <cstring>
#include <string>
#include <cmath>
#define rep(i,a,b) for (int i=a; i<=b; ++i)
#define per(i,b,a) for (int i=b; i>=a; --i)
#define mes(a,b)  memset(a,b,sizeof(a))
#define mp make_pair
#define ll long long
#define pb push_back
#define pii pair<int,int>
#define pll pair<ll,ll>
#define ls (rt<<1)
#define rs ((rt<<1)|1)
#define isZero(d)  (abs(d) < 1e-8)
using namespace std;
const int mod=1e9+7;

int t;
int n,a,b;
ll inv2,inv6;

int qpow(int base,int n,int mod)
{
    int ret=1;
    while(n)
    {
        if(n&1)
            ret=1ll*ret*base%mod;
        base=1ll*base*base%mod;
        n>>=1;
    }
    return ret;
}

const int N=2e6;
int prime[N+10],phi[N+10],visit[N+10];
int cnt;

void init()
{
    cnt=0,phi[1]=1;
    for(int i=2;i<=N;++i)
    {
        if(!visit[i])
            prime[++cnt]=i,phi[i]=i-1;
        for(int j=1;j<=cnt&&i*prime[j]<=N;++j)
        {
            visit[i*prime[j]]=1;
            if(i%prime[j]==0)
            {
                phi[i*prime[j]]=phi[i]*prime[j];
                break;
            }
            phi[i*prime[j]]=phi[i]*phi[prime[j]];
        }
    } 
    for(int i=1;i<=N;++i)
    {
        phi[i]=1ll*phi[i]*i%mod;
        phi[i]=(phi[i]+phi[i-1])%mod;
    }
}

unordered_map<int,int> mphi;

int pre_x3(int n)
{
    return 1ll*n*(n+1)%mod*(2*n+1)%mod*inv6%mod;
}

int pre_phi(int n)
{
    if(n<=N)
        return phi[n];
    if(mphi[n])
        return mphi[n];
    
    int ans=pre_x3(n);
    int l=2,r;
    while(l<=n)
    {
        r=n/(n/l);
        ans=(1ll*ans-1ll*pre_phi(n/l)*(r-l+1)%mod*(l+r)%mod*inv2%mod+mod)%mod;
        l=r+1;
    }
    return mphi[n]=ans;
}

int abss(int x)
{
    if(x<0)
        x+=mod;
    return x;
}

int main()
{
    inv2=qpow(2,mod-2,mod);
    inv6=qpow(6,mod-2,mod);
    init();
    scanf("%d",&t);
    while(t--)
    {
        scanf("%d%d%d",&n,&a,&b);
        int ans=abss(1ll*pre_phi(n)*inv2%mod-inv2)%mod;
        printf("%d\n",ans);        
    }
    return 0;
}