大致题意: 求(sum_{i=1}^nsum_{j=1}^mphi(ij))。
推式子
好妙的一道题。。。
考虑(n)和(m)的大小相差很多,因此我们可以考虑枚举(n),即令:
[S_n(m)=sum_{i=1}^mphi(ni)
]
这样有什么好处呢?最明显的一个就是,我们可以把原本无比棘手的含有两个变量的(phi)给拆开。
假设(n=prod P_i^{k_i}),则(phi(n)=phi(prod P_i) imes prod P_i^{k_i-1})。
不妨令(p=prod P_i^{k_i-1},q=prod P_i),则:
[S_n(m)=p imes sum_{i=1}^mphi(qi)=p imes sum_{i=1}^mphi(frac q{gcd(q,i)})phi(i) imes gcd(q,i)
]
此时式子中虽然有一个(gcd),但并不能无脑上反演(说不定可以,反正我不会)。
看到式子中有这么多(phi),且(phi(frac q{gcd(q,i)}))的分母中也刚好有(gcd(q,i)),因此我们可以利用(phi*I=Id),得到:
[S_n(m)=p imes sum_{i=1}^mphi(frac q{gcd(q,i)})phi(i) imessum_{d|gcd(q,i)}phi(d)
]
把第一个(phi)和最后一个(phi)合并,也就得到:
[S_n(m)=p imessum_{i=1}^mphi(i)sum_{d|gcd(q,i)}phi(frac qd)
]
套路地枚举(d)得到:
[S_n(m)=p imes sum_{d|q}phi(frac qd)sum_{i=1}^{lfloorfrac md
floor}phi(di)=p imes sum_{d|q}phi(frac qd)S_d(lfloorfrac md
floor)
]
这个式子可以递归求解,于是就可以上杜教筛了。
代码
#pragma GCC optimize(2)
#pragma GCC optimize("inline")
#include<bits/stdc++.h>
#define Tp template<typename Ty>
#define Ts template<typename Ty,typename... Ar>
#define Reg register
#define RI Reg int
#define Con const
#define CI Con int&
#define I inline
#define W while
#define N 100000
#define X 1000000007
#define Inc(x,y) ((x+=(y))>=X&&(x-=X))
using namespace std;
int n,m,sm;
class DuSieve//杜教筛
{
private:
#define SZ 100000
#define V(x) ((x)<=sm?f1[x]:f2[m/(x)])
int Pt,P[SZ+5],phi[SZ+5],sphi[N+5],f1[SZ+5],f2[SZ+5];map<int,int> s[SZ+5];
I int SPhi(CI x)//求Phi
{
if(x<=SZ) return sphi[x];if(V(x)) return V(x);RI l,r,t=1LL*x*(x+1)/2%X;
for(l=2;l<=x;l=r+1) r=x/(x/l),t=(t-1LL*(r-l+1)*SPhi(x/l)%X+X)%X;return V(x)=t;
}
public:
I DuSieve()//预处理
{
phi[1]=1;for(RI i=2,j,t;i<=SZ;++i)
for(!P[i]&&(phi[P[++Pt]=i]=i-1),j=1;j<=Pt&&(t=i*P[j])<=SZ;++j)
if(P[t]=1,i%P[j]) phi[t]=phi[i]*(P[j]-1);else {phi[t]=phi[i]*P[j];break;}
for(RI i=1;i<=SZ;++i) sphi[i]=(sphi[i-1]+phi[i])%X;
}
I int S(RI n,CI m)//求S
{
if(!m) return 0;if(n==1) return SPhi(m);if(m==1) return phi[n];//边界
if(s[n][m]) return s[n][m];RI i,k=n,p=1,q=1,t=0;//记忆化
for(i=1;P[i]*P[i]<=k;++i) if(!(k%P[i])) for(k/=P[i],q*=P[i];!(k%P[i]);k/=P[i],p*=P[i]);//质因数分解
for(k^1&&(q*=k),i=1;i*i<=q;++i) !(q%i)&&//暴力枚举q的因数
(t=(1LL*phi[q/i]*S(i,m/i)+t)%X,i^(q/i)&&(t=(1LL*phi[i]*S(q/i,m/(q/i))+t)%X));
return s[n][m]=1LL*p*t%X;//注意最后乘上p
}
}D;
int main()
{
RI i,j,t=0;for(scanf("%d%d",&n,&m),sm=sqrt(m),i=1;i<=n;++i) Inc(t,D.S(i,m));//暴力枚举n
return printf("%d
",t),0;
}