@
参考文献
咕咕日报：https://www.luogu.org/blog/lx-2003/mobius-inversion
各题题解
定义
对于此篇博客的一些同样的定义，在此给出：
(a|b)表示(a\%b=0)，(a⊥b)表示(gcd(a,b)=1)，([P])表示的是当(P)为真时，值为(1)，(P)为假时，值为(0)。
且对于下面的文章，两个数论函数(f,g)的狄利克雷卷积形式写成(f*g)。
狄利克雷卷积以及各种性质
数论函数
卷积？
卷积难道不是这样的形式吗（这里的(f(n))表示的是次数为(n)的系数）：
(f(n)=sumlimits_{i=0}^{n}a(i)b(n-i))
当然，这是普通卷积，适用于所有函数。
而狄利克雷卷积不一样，他是作用于数论函数的。
数论函数是什么呢？就是定义域在正整数域，值域在实数域的函数。
比如斐波那契数列、(φ)等。
形式
对于三个数论函数（(f(n))表示的是第(n)项的值）：
(f(n)=sumlimits_{d|n}a(d)b(frac{n}{d}))
这个形式就是狄利克雷卷积形式。
性质
-
交换律，即：(f*g=g*f)，这个显然。
-
结合律，(f*g*t=f*(g*t))。
证明：(sumlimits_{i*j*k}f(i)g(j)t(k)=sumlimits_{i*(j*k)=n}f(i)(g(j)t(k)))
-
分配率，((f+g)*h=f*h+g*h)，这个的话你只要把加法拆开，想想就知道了。
-
系数提出，也就是两个函数想乘，如果一个函数的所有值域都有一个公因子的话，那么可以将其提出，形式写为：((xf)*g=x(f*g))。
积性函数
也许你经常会看到某某某函数为完全积性函数。
但是却不知道他是什么，对吧。
这里我们给出定义：
当(f)是数论函数且(x⊥y)时(f(xy)=f(x)f(y))，那么(f)就是积性函数。
性质1
当(f)是积性函数时，(f(1)=1)否则(f)的值全部为(0)。
证明：如果(f(x))的值不为(0)，且(f(1)≠1)，因为(f(1)*f(1)=f(1))，所以(f(1)=0)，而(f(x)*f(1)=f(x))，当(f(x)>0)且(f(1)=0)，则(f(x)*f(1)=0)，不成立，所以(f(1)=1)。
性质2
两个积性函数的狄利克雷卷积也是积性函数。
下面证明来自咕咕日报（注意要利用到一个定理，就是：(n⊥m,a|n,b|m)，那么(a⊥b)）。
注：但是两个完全积性函数的狄利克雷卷积不一定是完全积性函数。
一些完全积性函数
这里我们给出一些常见的完全积性函数。
(ϵ)函数，又叫单位元。
(ϵ(i)=[i=1])（也就是只有(ϵ(1)=1)，其他都为(0)）
（性质：(ϵ*f=f)，这个想想就知道了吧QMQ）
(id)函数，(id(i)=i)。
(1)函数，(1(i)=1)。
逆函数
逆函数也是数论函数。（理论上来讲数论函数的除法就是乘上逆函数。）
对于一个数论函数(f)，其逆函数是(g)，那么(f*g=ϵ)。
相当于数论函数中的逆元。
但是问题又来了，如何求逆函数呢？
我们这里给出公式（还是咕咕日报的）：
(g(n)=frac{1}{f(1)}([n=1]-sumlimits_{i|n,i≠1}f(i)g(frac{n}{i})))
可以理解为是一种递推式。
把这个定义套进去：
(sumlimits_{i|n}f(i)g(frac{n}{i}))
(=f(1)g(n)+sumlimits_{i|n,i≠1}f(i)g(frac{n}{i}))
(=[n=1])
其实可以证明一个函数的逆函数只有一个（下面的证明不包括那个全部为(0)的数论函数，那玩意没有逆函数吧！！！QAQ）。（应该只有一个吧QAQ）
对于(f(1))，(g[1])是固定的，而对于(g[1])~(g[i-1])固定的话，(g[i])也是固定的。（这个就是不用公式也可以明白吧。）
所以逆函数是固定唯一的。
性质1
积性函数的逆函数也是积性函数。
来自咕咕日报：
这里解释一下。
他没有证(0)函数（全部为(0)的那个函数，以后我们讨论逆函数都不讨论他）。
同时为什么他可以直接化成这个式子呢，因为除了(0)函数的积性函数的(f(1)=1)，而对于(nm>1)的情况，([nm=1])肯定为(0)，所以直接化成了那个样子。
莫比乌斯反演
因子包含
对于(f)函数，我们定义(F=f*1)，即
(F(n)=sumlimits_{i|n}f(i))。
那么(f(n)=sumlimits_{i|n}mu(frac{n}{i})*F(i))。
这就是因子包含的反演公式，你也许会觉得这个公式很废柴，但是下面会讲的。（这个公式其实是真的废柴，反正我做了这么久的题目都很少套用这个，不过如果你不会这个的证明的话莫比乌斯反演还是有点难度的。）
狄利克雷卷积方式证明因子包含公式（适合较理性的同学证明）
终于到了正题了。
我们定义(mu)为(1)函数的逆。
那么(g=f*1)，那么就有(f=f*1*mu=g*mu)。
换成式子形式就是：
(f(n)=sumlimits_{i|n}g(i)*mu(frac{n}{i}))
证毕。
真的感觉狄利克雷卷积的证明方式简单易懂。
杨辉三角形与容斥原理
下面证明需要。
对于杨辉三角形第(i)行第(j)列就是(f(i,j))，根据定义就是(f(i,j)=f(i-1,f)+f(i-1,j-1))，我们默认(f(i,1)=f(i,i)=1)。
我们发现其实(f(i,j)=C_{i-1}^{j-1})。
用数学归纳法来证明，对于(i-1)行都满足这个式子，那么(i)行满不满足？
(C_{i-1}^{j}+C_{i-1}^{j-1}=frac{(i-1)!}{j!(i-j-1)!}+frac{(i-1)!}{(j-1)!(i-j)!})
(=frac{(i-1)!(i-j)}{j!(i-j)!}+frac{(i-1)!j}{j!(i-j)!}=frac{i!}{i!(i-j)!}=C_{i}^j)
那么我们知道了杨辉三角形的性质有什么用呢？
对于杨辉三角形第(i)层，我们选出下标是奇数的，发现他们映射到上一层刚好覆盖。
偶数的也一样。
那么只要奇数项为(+)，然后偶数项为(-)，一加就为(0)了。（除了三角形顶端）
μ函数的定义
根据逆函数的定义，我们易知(mu(1)=1,mu(2)=-1,mu(3)=-1...)。
仔细观察，(mu(i))当(i)为质数时为(-1)，(i=1)时为(1)。
我们还要得到更普遍的定义。
仔细一看这个式子（(n>1)）：
(g(n)=-sumlimits_{i|n,i≠1}g(frac{n}{i}))
(g(n)=-sumlimits_{i|n,i≠n}g(i))
这个不就是查自己所有因子和(1)除自己以外的(g)值和的负数吗。
那么有(2)质数因子的就是(1)，那么(3)个质数因子用容斥就可以知道是(-1)。
那么我们用数学归纳法来证明。
对于有(i(i\%2=0))个质数因子的数（假设前面的奇数个质数因子的数都为(-1)，偶数个质数因子的数都为(1)。）
(C_{i}^{0}-C_{i}^{1}+C_{i}^{2}-C_{i}^{3}+...-C_{i}^{i-1}=(C_{i}^{0}-C_{i}^{1}+C_{i}^{2}-C_{i}^{3}+...-C_{i}^{i-1}+C_{i}^{i})-C_{i}^{i}=0-C_{i}^{i}=-1)
反之偶数是(1)（因为偶数的话，最后一项是(-)，要(+)回来）。
以上证明利用杨辉三角形性质。
那么如果有的因子次数是(2)或以上呢，我们先证数字是质数的二次幂的数字，如：(p_0^{2})。
只有两个因子：(1,p_0)，一加为(0)。
那么数学归纳法走起。
证明(p_{0}^{i}=0)，当(p_{0}^{2})~(p_0^{i-1})都是(0)
那么真实有效的就只有(1)和(p_0)，那么就是(0)啦。
证毕。
而且因为(1)是积性函数，所以(mu)也是积性函数。
所以任何因子有质数次幂的数字都可以表示为一个数字乘以质数次幂。
那么也就都是(0)了。
至此，我们得到了(mu)的值。
(mu(i)=left{egin{matrix} 1（当i有奇数个质数因子）\ -1（当i有偶数个质数因子）\ 0（当i的某个质数因子的质数大于1） end{matrix} ight.)
容斥原理证明
容斥原理证明原理有个缺陷就是很难推式子，但是证明正确性确实挺好的。
不过经过这次证明我的容斥进步了不少
我们再把式子看一遍：
(F(n)=sumlimits_{i|n}f(i))
(f(n)=sumlimits_{i|n}mu(frac{n}{i})*F(i))
那么我们容斥就是通过容斥把(F(n))中的(f(n))保留，同时把(f(i)(i≠n))全部删去。
而我们继续观察下面的式子，这个式子我们再把其转化一下：
(f(n)=sumlimits_{i|n}mu(i)*F(frac{n}{i}))
对于(f(i)(i|n))，我们设(frac{n}{i})拆成(p_0^{a_0}p_1^{a_1}...p_i^{a_i})（老样子，(p)都是质数），那么这个式子的本质就是减去(0)个质数的情况，加上(1)个质数的情况，减去(2)个质数的情况，加上(3)个质数的情况，其实也就是讨论(p_0p_1p_2...p_i)的这个数字，因为如果质数质数大于(1)的话(mu)为(0)，所以其实是可以把指数去掉。
再仔细观察式子不就是求(f(i))的系数吗，(f(i))的系数加减可以看成是(C_{i}^0-C_i^1+...=0)。
当然，有个例外，(f(n)=C_0^0=1)，所以证毕。
μ函数的性质
由于(mu*1=ϵ)
(sumlimits_{i|n}mu(i)=[n=1])
这个性质十分重要！
倍数包含
(F(i)=sumlimits^{n}_{i|j}f(j))
(f(i)=sumlimits^{n}_{i|j}mu(frac{j}{i})f(j))
这个你只要理解了因子包含的容斥证明就可以证了。
数论分块
对于(frac{n}{i}(1<=i<=n,i∈Z^+))
那么这个有多少取值呢。
可以证明最多(2left lceil sqrt{n} ight ceil)。
当(i<sqrt{n})，有(sqrt{n})个取值
而(i=sqrt{n},frac{n}{i}=sqrt{n})，所以(i>=sqrt{n})也只有(sqrt{n})个取值。
所以理论上(frac{n}{i})的值应该有一段的(i)相同。
但是如果求出这一段的(r)呢？就是(left lfloorfrac{n}{left lfloor frac{n}{i} ight floor} ight floor)。
我们可以证一下。
对于(left lfloor frac{n}{i} ight floor)，我们可以看成是找到最大的整数(x)满足(x*i<=n)，同时(left lfloor frac{n}{x} ight floor=i)
为什么
对于(x=left lfloor frac{n}{i} ight floor)。
那么((x+1)*i=x*i+i)，由于(x=left lfloor frac{n}{i} ight floor)，所以(x*i<=n)，且(n-x*i<i)，所以((x+1)*i>n)，所以(leftlfloorfrac{n}{x} ight floor=i)
那么对于我们要看(left lfloor frac{n}{i} ight floor=y)的最大的(i)，那么其实就可以化成(left lfloor frac{n}{y} ight floor=i)，也就是最大的(i)了，而且由上面的式子可以得到(left lfloor frac{n}{i+1} ight floor<y)。
所以就是(left lfloorfrac{n}{left lfloor frac{n}{i} ight floor} ight floor)
那么这个有什么用后面的第一题就会讲到了。
练习
1
我们发现([a,b],[c,d])其实可以拆成四个区间。
(([1,b],[1,d])-([1,a-1],[1,d])-([1,b],[1,c-1])+([1,a-1],[1,c-1]))
开始推公式吧。其实这应该是例题来的
(sumlimits_{i=1}^{b}sumlimits_{j=1}^{d}[gcd(i,j)=k])
我们整体除以(k)，那么就是
(sumlimits_{i=1}^{frac{b}{k}}sumlimits_{j=1}^{frac{d}{k}}[gcd(i,j)=1])
(=sumlimits_{i=1}^{frac{b}{k}}sumlimits_{j=1}^{frac{d}{k}}[gcd(i,j)=1])
咦后面那个是不是可以换成(mu)函数呀。
(=sumlimits_{i=1}^{frac{b}{k}}sumlimits_{j=1}^{frac{d}{k}}sumlimits_{t|gcd(i,j)}mu(t))
那怎么求。
莫比乌斯反演推公式最重要的一步之一。
调换加法顺序。
我们先设数字，不然有点难写。
(b’=frac{b}{k},d'=frac{d}{k})。
下面证明用(b,d)代替(b',d')。
(sumlimits_{t=1}^{min(b,d)}mu(t)(sumlimits_{i=1}^{frac{b}{t}}sumlimits_{j=1}^{frac{d}{t}}1))
我们突然发现后面的部分可以(O(1))求出来，
但是前面循环怎么办呢。
爆算岂不是(O(n^2))
我们突然发现如果(frac{b}{t},frac{d}{t})不变，后面就不变，所以我们只要对(mu)做一遍前缀和就行了。
而(frac{b}{t},frac{d}{t})合起来也只有(O(sqrt{n}))和取值。
所以时间复杂度就是(O(nsqrt{n}))
其实即使不拆应该也可以算，只不过应该要处理(ceil)和(floor)，太麻烦了。。。
#include<cstdio>
#include<cstring>
#define N 51000
using namespace std;
typedef long long LL;
LL u[N];
int pr[N],top;bool v[N];
void first_do(int x)//u是积性函数
{
u[1]=1;
for(int i=2;i<=x;i++)
{
if(!v[i])v[i]=1,u[i]=-1,pr[++top]=i;
for(int j=1;j<=top && pr[j]*i<=x;j++)
{
v[pr[j]*i]=1;
u[pr[j]*i]=-u[i];
if(i%pr[j]==0)
{
u[pr[j]*i]=0;
break;
}
}
}
for(int i=1;i<=x;i++)u[i]+=u[i-1];//前缀和
}
inline LL mymin(LL x,LL y){return x<y?x:y;}
inline LL mymax(LL x,LL y){return x>y?x:y;}
LL calc(LL n,LL m)
{
n>m?(n^=m^=n^=m):0;
LL ans=0;
int last=0;
for(int i=1;i<=n;i=last+1)
{
last=mymin(n/(n/i),m/(m/i));
ans+=(u[last]-u[i-1])*(n/i)*(m/i);//数论分块
}
return ans;
}
LL solve(int a,int b,int c,int d,int k)
{
a--;c--;a/=k;b/=k;c/=k;d/=k;
return calc(b,d)-calc(a,d)-calc(b,c)+calc(a,c);//四个区间
}
int main()
{
first_do(50000);
int T;scanf("%d",&T);
while(T--)
{
int a,b,c,d,k;scanf("%d%d%d%d%d",&a,&b,&c,&d,&k);
printf("%lld
",solve(a,b,c,d,k));
}
return 0;
}
2
那么这道题目又是什么意思，一开始我们以为很难，但是仔细一看，秒呀，原来(1≤x,y≤n)，原来是同一个数字呀。
注：(prim)是质数集。
那么不就可以套用(φ)函数了吗，(φ(i))的定义是与([1,i])内与(i)互质的数的个数。
化式子((n<=m))。
(sumlimits_{d∈prim}^{n}sumlimits_{i=1}^{frac{n}{d}}sumlimits_{j=1}^{frac{n}{d}}[gcd(i,j)=1])
其实这个后面也有一道加强版
但是其实我们是可以用(φ)函数(O(1))求出后面的数字，就是用(φ)做一遍前缀和(g)。
那么就是(2g(frac{n}{d})-1)，后面那个(-1)是怎么回事呢，我们是把((i,j)(i<j))扩展到了所有数对，但是((1,1))不一样啊。
所以我们可以爆枚素数。
不就是一遍线性筛的事情吗。
(O(n))。
#include<cstdio>
#include<cstring>
#define N 11000000
#define M 1100000
using namespace std;
typedef long long LL;
LL phi[N];bool v[N];
int pr[M],top;
void first_do(int x)
{
phi[1]=1;
for(int i=2;i<=x;i++)
{
if(!v[i])pr[++top]=i,phi[i]=i-1;
for(int j=1;j<=top && pr[j]*i<=x;j++)
{
v[pr[j]*i]=1;
if(i%pr[j]==0){phi[i*pr[j]]=phi[i]*pr[j];break;}
phi[i*pr[j]]=phi[i]*phi[pr[j]];
}
}
for(int i=2;i<=x;i++)phi[i]+=phi[i-1];//前缀和
for(int i=2;i<=x;i++)phi[i]=phi[i]*2-1;//处理
}
int n;
int main()
{
scanf("%d",&n);
first_do(n);
LL ans=0;
for(int i=1;i<=top;i++)
{
int x=pr[i];
ans+=phi[n/x];
}//O(n)出奇迹
printf("%lld
",ans);
return 0;
}
3
完全平方数？
这道题目应该是考你容斥的。
设(i)为(p_0^{a_0}p_1^{a_1}p_2^{a_2}...p_i^{a_{i}})。
然后把(a>1)的提出来化成(i)的一个因子(b_0^2,b_1^2b_2^2...b_j^2)。
那么我们就是想让统计区间内提出来后因子为(1)的保留，因子大于(1)的消掉。
那么就是(sumlimits_{i=1}^{i^2<=n}mu(i)frac{n}{i^2})。
预处理(mu)。
那么现在问题就是我们怎么知道(n)呢？二分枚举出奇迹！！！
时间复杂度：(O(Tsqrt{n}logn))
#include<cstdio>
#include<cstring>
#include<cmath>
#define N 51000
using namespace std;
typedef long long LL;
int u[N];bool v[N];int pr[N],top;
void first_do(int x)
{
u[1]=1;
for(int i=2;i<=x;i++)
{
if(!v[i])pr[++top]=i,u[i]=-1;
for(int j=1;j<=top && pr[j]*i<=x;j++)
{
v[pr[j]*i]=1;
if(i%pr[j]==0){u[pr[j]*i]=0;break;}
u[pr[j]*i]=-u[i];
}
}
}
int check(int x)
{
int ed=sqrt(x)+1,ans=0;
for(int i=1;i<=ed;i++)ans+=u[i]*(x)/(i*i);
return ans;
}
int main()
{
first_do(50000);
int T;scanf("%d",&T);
while(T--)
{
int n;scanf("%d",&n);
int l=1,r=2000000000,mid,ans;
while(l<=r)
{
mid=(r-l)/2+l;
if(check(mid)>=n)ans=mid,r=mid-1;
else l=mid+1;
}
printf("%d
",ans);
}
return 0;
}
4
这道题目我们还是继续化简式子（一下证明(n≤m)）。
(sumlimits_{t∈prim}^{n}sumlimits_{i=1}^{n}sumlimits_{j=1}^{m}[gcd(i,j)=d])
(sumlimits_{t∈prim}^{n}sumlimits_{i=1}^{frac{n}{d}}sumlimits_{j=1}^{frac{m}{d}}sumlimits_{k|gcd(i,j)}mu(k))
调换顺序
(sumlimits^{n}_{t∈prim}sumlimits_{k=1}^{frac{n}{d}}mu(k)sumlimits^{frac{n}{kt}}_{i=1}sumlimits^{frac{m}{kt}}_{j=1}1)
咦，后面的是不是可以(O(1))算。
但是你是不是还是要枚举(t)？看看多组数据，直接爆炸！！！Boom！
所以我们这里引入这种莫反的另外一个骚操作！
设(T=tk)。
(sumlimits_{T=1}^{n}sumlimits_{t∈prim,t|T}mu(frac{T}{t})sumlimits_{i=1}^{frac{n}{T}}sumlimits_{j=1}^{frac{n}{T}}1)
这不是一样吗QAQ。
一样吗-_-
(sumlimits_{T=1}^{n}sumlimits_{i=1}^{frac{n}{T}}sumlimits_{j=1}^{frac{n}{T}}sumlimits_{t∈prim,t|T}mu(frac{T}{t}))
原来能这么化吗！！！！！！！！！！
我们只要预处理后面那个就行了！！！！！
而后面那个的话预处理应该是差不多还行吧QAQ，我也不知道也。
不过应该是接近一亿的。
然后前面的数论分块就行了。
时间复杂度(O(?+Tsqrt{n}))
#include<cstdio>
#include<cstring>
#define N 11000000
#define M 2100000
using namespace std;
typedef long long LL;
LL u[N],g[N];
int pr[M],top;bool v[N];
void first_do(int x)
{
u[1]=1;
for(int i=2;i<=x;i++)
{
if(!v[i])pr[++top]=i,u[i]=-1;
for(int j=1;j<=top && pr[j]*i<=x;j++)
{
v[pr[j]*i]=1;
if(i%pr[j]==0)
{
u[i*pr[j]]=0;
break;
}
u[i*pr[j]]=-u[i];
}
}
for(int i=1;i<=top;i++)
{
for(int j=1;pr[i]*j<=x;j++)g[pr[i]*j]+=u[j];
}
for(int i=1;i<=x;i++)g[i]+=g[i-1];
}
inline int mymin(int x,int y){return x<y?x:y;}
LL calc(int x,int y)
{
x>y?(x^=y^=x^=y):0;
LL ans=0;
int last=0;
for(int i=1;i<=x;i=last+1)
{
last=mymin(x/(x/i),y/(y/i));
ans+=(g[last]-g[i-1])*(LL)(x/i)*(LL)(y/i);
}
return ans;
}
int main()
{
first_do(10000000);
int T;scanf("%d",&T);
while(T--)
{
int n,m;scanf("%d%d",&n,&m);
printf("%lld
",calc(n,m));
}
return 0;
}
5
好题。
一道很成功的将数论和数据结构结合的题目（设(n<=m)）。
设(d(i))为一个数字的约数和。
那么这个函数很明显是(id*1)，所以是一个积性函数，可以用线性筛筛。
然后开始列式子，我们假设没有(a)的限制。
(sumlimits_{i=1}^{n}sumlimits_{j=1}^{m}d(gcd(i,j)))
(sumlimits_{d=1}^{n}sumlimits_{i=1}^{frac{n}{d}}sumlimits_{j=1}^{frac{m}{d}}d(d)[gcd(i,j)=1])
(sumlimits_{d=1}^{n}d(d)sumlimits_{i=1}^{frac{n}{d}}sumlimits_{j=1}^{frac{m}{d}}sumlimits_{k|gcd(i,j)}mu(k))
(sumlimits_{d=1}^{n}d(d)sumlimits_{k=1}^{frac{n}{d}}mu(k)frac{n}{kd}*frac{m}{kd})
我们再设(T=kd)
(sumlimits_{T=1}^{n}sumlimits_{k|T}d(frac{T}{k})mu(k)(frac{n}{T})(frac{m}{T}))
(sumlimits_{T=1}^{n}(frac{n}{T})(frac{m}{T})sumlimits_{k|T}d(frac{T}{k})mu(k))
仔细一看，前面可以整除分块了呢！！！
但是后面QAQ。
由于我们需要数论分块，我们需要后面这坨东西的前缀和。
但是我们发现，只有当(frac{T}{k}>=a)的时候，(d(frac{T}{k})mu(k))才可以造成影响。
所以我们应当把询问拿(a)排个序。
然后对于目前的(a)，暴力查哪个(d)又可以了（这个按(d)值排个序），然后继续暴力队可以更新的(T)拿去更新，时间复杂度：(O(nlogn))，完美！！！
不对，你前缀和怎么查。。。
所以我们用树状数组解决。
时间复杂度：(O(Qsqrt{n}+nlog^2n))。
#include<cstdio>
#include<cstring>
#include<algorithm>
#define N 110000
using namespace std;
int a[N],mm=100000/*max*/,mod;
inline int lowbit(int x){return x&-x;}
int getsum(int x)
{
int ans=0;
while(x)
{
ans+=a[x];
x-=lowbit(x);
}
return ans;
}
void change(int x,int k)
{
while(x<=mm)
{
a[x]+=k;
x+=lowbit(x);
}
}
//树状数组
struct node
{
int x,y;/*约数和*/
}f[N];
inline bool cmp1(node x,node y){return x.x<y.x;}
int u[N],low[N];
bool v[N];int pr[N],top;
void ksc(int x,int y)//x^k<=y
{
int now=1;
while(now<=y/x)now*=x,f[now].x=f[now/x].x*x+1;
}
void first_do(int x)
{
u[1]=1;f[1].x=1;
for(int i=2;i<=x;i++)
{
if(!v[i])pr[++top]=i,ksc(i,x),low[i]=i,u[i]=-1;
for(int j=1;pr[j]*i<=x && j<=top;j++)
{
v[pr[j]*i]=1;
if(i%pr[j]==0)
{
u[i*pr[j]]=0;
f[i*pr[j]].x=f[i/low[i]].x*f[pr[j]*low[i]].x;
low[i*pr[j]]=low[i]*pr[j];
break;
}
u[i*pr[j]]=-u[i];
f[i*pr[j]].x=f[i].x*f[pr[j]].x;
low[i*pr[j]]=pr[j];
}
}
}//线性筛
int T;
struct question
{
int n,m,a,id;
}qu[N];int zans[N];
bool cmp(question x,question y){return x.a<y.a;}
inline int mymin(int x,int y){return x<y?x:y;}
int solve(int x,int y)
{
x>y?(x^=y^=x^=y):0;
int last;
int ans=0;
for(int i=1;i<=x;i=last+1)
{
last=mymin(x/(x/i),y/(y/i));
ans+=(getsum(last)-getsum(i-1))*(x/i)*(y/i);
}
return ans;
}
int main()
{
first_do(mm);
for(int i=1;i<=mm;i++)f[i].y=i;
sort(f+1,f+mm+1,cmp1);
int T;scanf("%d",&T);
for(int i=1;i<=T;i++){scanf("%d%d%d",&qu[i].n,&qu[i].m,&qu[i].a);qu[i].id=i;}
sort(qu+1,qu+T+1,cmp);//排序
int j=1;
for(int i=1;i<=T;i++)
{
while(f[j].x<=qu[i].a && j<=mm)
{
for(int k=f[j].y;k<=mm;k+=f[j].y)change(k,f[j].x*u[k/f[j].y]);
j++;
}
zans[qu[i].id]=solve(qu[i].n,qu[i].m);
}
for(int i=1;i<=T;i++)printf("%d
",zans[i]&(~(1<<31)));
return 0;
}
6
下面都是默认(n<=m)
这道题目是(lcm(i,j))？
我有点慌张。
但是一想(lcm(i,j)=frac{ij}{gcd(i,j)})。
于是又开始做了起来。
(sumlimits_{d=1}^{n}frac{1}{d}sumlimits_{i=1}^{frac{n}{d}}sumlimits_{i=1}^{frac{m}{d}}i*j*d^2[gcd(i,j)=1])
注：(i,j)除了(d)，所以要乘回来，变成了(d^2)。
(sumlimits_{d=1}^{n}dsumlimits_{i=1}^{frac{n}{d}}sumlimits_{i=1}^{frac{m}{d}}i*jsumlimits_{k|gcd(i.j)}mu(k))
(sumlimits_{d=1}^{n}dsumlimits_{k=1}^{frac{n}{d}}mu(k)k^2sumlimits_{i=1}^{frac{n}{dk}}isumlimits_{j=1}^{frac{m}{dk}}j)
然后设(T=dk)
(sumlimits_{T=1}^{n}sumlimits_{d|T}dfrac{T^2}{d^2}mu(frac{T}{d})sumlimits_{i=1}^{T}isumlimits_{j=1}^{T}j)
设(g(n)=sumlimits_{i=1}^{n}i=frac{n(n+1)}{2})
(sumlimits_{T=1}^{n}T^2g(frac{n}{T})g(frac{m}{T})sumlimits_{d|T}frac{mu(frac{T}{d})}{d})
这个式子确确实实是对的，但是(nlogn)预处理太顶了。
所以我们需要观察我们在那里漏了什么东西。
其实我们发现对于
(sumlimits_{d=1}^{n}dsumlimits_{k=1}^{frac{n}{d}}mu(k)k^2sumlimits_{i=1}^{frac{n}{dk}}isumlimits_{j=1}^{frac{m}{dk}}j)
我们设(sum(n,m)=sumlimits_{d=1}^{n}mu(d)d^2sumlimits_{i=1}^{frac{n}{d}}isumlimits_{j=1}^{frac{m}{d}}j)。
这个式子，只要我们处理除了(mu)，就可以数论分块(O(sqrt{n}))解决。
式子化为
(sumlimits_{d=1}^{n}d*sum(frac{n}{d},frac{m}{d}))
也是数论分块。
那么就可以(O(n))做出来啦！
#include<cstdio>
#include<cstring>
#define N 11000000
#define NN 3100000
using namespace std;
typedef long long LL;
LL u[N],mod=20101009;
bool v[N];
int pr[NN],top;
void pre_do(int x)
{
u[1]=1;
for(int i=2;i<=x;i++)
{
if(!v[i])u[i]=-1,pr[++top]=i;
for(int j=1;j<=top && i<=x/pr[j];j++)
{
v[i*pr[j]]=1;
if(i%pr[j]==0)
{
u[i*pr[j]]=0;
break;
}
u[i*pr[j]]=-u[i];
}
}
for(int i=1;i<=x;i++)u[i]=(u[i-1]+((LL)i*(LL)i)%mod*u[i])%mod,u[i]=(u[i]+mod)%mod;
}
template <class T>
inline T mymin(T x,T y){return x<y?x:y;}
LL calc(LL x,LL y)
{
x>y?x^=y^=x^=y:0;
int last=0;
LL ans=0;
for(int i=1;i<=x;i=last+1)
{
last=mymin(x/(x/i),y/(y/i));
LL xx=x/i,yy=y/i;
ans+=((((xx*(xx+1)/2)%mod)*((yy*(yy+1)/2)%mod)+mod)%mod)*(u[last]-u[i-1]);
ans%=mod;
}
return (ans+mod)%mod;
}
LL solve(LL x,LL y)
{
x>y?x^=y^=x^=y:0;
int last=0;
LL ans=0;
for(int i=1;i<=x;i=last+1)
{
last=mymin(x/(x/i),y/(y/i));
ans+=(((LL)(last-i+1)*(last+i)/2)%mod)*calc(x/i,y/i);
ans%=mod;
}
return ans;
}
int main()
{
pre_do(10000000);
int n,m;scanf("%d%d",&n,&m);
printf("%lld
",(solve(n,m)+mod)%mod);
return 0;
}
7
这道题目也是一道套路题呀。
(n<=m)
(sumlimits_{d=1}^{n}sumlimits_{i=1}^{n/d}sumlimits_{j=1}^{m/d}[gcd(i,j)==1]φ(i)φ(j))
(sumlimits_{d=1}^{n}sumlimits_{i=1}^{n/d}sumlimits_{j=1}^{m/d}φ(i)φ(j)sumlimits_{k|gcd(i,j)}mu(k))
(sumlimits_{d=1}^{n}sumlimits_{k=1}^{n/d}mu(k)sumlimits_{k|i}^{n/d}φ(i)sumlimits_{k|j}^{m/d}φ(j))
(sumlimits_{d=1}^{n}sumlimits_{k=1}^{n/d}mu(k)sumlimits_{i=1}^{n/dk}φ(ik)sumlimits_{j=1}^{m/dk}φ(jk))
然后我们设(g(n,k))表示的是(sumlimits_{i=1}^{n}φ(ik))。
那么这个是可以(O(nlogn))的时间和空间预处理出来的。
(sumlimits_{d=1}^{n}sumlimits_{k=1}^{n/d}mu(k)g(n/dk,k)g(m/dk,k))
设(T=dk)
(sumlimits_{T=1}^{n}sumlimits_{k|T}mu(k)g(n/T,k)g(m/T,k))
我们再设(sum(x,y,T)=sumlimits_{k|T}mu(k)g(x,k)g(y,k))。
但是猛然发现这个玩意微积分完后的时间和空间复杂度是(O(n^2))的，怎么用的起呀！！！！
但是这个时候有人灵光一闪，如果如果我们只处理(x,y<=B)的数字呢？（其实我们发现(k>=frac{n}{B})）
我们设(B=sqrt{n})。
那么就是(O(nsqrt{n}))的空间了！！！
用一个神奇的网站https://www.wolframalpha.com/，输入∫_sqrt{n}^{n}(frac{n}{i})^2di
但是预处理怎么办呢？
我们就爆枚吧QMQ。
爆枚时间(O(nsqrt{n}logn))，可以证明一个数字的约数个数的期望是(log)的。
然后我们就可以对于(k>=frac{n}{B})可以数论分块，而对于(k<=frac{n}{B})爆枚！！！
时间复杂度(O((T+n)sqrt{n}logn))
当然，这道题目很神奇，(B)取小一点反而能AC，估计又是常数问题。
#include<cstdio>
#include<cstring>
#include<vector>
#define N 51000
#define SN 210
using namespace std;
int mod=23333;
const int B=190;
vector<int>g[N];//g[k][n]
int pr[N],top;bool v[N];
int inv[N],u[N];
vector<int>r[SN][SN];
void first_do(int x)
{
inv[1]=1;u[1]=1;
for(int i=2;i<=x;i++)
{
if(!v[i])pr[++top]=i,inv[i]=i-1,u[i]=-1;
inv[i]%=mod;
for(int j=1;j<=top && i*pr[j]<=x;j++)
{
v[i*pr[j]]=1;
if(i%pr[j]==0)
{
inv[i*pr[j]]=inv[i]*pr[j];
u[i*pr[j]]=0;
break;
}
inv[i*pr[j]]=inv[i]*inv[pr[j]];u[i*pr[j]]=-u[i];
}
}
for(int i=0;i<=x;i++)g[0].push_back(0);
for(int i=1;i<=x;i++)
{
g[i].push_back(0);
for(int j=i;j<=x;j+=i)g[i].push_back((g[i-1][j/i]+inv[j])%mod);
}
for(int i=1;i<=B;i++)
{
for(int j=i;j<=B;j++)
{
int lim=x/j;r[i][j].resize(lim-B+1);
for(int k=1;k<=lim;k++)
{
//注意一个事情！o>B
for(int o=(B/k+1)*k;o<=lim;o+=k)r[i][j][o-B]+=u[k]*g[i][k]*g[j][k],r[i][j][o-B]=((r[i][j][o-B])%mod+mod)%mod;
}
}
}
for(int i=1;i<=B;i++)
{
for(int j=i;j<=B;j++)
{
int lim=x/j;
for(int k=B+1;k<=lim;k++)r[i][j][k-B]+=r[i][j][k-1-B],r[i][j][k-B]%=mod;
}
}
//预处理完毕
}
inline int mymin(int x,int y){return x<y?x:y;}
int main()
{
first_do(50000);
int T;scanf("%d",&T);
while(T--)
{
int n,m;scanf("%d%d",&n,&m);
n>m?(n^=m^=n^=m):0;
int ans=0;
int last;
for(int i=50000/B+1;i<=n;i=last+1)
{
last=mymin(n/(n/i),m/(m/i));
ans+=r[n/i][m/i][last-B]-r[n/i][m/i][i-1-B];ans%=mod;
}
//之前预处理的部分
int ed=mymin(50000/B,n);
for(int i=1;i<=ed;i++)
{
for(int j=i;j<=ed;j+=i)
{
ans+=u[i]*g[n/j][i]*g[m/j][i],ans%=mod;
}
}
printf("%d
",(ans+mod)%mod);
}
return 0;
}
8
这道题目其实是道SB题，只不过因为很多人（包括我）都是第一次遇到(prod)！！！（讲真我一开始打算用(log)换成加法，但是发现小数次幂我解决不了。。。）
(prodlimits_{d=1}^{n}f(d)^{sumlimits_{i=1}^{n/d}sumlimits_{i=1}^{m/d}[gcd(i,j)==1]})
(prodlimits_{d=1}^{n}f(d)^{sumlimits_{i=1}^{n/d}sumlimits_{i=1}^{m/d}sumlimits_{k|gcd(i,j)}mu(k)})
(prodlimits_{d=1}^{n}f(d)^{sumlimits_{k=1}^{n/d}mu(k)sumlimits_{i=1}^{n/dk}sumlimits_{i=1}^{m/dk}})
要不是多组数据我就用第6题方法艹过去了
设(T=dk)（我就不会下一步怎么化然后看了题解）
(prod_{T=1}^{n}prod_{d|T}f(d)^{mu(frac{T}{d})(n/T)(m/T)})
后面那个((n/T)(m/T))用快速幂加数论分块爆搞。
但是我们要先处理出(prod_{d|T}f(d)^{mu(frac{T}{d})})的前缀积。。。
所以我们(O(nlogn))预处理一下，然后就可以做了。
时间复杂度：(O(nlog{mod}+nlogn+Tsqrt{n}log{mod}))
实际上而言。((n/T)(m/T))可以(\%(mod-1))（费马小定理）
#include<cstdio>
#include<cstring>
#define N 1100000
using namespace std;
typedef long long LL;
LL mod=1e9+7;
LL u[N],g[N],f[N],fn[N]/*逆元*/;
int pr[N],top;bool v[N];
LL ksm(LL x,int k)
{
if(x==0)return 1;
LL ans=1;
while(k>1)
{
if(k&1)ans=(ans*x)%mod;
x=(x*x)%mod;k>>=1;
}
return (ans*x)%mod;
}
void pre_do(int x)
{
f[1]=fn[1]=1;for(int i=2;i<=x;i++)f[i]=(f[i-1]+f[i-2])%mod,fn[i]=ksm(f[i],mod-2);
u[1]=1;
for(int i=2;i<=x;i++)
{
if(!v[i])pr[++top]=i,u[i]=-1;
for(int j=1;j<=top && pr[j]<=x/i;j++)
{
v[pr[j]*i]=1;
if(i%pr[j]==0)
{
u[i*pr[j]]=0;
break;
}
u[i*pr[j]]=-u[i];
}
}
for(int i=1;i<=x;i++)g[i]=1;
for(int i=1;i<=x;i++)
{
if(!u[i])continue;
for(int j=i;j<=x;j+=i)
{
g[j]*=(u[i]==-1?fn[j/i]:f[j/i]),g[j]=(g[j]%mod+mod)%mod;
}
}
for(int i=2;i<=x;i++)g[i]=((g[i]*g[i-1])%mod+mod)%mod;
}//预处理
inline int mymin(int x,int y){return x<y?x:y;}
int main()
{
pre_do(1000000);
int T;scanf("%d",&T);
while(T--)
{
int n,m;scanf("%d%d",&n,&m);
n>m?(n^=m^=n^=m):0;
int last;
LL ans=1;
for(int i=1;i<=n;i=last+1)
{
last=mymin(n/(n/i),m/(m/i));
ans*=ksm((g[last]*ksm(g[i-1],mod-2))%mod/*求逆元*/,((LL)(n/i)*(m/i))%(mod-1));
ans%=mod;
}
printf("%lld
",(ans+mod)%mod);
}
return 0;
}
9
还有！！！
肾都炸了！！！！
博主你是不是故意的！！！
最后一道题目放一个毒瘤之神？
毒瘤之神：
方法1：
(φ(i,j)=φ(frac{i}{gcd(i,j)})*φ(j*gcd(i,j)))
但是明眼人一看就发现错了。
(frac{i}{gcd(i,j)})不一定与(gcd(i,j))互质。
方法2：
(φ(ij)=frac{φ(i)φ(j)gcd(i,j)}{φ(gcd(i,j))})
(sumlimits_{d=1}^{n}sumlimits_{i=1}^{n}sumlimits_{i=1}^{m}φ(i)*φ(j)frac{d}{φ(d)}[gcd(i,j)==d])
(sumlimits_{d=1}^{n}sumlimits_{i=1}^{n/d}sumlimits_{i=1}^{m/d}sumlimits_{k|gcd(i,j)}frac{d}{φ(d)}mu(k)φ(id)*φ(jd))
(sumlimits_{d=1}^{n}sumlimits_{k=1}^{n/d}sumlimits_{i=1}^{n/dk}sumlimits_{i=1}^{m/dk}mu(k)φ(ikd)*φ(jkd)frac{d}{φ(d)})
(sumlimits_{T=1}^{n}sumlimits_{k|T}sumlimits_{i=1}^{n/T}sumlimits_{i=1}^{m/T}mu(k)φ(iT)φ(jT)frac{frac{T}{k}}{φ(frac{T}{k})})
(sumlimits_{T=1}^{n}sumlimits_{k|T}mu(k)frac{frac{T}{k}}{φ(frac{T}{k})}sumlimits_{i=1}^{n/T}φ(iT)sumlimits_{i=1}^{m/T}φ(jT))
(G(n,k)=sumlimits_{i=1}^{n}φ(ik))
(F(T)=sumlimits_{k|T}mu(k)frac{frac{T}{k}}{φ(frac{T}{k})})
(S(x,y,z)=sumlimits_{T=1}^{x}sumlimits_{k|T}mu(k)frac{frac{T}{k}}{φ(frac{T}{k})}sumlimits_{i=1}^{y}φ(iT)sumlimits_{i=1}^{z}φ(jT))
(S(x,y,z)=S(x-1,y,z)+F(x)G(y,x)G(z,x))
所以我们只要处理(x<=n,1<=y,z<=B)的，对于(y,z>B)的我们暴力处理，所以时间复杂度是(O(T(B+frac{n}{B})+nB^2))，而且空间复杂度是(O(nB^2))。
当(B=35)时这个算法跑的很快。
#include<cstdio>
#include<cstring>
#include<vector>
#define N 110000
using namespace std;
const int B=35;
typedef long long LL;
LL mod=998244353;
inline LL ksm(LL x,int k)
{
if(k==0)return 1;
LL ans=1;
while(k>1)
{
if(k&1)ans=(ans*x)%mod;
x=(x*x)%mod;k>>=1;
}
return (ans*x)%mod;
}
vector<LL> g[N];
vector<int/*TM开到LL貌似会炸*/> s[40][40];//f[?][i][j]
int pr[N],top;bool v[N];
LL u[N],inv[N],ni[N],f[N];
void pre_do(int x)
{
u[1]=1;inv[1]=ni[1]=1;
for(int i=2;i<=x;i++)
{
g[i].resize(x/i+1);
if(!v[i])pr[++top]=i,u[i]=-1,inv[i]=i-1;
ni[i]=ksm(inv[i],mod-2);
for(int j=1;j<=top && pr[j]*i<=x;j++)
{
v[pr[j]*i]=1;
if(i%pr[j]==0)
{
u[i*pr[j]]=0;
inv[i*pr[j]]=inv[i]*pr[j];
break;
}
u[i*pr[j]]=-u[i];
inv[i*pr[j]]=inv[i]*inv[pr[j]];
}
}
g[0].resize(x+1);g[1].resize(x+1);
for(int i=1;i<=x;i++)
{
for(int j=i;j<=x;j+=i)g[j/i][i]=(g[j/i-1][i]+inv[j])%mod;
}
for(int i=1;i<=x;i++)
{
LL y=((LL)i*ni[i])%mod;
for(int j=i;j<=x;j+=i)f[j]+=y*u[j/i],f[j]%=mod;
}
for(int i=1;i<=x;i++)f[i]=(f[i]+mod)%mod;
for(int i=1;i<=B;i++)
{
for(int j=i;j<=B;j++)
{
int lim=x/j;s[i][j].resize(lim+1);
for(int k=1;k<=lim;k++)s[i][j][k]=((LL)s[i][j][k-1]+((f[k]*g[i][k])%mod)*g[j][k])%mod;
}
}
}
inline int mymax(int x,int y){return x>y?x:y;}
inline int mymin(int x,int y){return x<y?x:y;}
int main()
{
pre_do(100000);
int T;scanf("%d",&T);
while(T--)
{
int n,m;scanf("%d%d",&n,&m);
n>m?(n^=m^=n^=m):0;
int last;LL ans=0;
/*
m/i>B
1/i>B/m
i<m/B
*/
for(int i=mymax(m/B+1,1);i<=n;i=last+1)
{
last=mymin(n/(n/i),m/(m/i));
ans+=s[n/i][m/i][last]-s[n/i][m/i][i-1];
ans=(ans%mod+mod)%mod;
}
for(int i=1;i<=m/B;i++)
{
ans=(ans+((f[i]*g[n/i][i])%mod)*g[m/i][i])%mod;
}
printf("%lld
",ans);
}
return 0;
}
小结&总结证明
-
为了循环内(int)不爆，(a*b<=c)尽量写成(a<=c/b)，但是有人会担心(c/b)不是向下取整吗，那么我们就看看(a=ceil(c/b))满不满足吧。
证明：
(a*b=ceil(c/b)*b>=c)
咦不是有(=)吗，但是很抱歉，等于的情况仅当(b|c)，也就是(ceil(c/b)=floor(c/b))。
-
莫比乌斯反演最重要的几步：套(mu)，调换顺序，数论分块，套上另外的函数初始化，设数字替换（如(T=dk)的操作）。