hdu6706 huntian oy

Posted 2022-08-22 encodetalker

hdu
好久没写数论函数题了，上一次写还是在纪中学min25筛的时候了，赶紧来一道补下手感

题面：求
\[ \sum_i=1^n\sum_j=1^igcd(i^a-j^a,i^b-j^b)[gcd(i,j)=1] \]
保证\(n,a,b\leq 10^9,gcd(a,b)=1\)

知道\((i^a-j^a)\)这个玩意因式分解后有\((i-j)\)，不妨大力猜一下\(gcd(i^a-j^a,i^b-j^b)=i-j\)，证明明天再补
那么\[ \beginaligned 原式=&\sum_i=1^n\sum_j=1^i(i-j)[gcd(i,j)=1]\=&\sum_i=1^n\sum_j=1^i(i-j)\sum_d|gcd(i,j)\mu(d)\=&\sum_d=1^n\mu(d)d\sum_i=1^\lfloor\fracnd\rfloor\sum_j=1^ii-j\\endaligned \]

后面那个你拿平方和公式算一下应该有\(\sum_i=1^n\sum_j=1^ii-j=\fracn^3-n6\),前面那一部分就是\(\mu·id\)，卷上\(id\)之后就是\(\epsilon\),杜教筛直接求

感觉我回忆杜教筛的时间都比前面推导的时间长。。。

#include<iostream>
#include<string.h>
#include<string>
#include<stdio.h>
#include<algorithm>
#include<vector>
#include<math.h>
#include<queue>
#include<set>
#include<map>
#include<unordered_map>
using namespace std;
typedef long long ll;
typedef long double db;
typedef pair<int,int> pii;
const int N=6000000;
const db pi=acos(-1.0);
#define lowbit(x) (x)&(-x)
#define sqr(x) (x)*(x)
#define rep(i,a,b) for (register int i=a;i<=b;i++)
#define per(i,a,b) for (register int i=a;i>=b;i--)
#define fir first
#define sec second
#define mp(a,b) make_pair(a,b)
#define pb(a) push_back(a)
#define maxd 1000000007
#define inv6 166666668
#define eps 1e-8
int pri[N+10],tot=0,mu[N+10];
int sum[N+10];
bool vis[N+10];
unordered_map<int,ll> has;

int read()

    int x=0,f=1;char ch=getchar();
    while ((ch<'0') || (ch>'9')) if (ch=='-') f=-1;ch=getchar();
    while ((ch>='0') && (ch<='9')) x=x*10+(ch-'0');ch=getchar();
    return x*f;


void sieve()

    mu[1]=1;
    rep(i,2,N)
    
        if (!vis[i]) pri[++tot]=i;mu[i]=-1;
        for (int j=1;j<=tot && i*pri[j]<=N;j++)
        
            vis[i*pri[j]]=1;
            if (i%pri[j]==0) break;
            else mu[i*pri[j]]-=mu[i];
        
    
    rep(i,1,N) sum[i]=(sum[i-1]+i*mu[i]+maxd)%maxd;


ll calc(int n)

    ll ans=1ll*n*n%maxd*n%maxd;
    ans=(ans-n+maxd)%maxd;
    ans=ans*inv6%maxd;
    return ans;


ll query(int n)

    if (n<=N) return sum[n];
    if (has[n]) return has[n];
    ll ans=1;int l,r;
    for (l=2;l<=n;l=r+1)
    
        r=n/(n/l);
        ll tmp=(1ll*(r+l)*(r-l+1)/2)%maxd;
        ans=(ans-tmp*query(n/l)%maxd+maxd)%maxd;
    
    has[n]=ans;
    return ans;

        
ll f(int n,int a,int b)

    ll ans=0;
    int l,r;
    for (l=1;l<=n;l=r+1)
    
        r=n/(n/l);
        ll nowsum=(query(r)-query(l-1)+maxd)%maxd;
        ans=(ans+nowsum*calc(n/l)%maxd)%maxd; 
    
    return ans;


int main()

    int T=read();
    sieve();
    while (T--)
    
        int n=read(),a=read(),b=read();
        printf("%lld\n",f(n,a,b));
    
    return 0;