莫比乌斯反演第一题。莫比乌斯反演入门
数论题不多BB,直接推导吧。
首先,发现题目所求(ans=sum_{i=1}^nsum_{j=1}^m [gcd(i,j)=prime])
考虑反演,我们令(f(d))为(gcd(i,j)(iin[1,n],jin[1,m])=d)的个数,(F(n))为(d|gcd(i,j))的个数
即:
[f(d)=sum_{i=1}^nsum_{j=1}^m [gcd(i,j)=d]
]
[F(s)=sum_{s|d}f(d)=lfloorfrac{n}{s}
floor lfloorfrac{m}{s}
floor
]
根据莫比乌斯反演定理,有:
[f(n)=sum_{n|d} mu(lfloorfrac{d}{n}
floor)F(d)
]
所以开始化简,原式:
[ans=sum_{pin prime} sum_{i=1}^n sum_{j=1}^m[gcd(i,j)=p]
]
后面那个就是我们设的(f(p)),所以我们带入得:
[ans=sum_{pin prime} f(p)
]
用上面的公式反演一下:
[ans=sum_{pin prime} sum_{p|d}mu(lfloorfrac{d}{p}
floor)F(d)
]
这个不好处理,我们考虑把枚举(p)变成枚举(lfloorfrac{d}{p} floor),则有:
[ans=sum_{pin prime} sum_{d=1}^{min(lfloorfrac{n}{p}
floor,lfloorfrac{m}{p}
floor)}mu(d)F(dcdot p)
]
再根据(F(s))的性质代进去就有:
[ans=sum_{pin prime} sum_{d=1}^{min(lfloorfrac{n}{p}
floor,lfloorfrac{m}{p}
floor)}mu(d)lfloorfrac{n}{dcdot p}
floorlfloorfrac{m}{dcdot p}
floor
]
然后我们用(T)换掉(dcdot p),则有:
[ans=sum_{T=1}^{min(n,m)} sum_{p|T,pin prime} mu(lfloorfrac{T}{p}
floor)lfloorfrac{n}{T}
floorlfloorfrac{m}{T}
floor
]
把(lfloorfrac{n}{T} floorlfloorfrac{m}{T} floor)提到前面来就有:
[ans=sum_{T=1}^{min(n,m)} lfloorfrac{n}{T}
floorlfloorfrac{m}{T}
floorsum_{p|T,pin prime} mu(lfloorfrac{T}{p}
floor)
]
然后我们发现这个式子就是(O(n))的了,一看题目,多组询问!
结果还是套路,看到那个(lfloorfrac{n}{T} floorlfloorfrac{m}{T} floor)发现可以除法分块搞。
后面的式子稍加分析可以发现在筛法结束后枚举每个素数,再枚举(lfloorfrac{T}{p} floor)来统计。最后做一个前缀和即可。
CODE
#include<cstdio>
#include<cctype>
#define RI register int
using namespace std;
const int P=10000000;
int t,n,m,prime[P+5],cnt,miu[P+5]; long long sum[P+5],ans; bool vis[P+5];
class FileInputOutput
{
private:
#define tc() (A==B&&(B=(A=Fin)+fread(Fin,1,S,stdin),A==B)?EOF:*A++)
#define pc(ch) (Ftop<S?Fout[Ftop++]=ch:(fwrite(Fout,1,S,stdout),Fout[(Ftop=0)++]=ch))
#define S 1<<21
char Fin[S],Fout[S],*A,*B; int Ftop,pt[25];
public:
FileInputOutput() { A=B=Fin; Ftop=0; }
inline void read(int &x)
{
x=0; char ch; int flag=1; while (!isdigit(ch=tc())) flag=ch^'-'?1:-1;
while (x=(x<<3)+(x<<1)+(ch&15),isdigit(ch=tc())); x*=flag;
}
inline void write(long long x)
{
if (!x) return (void)(pc(48),pc('
')); RI ptop=0;
while (x) pt[++ptop]=x%10,x/=10; while (ptop) pc(pt[ptop--]+48); pc('
');
}
inline void Fend(void)
{
fwrite(Fout,1,Ftop,stdout);
}
#undef tc
#undef pc
#undef S
}F;
#define Pi prime[j]
inline void Euler(void)
{
vis[1]=miu[1]=1; RI i,j; for (i=2;i<=P;++i)
{
if (!vis[i]) prime[++cnt]=i,miu[i]=-1;
for (j=1;j<=cnt&&i*Pi<=P;++j)
{
vis[i*Pi]=1; if (i%Pi) miu[i*Pi]=-miu[i]; else break;
}
}
for (j=1;j<=cnt;++j) for (i=1;i*Pi<=P;++i) sum[i*Pi]+=miu[i];
for (i=1;i<=P;++i) sum[i]+=sum[i-1];
}
#undef Pi
inline int min(int a,int b)
{
return a<b?a:b;
}
inline void swap(int &a,int &b)
{
int t=a; a=b; b=t;
}
int main()
{
//freopen("CODE.in","r",stdin); freopen("CODE.out","w",stdout);
for (Euler(),F.read(t);t;--t)
{
F.read(n); F.read(m); ans=0; if (n>m) swap(n,m);
for (RI l=1,r;l<=n;l=r+1)
{
r=min(n/(n/l),m/(m/l)); ans+=1LL*(n/l)*(m/l)*(sum[r]-sum[l-1]);
}
F.write(ans);
}
return F.Fend(),0;
}