HDU6607Easy Math Problem（min_25筛+杜教筛+拉格朗日插值）

HDU6607Easy Math Problem（min_25筛+杜教筛+拉格朗日插值）

题目大意

求
$$\sum _{i=1}^n\sum _{j=1}^{n}gcd(i,j{)}^{k}lcm(i,j)[gcd(i,j)\in prime]$$

其中$1\le n\le 10^{10},1\le k\le 100$

解题思路

$$\begin{array}{rl}\sum _{i=1}^n\sum _{j=1}^{n}gcd(i,j{)}^{k}lcm(i,j)[gcd(i,j)\in prime]& =\sum _p\sum _{i=1}^{\lfloor \frac{n}{p}\rfloor }\sum _{j=1}^{\lfloor \frac{n}{p}\rfloor }{p}^{k+1}\cdot (i\cdot j)[gcd(i,j)=1]\\ & =\sum _p{p}^{k+1}\sum _{i=1}^{\lfloor \frac{n}{p}\rfloor }\sum _{j=1}^{\lfloor \frac{n}{p}\rfloor }\cdot (i\cdot j)[gcd(i,j)=1]\end{array}$$

对于${\sum }_{i=1}^{\lfloor \frac{n}{p}\rfloor }{\sum }_{j=1}^{\lfloor \frac{n}{p}\rfloor }\cdot (i\cdot j)[gcd(i,j)=1]$单独分析
$$\sum _{i=1}^n\sum _{j=1}^{n}j\cdot i[gcd(i,j)=1]$$

可以对原本的区间进行拆分
$$[1\le i\le n][1\le j\le n]=[1\le i\le n][1\le j<i]+[1\le i<j][1\le j\le n]+[i=j][1\le i\le n]$$
对第一部分求和
$$\sum _{i=1}^{n}i\sum _{j=1}^{i-1}j[gcd(i,j)]$$
根据公式KaTeX parse error: \tag works only in display equations就有
$$\sum _{i=1}^n\frac{i[i=1]+i^2\phi (i)}{2}$$
而区间的第二部分和第一部分实质上是一样的，那么前面两部分就是${\sum }_{i=1}^{n}i[i=1]+i^2\phi (i)$而第三部分除了$i=j=1$的情况以外，其总为0，而当其为1时答案恰好为1，故最终的答案就是
$$\sum _{i=1}^n\sum _{j=1}^{n}j\cdot i[gcd(i,j)=1]=\sum _{i=1}^n{i}^2\phi (i)$$
而对于$i{d}^2\cdot \phi$这个函数可以通过与$i{d}^2$卷积
$$(i{d}^2\cdot \phi )\ast (i{d}^2)=\sum _{d|n}{d}^2\phi (d)\frac{n^2}{d^2}=n^2\sum _{d|n}\phi (d)=n^3$$
其中$n^3$和$n^2$都有公式可以$O(1)$求出，因此可以很方便地应用杜教筛求解

而对于${\sum }_p{p}^{k+1}$可以通过min25筛中的子方案求出所有的整除分块中需要用到的前缀和

这个问题的主要复杂度就是通过拉格朗日插值处理出前$\sqrt{n}$的${\sum }_{i=1}^n{i}^{k+1}$的前缀和，这部分的复杂度是$O(k\sqrt{n})$

加上min25筛加上杜教筛的复杂度，总的复杂度就是$O(k\sqrt{n}+\frac{n^{\frac{3}{4}}}{logn}+n^{\frac{2}{3}})$

AC代码

