Description

給定n,m，求\(\sum_{i=1}^{n}\sum_{j=1}^{m}\varphi(ij)\)模10^9+7的值。

Input

僅一行，兩個整數n,m。

Output

僅一行答案。

Sample Input

100000 1000000000

Sample Output

857275582

數據規模

1<=n<=10^5，1<=m<=10^9。

sol

%%%ranwen！！！

前置技能：

1. \(n=\sum_{d|n}\varphi(d)\)
2. \(\varphi(ij)=\varphi(\frac{i}{d})*\varphi(j)*d\quad[d=(i,j)]\)
3. \([\mu(x)=1||-1]\quad\varphi(\frac{x}{d})*\varphi(\frac{d}{e})=\varphi(\frac{x}{e})\quad[e|d]\)

證明：1不多說，2的意思就是讓ij互質然後直接分解成兩個phi，接著把gcd產生的倍數貢獻乘回去，3因為x沒有平方因子。除以d之後就不會有d 的因子了，所以與\(\frac{d}{e}\)互質，滿足積性函數性質，乘起來即可。

解法：

首先觀察數據範圍可知，n的範圍較小，可以進行枚舉，而m的範圍極大，已經超過了線性篩的範圍。

我們考慮枚舉i，然後推一波式子：

設\(s(n,m)=\sum_{i=1}^{m}\varphi(ni)\)

設w是n所有質因子一次方的乘積，v=n/w，則：

\(s(n,m)=v*\sum_{i=1}^{m}\varphi(iw)\)

設\(d=(i,w)\)，然後用公式2得：

\(s(n,m)=v*\sum_{i=1}^{m}\varphi(i)*\varphi(\frac{w}{d})*d\)

用公式1，得：

\(s(n,m)=v*\sum_{i=1}^{m}\varphi(i)*\varphi(\frac{w}{d})*\sum_{e|d}\varphi(\frac{d}{e})\)

用公式3，得：

\(s(n,m)=v*\sum_{i=1}^{m}\varphi(i)*\sum_{e|i,e|w}\varphi(\frac{w}{e})\)

因為 \(d=(i,w)\)，所以：

\(s(n,m)=v*\sum_{e|w}\varphi(\frac{w}{e})*\sum_{i=1}^{\lfloor\frac{m}{e}\rfloor}\varphi(ie)\)

根據\(s(n,m)\)的定義得：

\(s(n,m)=v*\sum_{e|w}\varphi(\frac{w}{e})*s(e,\frac{m}{e})\)

當n等於1的時候，可以直接使用杜教篩計算。

直接記憶化搜索即可。

因為每一步的m都是\(\lfloor\frac{x}{y}\rfloor\)的形式，所以可以使用下底分塊法來計算。

時間復雜度\(O(n\sqrt{m}+n^{\frac{2}{3}})\)。

代碼

#include <bits/stdc++.h>
using namespace std;
int n,m,tot,ans,phi[1000005],sum[1000005],pri[1000005],low[1000005],vis[1000005],P=1e9+7;
map<pair<int,int>,int>a;map<int,int>b;
int djs(int x)
{
    if(x<=1e6) return sum[x];
    if(b.count(x)) return b[x];
    int tp=1ll*x*(x+1)/2%P,last;
    for(int i=2;i<=x;i=last+1) last=x/(x/i),tp=(tp-1ll*(last-i+1)*djs(x/i)%P+P)%P;
    b.insert(make_pair(x,tp));return tp;
}
int solve(int x,int y)
{
    if(!x||!y) return 0;
    if(x==1) return djs(y);
    if(y==1) return phi[x];
    if(a.count(make_pair(x,y))) return a[make_pair(x,y)];
    int w=low[x],v=x/w,lim=floor(sqrt(w)+0.1),tp=0;
    for(int i=1;i<=lim;i++) if(w%i==0)
    {
        tp=(tp+1ll*phi[w/i]*solve(i,y/i)%P)%P;
        if(i!=w/i) i=w/i,tp=(tp+1ll*phi[w/i]*solve(i,y/i)%P)%P,i=w/i;
    }
    tp=1ll*tp*v%P;a.insert(make_pair(make_pair(x,y),tp));
    return tp;
}
int main()
{
    phi[1]=low[1]=sum[1]=1;
    for(int i=2;i<=1000000;sum[i]=(sum[i-1]+phi[i])%P,i++)
    {
        if(!vis[i]){pri[++tot]=i;phi[i]=i-1;low[i]=i;}
        for(int j=1;j<=tot&&i*pri[j]<=1000000;j++)
        {
            vis[i*pri[j]]=1;
            if(i%pri[j]==0){phi[i*pri[j]]=phi[i]*pri[j];low[i*pri[j]]=low[i];break;}
            phi[i*pri[j]]=phi[i]*(pri[j]-1),low[i*pri[j]]=low[i]*pri[j];
        }
    }
    scanf("%d%d",&n,&m);
    if(n>m) swap(n,m);
    for(int i=1;i<=n;i++) ans=(ans+solve(i,m))%P;
    printf("%d\n",ans);
}

【bzoj3512】DZY Loves Math IV 杜教篩+記憶化搜索+歐拉函數