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

Solution

神題。。

考慮\(n\)比較小，試著枚舉\(n\)，處理\(m\)。

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

那麽考慮拆解\(n\)，設\(n=\prod_{i=1}^{k}a_i^{p_i}\)，對於\(n\)的每一個質因子，提一個出來，剩下的可以直接乘上去，即：

設\(y=\prod_{i=1}^k a_i^{p_i-1}\)，\(w=\prod_{i=1}^ka_i\)，式子就變成：
\[ \begin{align} s(n,m)&=y\cdot\sum_{i=1}^{m}\varphi(wi)\&=y\cdot \sum_{i=1}^{m}\varphi(\frac{w}{(i,w)})\cdot \varphi(i)\cdot (i,w) \end{align} \]
其中\((i,j)\)表示\(\gcd(i,j)\)。

然後由於：
\[ \sum_{d|n}\varphi(d)=n \]
可以反過來把後面的\(\gcd\)化一下，即：

當\(n=1\)時這個東西就變成了\(\sum_{i=1}^{m}\varphi(i)\)，杜教篩一下就好了。

#include<bits/stdc++.h>
//#include<tr1/unordered_map>
using namespace std;
//using namespace std :: tr1;

#define ll long long 
//#define ONLINE_JUDGE

#ifdef ONLINE_JUDGE
#define getchar() ((p1==p2&&(p2=(p1=buf)+fread(buf,1,1<<21,stdin)),p1==p2)?EOF:*p1++)
#endif

namespace fast_IO {
    char buf[1<<21],*p1=buf,*p2=buf;

    template <typename T> inline void read(T &x) {
        x=0;T f=1;char ch=getchar();
        for(;!isdigit(ch);ch=getchar()) if(ch=='-') f=-f;
        for(;isdigit(ch);ch=getchar()) x=x*10+ch-'0';x*=f;
    }
    template <typename T,typename... Args> inline void read(T& x,Args& ...args) {
        read(x),read(args...);
    }

    char buf2[1<<21],a[80];int p,p3=-1;

    inline void flush() {fwrite(buf2,1,p3+1,stdout),p3=-1;}
    template <typename T> inline void write(T x) {
        if(p3>(1<<20)) flush();
        if(x<0) buf2[++p3]='-',x=-x;
        do {a[++p]=x%10+48;} while(x/=10);
        do {buf2[++p3]=a[p];} while(--p);
        buf2[++p3]='\n';flush();
    }
    template <typename T,typename... Args> inline void write(T x,Args ...args) {
        write(x),write(args...);
    }
}

using fast_IO :: read;
using fast_IO :: write;
using fast_IO :: flush;

const int maxn = 3e6+10;
const int mod = 1e9+7;
const int N = 1e5+10;

int phi[maxn],pri[maxn],vis[maxn],tot,w[maxn],varphi[maxn];

void sieve() {
    phi[1]=1;w[1]=1;
    for(int i=2;i<maxn;i++) {
        if(!vis[i]) pri[++tot]=i,phi[i]=i-1,w[i]=i;
        for(int j=1;j<=tot&&i*pri[j]<maxn;j++) {
            vis[i*pri[j]]=1;
            if(i%pri[j]==0) {phi[i*pri[j]]=phi[i]*pri[j],w[i*pri[j]]=w[i];break;}
            phi[i*pri[j]]=phi[i]*(pri[j]-1);
            w[i*pri[j]]=w[i]*pri[j];
        }
    }
    for(int i=1;i<maxn;i++) varphi[i]=phi[i],phi[i]=(phi[i]+phi[i-1])%mod;
}

map<int,int> Phi,S[N];

int sum_phi(int x) {
    if(x<maxn) return phi[x];
    if(Phi[x]) return Phi[x];
    int res=1ll*x*(x+1)/2%mod,i=2;
    while(i<=x) {
        int pre=i;i=x/(x/i);
        res=(res-1ll*(i-pre+1)*sum_phi(x/i)%mod)%mod;i++;
    }return Phi[x]=(res+mod)%mod;
}

vector <int > Div;

int s(int n,int m) {
    if(!n) return 0;
    if(!m) return 0;
    if(n==1) return sum_phi(m);
    if(S[n][m]) return S[n][m];
    int y=n/w[n],res=0;
    for(int i=1;i*i<=w[n];i++) {
        if(w[n]%i) continue;
        int j=w[n]/i;
        res=(res+1ll*varphi[j]*s(i,m/i)%mod)%mod;
        if(i!=j) res=(res+1ll*varphi[i]*s(j,m/j)%mod)%mod;
    }
    return S[n][m]=1ll*res*y%mod;
}

int main() {
    sieve();
    int n,m;read(n,m);int ans=0;
    for(int i=1;i<=n;i++) ans=(ans+s(i,m))%mod;
    write(ans);
    flush();
    return 0;
}

[bzoj3512] DZY Loves Math IV