#include<bits/stdc++.h>
using namespace std;
#define int long long
const int size=2e5+5;
const int mod=1e9+7;
typedef long long LL;
bool prime[size]; int p[size];
int tot=0;
int inv6,inv2;
int phi[size];
int sumiiphi[size];
int quick_pow(int a,int b)
{
    register int ans=1;
    while(b)
    {
        if(b&1) ans=1LL*ans*a%mod;
        a=1LL*a*a%mod;
        b>>=1;
    }
    return ans;
}
int inv[205];
void init()
{
    tot=0;
    phi[1]=1;
    for(int i=1;i<size;i++) prime[i]=true;
    for(int i=2;i<size;i++)
    {
        if(prime[i])
            p[++tot]=i,phi[i]=i-1;
        for(int j=1;j<=tot&&p[j]*i<size;j++) 
        {
            prime[i*p[j]]=false;
            if(i%p[j]==0) 
            {
                phi[i*p[j]]=phi[i]*p[j];
                break;
            }
            else
            phi[i*p[j]]=phi[i]*phi[p[j]];                
        }
    }
    sumiiphi[0]=0;
    for(int i=1;i<size;i++) sumiiphi[i]=(sumiiphi[i-1]+1LL*i*i%mod*phi[i]%mod)%mod;
    inv6=quick_pow(6,mod-2);
    inv2=quick_pow(2,mod-2);
    inv[0]=inv[1]=1;
    for(int i=2;i<=200;i++)
	inv[i] = (mod - (mod / i)) * inv[mod % i] % mod;
	for (int i = 1; i <= 200; ++i)inv[i] = inv[i] * inv[i - 1] % mod;
}
int tol;
int hk[size<<1];LL num[size<<1];
int pre[size];
int coeff[205];
void pre_lage(int k)
{
    for(int i=1;i<=k+2;i++) coeff[i]=quick_pow(i,k);
    coeff[0]=0;
    for(int i=1;i<=k+2;i++) coeff[i]=(coeff[i-1]+coeff[i])%mod;
}
LL suf[205],bef[205];
int lage(int n,int k)
{
    if(n<=k+2) return coeff[n];
    bef[0] = suf[0] = 1;
	for (int i = 1; i <= k + 2; ++i) {
		bef[i] = bef[i - 1] * ((n - i) % mod) % mod;
		suf[i] = suf[i - 1] * ((n + i - k - 3) % mod) % mod;
	}
	LL res = 0;
	for (int i = 1; i <= k + 2; ++i) {
		LL s = coeff[i] * bef[i - 1] % mod * suf[k - i + 2] % mod * inv[i - 1] % mod * inv[k + 2 - i] % mod;
		if ((k + 2 - i) & 1)s = -s;
		res += s;
		res = (res % mod + mod) % mod;
	}
	return res;
}
int s;LL n;
int id1[size],id2[size];
inline int ID(LL x){return x<=s?x:tol-n/x+1;}
void get_h(int k)
{
    s=sqrt(n);
    while(1LL*s*s<=n) s++;
    while(1LL*s*s>n) s--;
    pre[0]=0;
    for(register int i=1;p[i]<=s;i++)
    {
        pre[i]=(1LL*pre[i-1]+quick_pow(p[i],k+1))%mod;
    }
    tol=0;
    for(int i=1;i<=s;i++) num[++tol]=i;
    for(int i=s;i>=1;i--) if(n/i>s) num[++tol]=n/i;
    pre_lage(k+1);
    for(register int i=1;i<=tol;i++) hk[i]=(lage(num[i]%mod,k+1)-1+mod)%mod;
    hk[0]=0;
    int x=1;
    for(int j=1;j<=tot&&p[j]<=s;j++)
    {
        while(num[x]<p[j]*p[j]) x++;
        for(int i=tol;i>=x;i--)
        {
            hk[i]=((hk[i]-1LL*(pre[j]-pre[j-1]+mod)%mod*(hk[ID(num[i]/p[j])]-pre[j-1]))%mod+mod)%mod;
        }
    }
}
unordered_map<LL,int> mp;
inline int sum2(int n)
{
    return 1LL*n*(n+1)%mod*(2*n+1)%mod*inv6%mod;
}
inline int sum3(int n)
{
    int ans=1LL*n*(n+1)%mod*inv2%mod;
    return 1LL*ans*ans%mod;
}
inline int S(LL n)
{
    if(n<size) return sumiiphi[n];
    if(mp.count(n)) return mp[n];
    LL ans=sum3(n%mod);//n取模 
    for(LL l=2,r;l<=n;l=r+1)
    {
        r=n/(n/l);
        ans=((ans-1LL*S(n/r)*(sum2(r%mod)-sum2((l-1)%mod)))%mod+mod)%mod;
    }
    return mp[n]=ans;
}
int32_t main()
{
    int t;
    scanf("%lld",&t);
    int k;
    init();
    while(t--)
    {
        scanf("%lld%lld",&n,&k);
        mp.clear();
        get_h(k);
        int ans=0;
        for(int i=1;i<=tol;i++)
        {
            ans=(ans+1LL*(hk[i]-hk[i-1])%mod*S(n/num[i])%mod+mod)%mod;
        }
        printf("%lld\n",ans);
    }
}