Number Theory
202404
0 前言
离散数学是本书的重点,而整数又是离散数学的中心议题。数论 是讨论整数性质的重要数学分支,因此我们要来探索它。
——《具体数学》第四章
标有 *
的为拓展内容,或者说比较难的题目,但它们都非常有趣。
部分题目的代码是洛谷的提交记录,阅读这篇文章的大多数同学都该能看(若有需要可以找我要)。
没看懂的地方直接问,不排除我会存在手误。
1 基本概念
1.1 整除
两个整数相除,我们自然想去讨论商是否是一个整数。小学生常常期待这个,因为整数比分数写起来好看。
这引出了数论中的基础概念:整除。
如果
同样的,如果
1.2 整值函数
整数是离散数学的支柱,我们常常需要将分数或者任意的实数转换到整数。 ——《具体数学》第三章
我们首先来讨论 底(floor,最大整数)函数和 顶(ceiling,最小整数)函数,对于所有实数
俗称 下取整 和 上取整,底和顶也会在数论中有许多应用,在后续的部分内容中我们也会用到。
我们有时称
尝试在这里就证明一个非常有趣且有用的事实:如果
设
于是只要
画出图来容易发现这是显然的,这里不做证明。
我们发现其实
1.3 带余除法
当
告诉我们,可以将
这里的
我们约定
容易从定义证明这个法则,因为如果
且模数为零时该式显然也为真(两边都是
关于
。 。 。
前两条定义了
这些性质都是耳熟能详的,这里也不做证明。
1.4 同余关系
模运算是数论提供的一种重要工具,我们在上面把它当成二元运算使用,而在这里我们也可以把
由于
同余符号
乘法同样有效,只要处理的对象是整数:
证明直接作差:
这样一来,我们对方程所习惯做的大多数代数运算对同余式都可以运用,但并不是所有运算。除法运算显然不在其中。
以下部分可能会用到之后要讲的内容。
如果
然而在
它的证明是直接在式子两边乘上
进一步观察改变模的想法,我们可以得到另外一些式子:
反过来,如果我们知道对于两个小的模数有
因为如果
注意关注这个式子的特殊情形
2 欧几里得
有着我们熟悉的
2.1 欧几里得算法(gcd)
欧几里得算法,又称辗转相除法,常用递归得到最大公因数。
而这是容易简单证明的:
设
于是
由于
这样一来,
CF1806F2 GCD Master (hard version)
给定
和一个长度为 的序列 。 定义一次对一个长度为
的序列的操作为,选择序列中两个下标 ,删去 与 ,然后在序列末端加入 。 例如,对于
,一次操作可以选择下标 与 ,这样操作后,序列变成 。 给定
,求对序列 执行 次操作后得到序列中的数的和的最大值。
。
感谢🐨的供题。
首先我们考虑有相同的数如何处理?
发现我们对相同的数操作,相当于删去一个,进而在后面的讨论中我们发现其实相同的数并无意义,所以我们最后枚举删了几个相同的数就可以了,于是我们就只需要对不相同的数进行操作即可。
考虑每一次操作,删除两个数
那么这个问题就被转化成把数分成若干组,删去这一组并加上它们的
我们如何让分组最优?不难猜到把所有的分成一组。
假设现在我们存在两组
我们尝试合并两组,也就是把
由于
一定是成立的,所以合并两组一定更优,于是我们就得到了这一结论。
这样我们就可以取枚举
但在这道题中,依次枚举
假设我们现在加入一个没有选的最小的数
因为元素互不相同,所以
于是这样的调整一定是更优的。
这就意味这什么呢?
假设我们选了
也就是说我们只需要枚举选了多少个数和最大的数选的是什么就可以了。但这样还是会 T。
考虑把
双倍经验:CF1806F1 GCD Master (easy version)。
2.2 扩展欧几里得算法(exgcd)
扩展欧几里得算法,即 exgcd,旨在求不定方程
关于一个不定方程
如果
均为整数,则有整数 满足 。
证明是简单的,我们设
容易发现,如果
同时,根据 裴蜀定理 我们可以得到一个推论:
两个整数
互质,当且仅当 存在解 。
那如何求这组方程
我们考虑欧几里得算法的递归过程
当递归结束时
所以我们可以尝试先递归去求
那么
于是这就可以在欧几里得算法中解决了。
int exgcd(int a,int b,int &x,int &y){
if(b==0) return x=1,y=0,a;
int d=exgcd(b,a%b,y,x);
y-=a/b*x;return d;
}
那么我们如何从一组解推广到更多组解呢?(这样可以方便找到一组最小正整数解)
我们用 exgcd 得到了一组解
对于
所以
通过这样的调整,我们就可以找到想要的解了。
而回到最初,如果
P1082 [NOIP2012 提高组] 同余方程
求关于
的同余方程 的最小正整数解。
整理一下方程,即将同余拆开,我们可以得到
那么下面这个东西是直接可以用 exgcd 解决的,于是你就做完啦!
那如何找到最小整数解呢?
根据上面的转化,其实就是找到最小的
然后你就会求逆元了!
ABC340F S = 1
给你整数
和 ,它们至少满足 和 中的一个。 请找出一对满足以下所有条件的整数
。如果不存在这样的一对,请报告。
- 在
平面上,顶点为 的三角形的面积为 。
。
机翻的,还能看。一场非常近的 ABC,有些人在场上。
考虑如何计算这三个点组成的三角形面积呢?
用四边形面积直接减去三个角上的三角形,于是可以得到
而由于我们只考虑了一种象限的情况,所以
题目要求
2.3 类欧几里得算法
设
其中
这个式子和《具体数学》3.5 中的例三相当像啊,也就是说当
但是先不急,我们先把类欧几里得算法推完再来讨论这个问题。
如果
我们就把
接下来,用一些处理和式和处理底和顶的技巧
设
这就是一个递归的式子了,交换了
这也就是欧几里得算法的辗转相除过程,于是类欧几里得名字就是这样来的。
接下来我们对于类欧几里得算法进行一个推广,尝试去求另外两个和式:
首先来推导
接下来考虑
设
接着,我们推导
考虑
考虑
那么
还是设
三个函数一起进行递归,于是你就可以通过 P5170 【模板】类欧几里得算法 啦!代码。
const ll mod=998244353,i2=499122177,i6=166374059;
ll T,a,b,c,n;
struct node{
node () {f=g=h=0;}
ll f,g,h;
}ans;
node calc(ll a,ll b,ll c,ll n){
ll ac=a/c,bc=b/c,m=(a*n+b)/c,in=n*(n+1)%mod*(2*n+1)%mod*i6%mod;
node res;
if(a==0){
res.f=bc*(n+1)%mod;
res.g=bc*n%mod*(n+1)%mod*i2%mod;
res.h=bc*bc%mod*(n+1)%mod;
return res;
}
if(a>=c||b>=c){
res.f=n*(n+1)%mod*i2%mod*ac%mod+(n+1)*bc%mod;res.f%=mod;
res.g=in*ac%mod+n*(n+1)%mod*i2%mod*bc%mod;res.g%=mod;
res.h=ac*ac%mod*in%mod+bc*bc%mod*(n+1)%mod+ac*bc%mod*n%mod*(n+1)%mod;res.h%=mod;
node e=calc(a%c,b%c,c,n);
res.h+=e.h+2*bc%mod*e.f%mod+2*ac%mod*e.g%mod;
res.g+=e.g;res.f+=e.f;
res.h%=mod;res.f%=mod;res.g%=mod;
return res;
}
node e=calc(c,c-b-1,a,m-1);
res.f=(n*m%mod-e.f+mod)%mod;
res.g=(m*n%mod*(n+1)%mod-e.h-e.f+mod+mod)%mod*i2%mod;
res.h=(n*m%mod*(m+1)%mod-2*e.g%mod-2*e.f%mod-res.f%mod+3*mod)%mod;
return res;
}
P5179 Fraction
给你四个正整数
,求一个最简分数 满足 。 若有多组解,输出
最小的一组,若仍有多组解,输出 最小的一组。
。
一道用到类欧类似思想的题目。
什么时候我们可以不加思考地直接算出这个目标分数?
发现当
举一个例子:
对其取倒数,可以得到
把两边
我们把这两边的分数变成更小的真分数,继续循环:
从而一步一步带回即可。代码。
这道题中我们把假分数转成真分数的操作和类欧类似。
ABC283Ex Popcount Sum
组询问,每组求 内, 余 的数在二进制下 的个数的和。
。
对 (x>>i)&1
那么转化成数学表达式就是
而由于我们只统计
这个东西就是容易用类欧几里得算法在
*当 时的封闭形式
我们来尝试求
的封闭形式,也就是
前置知识:来自《具体数学》3.25 和 3.26。
如何证明这个东西?
你可以看成把
个物品分成 组,那么一定有一些组是 个物品,另外一些就是 个。 而从小到大排好序后,第
组刚好有 个元素。
同时,我们还可以对这个东西进行推广,用
替换 ,就可以得到一个对所有实数 都成立的恒等式:
首先可以考虑打表,或许你可以发现一些东西,具体可以去看《具体数学》3.5。
而推导时我们可以用刚才推导类欧几里得类似的方法:对
把它拆成整数部分和小数部分就是
考虑把三项分开来求和,对于第一项和第三项,我们都涉及到了
那么证明处理这个神秘东西呢?
容易发现这是简单的,设
于是它所得到的就是把按照某种次序排列的数
所以这个式子的第一部分的和就是
最后一步就用到了上面前置知识中的东西。
而对于第三部分,我们发现它其实就是一个等差数列重复了
中间的第二部分,也是非常简单的,同样是一个等差数列,它的和就是
于是我们要求的封闭形式就找出来了
如果对封闭形式稍作处理,就可以让它关于
就得到这样的互反律
类欧几里得算法,完。
3 同余关系
在前置知识中提到了同余关系的相关法则,那么在这里我们希望对它的理解更为深刻一些。
3.1 费马小定理
费马小定理:
我们知道,
这个式子是在
由于
同时,费马小定理还有一个更通用的表达方式
3.2 逆元
如果一个线性同余方程
容易发现,这个东西可以用上面我们讨论过的 exgcd 直接求解,同时当
这样两种方法都可以做到单次求逆元
而根据前面我们对 exgcd 是否有解的讨论,容易发现
可是有些时候每次都单独求逆元并不能满足我们所需要的时间复杂度,所以这里就有了另外两种求多个数逆元的方法:
-
线性求逆元:
首先我们知道
,我们尝试线性推出 。那么令
,于是 有 。我们将式子两边同时乘上
可以得到:这样就可以做到线性递推求逆元了。时间复杂度
,就可以通过 P3811 【模板】模意义下的乘法逆元 啦! -
线性求
个数 的逆元:利用前缀积求逆元设
是 的前缀积,那么如果我们知道了 的逆元,于是可以从 倒推回去: 。这样每个数的逆元就是
,时间复杂度 ,这种方法在求阶乘的逆元中常常用到。
3.3 二次探测定理
二次探测定理:
如果
我们先不急着讨论这个东西,先来讨论
首先我们应该考虑
所以
这就意味着当
现在,
每个素数都是独立于其他素数的,对于
于是如果
在这个推导过程中,我们已经证明了 二次探测定理 了!
3.4 威尔逊定理
接下来引出一个逆元的应用——威尔逊定理,他会在之后有关素数的判定中有所应用。
这个定理的一半是平凡的:如果
具体来说,在 1.4 我们引入了一个改变模数的式子,也就是如果
威尔逊定理的另一半说的是,
因为
这里
那哪些数的逆元就是它本身呢?
容易发现其实就是问
当
而
3.5 有趣的遗留问题
在上面推类欧几里得算法的一个拓展时,我们希望证明
按照某种次序恰好组成了
的
或许上面已经给出了一些不太严谨的简要证明,我们在这里希望把它写得更完整一些。
首先证明得到前面
从而当
那么现在我们必须指出,那
我们记
而这种情形就是
根据 鸽巢原理,我们就可以轻松证明出这
它可以被理解成在一个环上面每次跳
4 线性同余方程组
求解形如
的 线性同余方程组,在 OI 中有着广泛的应用。
4.1 中国剩余定理(CRT)
中国剩余定理全称 Chinese Remainder Theorem,简称 CRT。它用于求解 模数两两互质 的线性同余方程组,即对于任意
首先我们来看一个例子:
在《孙子算经》中有这样一个问题:“今有物不知其数,三三数之剩二(除以
余 ),五五数之剩三(除以 余 ),七七数之剩二(除以 余 ),问物几何?
CRT 的步骤是这样的:
-
找出三个数:
从
和 的公倍数中找出被 除余 的最小数 ;从
和 的公倍数中找出被 除余 的最小数 ;从
和 的公倍数中找出被 除余 的最小数 。 -
用
( 为最终除以 的余数),用 ( 为最终除以 的余数),用 。 -
把以上三个数相加得到
,再用它除以 的最小公倍数 ,得到余数 即为最小满足条件的数。
它为什么对的呢?
我们假设
从
同理我们可以得到:
除以 余 ,且是 和 的公倍数。 除以 余 ,且是 和 的公倍数。 除以 余 ,且是 和 的公倍数。
所以 CRT 用数学表示出来就是:
设
因为
#include<bits/stdc++.h>
using namespace std;
#define ll long long
const int N=15;
int n;
ll p,a[N],b[N],m,x,y,ans=0;
void exgcd(ll a,ll b ,ll &x,ll &y){
if(b==0){x=1,y=0;return;}
exgcd(b,a%b,y,x);y-=a/b*x;
}
int main(){
scanf("%d",&n);m=1ll;
for(int i=1;i<=n;i++){
scanf("%lld%lld",&a[i],&b[i]);
m*=a[i];
}
for(int i=1;i<=n;i++){
p=m/a[i];
exgcd(p,a[i],x,y);
ans+=b[i]*p*(x<0?x+a[i]:x);
}
printf("%lld\n",ans%m);
return 0;
}
4.2 扩展中国剩余定理(exCRT)
如果就是单纯的 线性同余方程组 怎么解呢?也就是
这就是 exCRT 所解决的问题,它能解决的东西更为普遍,故使用范围也会更广。
先考虑只有两个方程的简单情况。
容易发现
那么这就是一个可以用 exgcd 求解的方程了。
如果
否则,我们可以得到新的方程,也就是合并了两个方程:
推广到
我们就成功求解了 线性同余方程组!
P4777 【模板】扩展中国剩余定理(EXCRT) 的代码。
ll mul(ll a,ll b,ll m){
ll res=0ll;
while(b){if(b&1) res=(res+a)%m;a=(a+a)%m;b>>=1;}
return res;
}
ll exgcd(ll a,ll b, ll &x,ll &y){
if(b==0){x=1,y=0;return a;}
ll d=exgcd(b,a%b,y,x);y-=(a/b)*x;return d;
}
ll wrk(){
ll ans=a[1],M=b[1];
for(int i=2;i<=n;i++){
ll A=M,B=b[i],C=(a[i]-ans%B+B)%B;
ll gcd=exgcd(A,B,x,y),bg=B/gcd;
if(C%gcd!=0) return -1;
x=mul(x,C/gcd,bg);
ans+=x*M;M*=bg;
ans=(ans%M+M)%M;
}
return (ans%M+M)%M;
}
5 素数
接下来我们来讨论素数。
5.1 素数的例子
如果一个正整数
而其他那些还有非平凡因子的数都称为 合数,每一个大于
任何正整数
而 算数基本定理 表示:仅有一种方法将
关于素数的其他定理:
-
第
个素数 大约是 的自然对数的 倍: -
对于不超过
的素数个数 :
证明过于复杂,这里不给出证明,大家可以自行了解。
还有一种与素数相关的数:梅森数。
形如
的数,其中
如果
但是当
5.2 互素关系
当
我们写成
而一个分数是最简分数,当且仅当
那我们有没有什么好的方法来构造满足
答案是肯定的,他被成为 Stern-Brocot 树,其思想是从两个分数
- 在两个相邻的分数
和 之间插入 。
新的分数就是
接下来就是
于是我们把这个分数阵列可以看成是一棵无线的二叉树构造,就是这样
每一个分数都是
为什么这种构造是对的呢?为什么每次的
因为我们发现,在这个构造中任意一个阶段的相邻分数都满足
首先这个关系在
括号展开就一样了,而这也让我们发现
一个中位分数不一定在原先两个的中间,但它的确位于它们两个中间的某个地方,于是我们就不可能得到相同的分数。
但是我们现在仍然存在一个问题——会遗漏么?
这个问题时简单的,由于它是无限的,所以我们会不断用上面的单调性去逼近它,这样就不会漏。
这里我们再引入一个概念,阶为
容易发现它是 S-B 树 的子树,好像也没什么好讲的。
事实上,我们可以把 S-B 树 看成一个表示有理数的 数系,因为每一个正的最简分数都恰好出现一次。
我们用字母
我们记当前的字符串为
例如
这个东西和矩阵乘法的构造是比较像的,也就是我们尝试去构造矩阵
那么我们知道
然后我们就发现其实
于是我们就可以通过一个
还有一种方法就是改变
因为
类似的性质对
这样的定位就可以不用矩阵乘法。
关于互素的讨论到此结束,更多可以参见《具体数学》!
P8058 [BalkanOI2003] Farey 序列
把所有分子和分母都
的 最简真分数 从小到大排成一行,形成的序列称为 Farey 序列。 求出
所对应的 Farey 序列中第 小的数。
符合条件分数的个数。
sto smb orz
给出一个 smb 的亚线性做法。
有了上面的一些简单铺垫,我们不难想到在 S-B 树 上面进行二分。
那么对于一个最简分数
容易列出式子
对这个式子进行莫反,可以得到
我们发现前面是可以用杜教筛处理的(其实这道题中不用),而后面这一块是类欧几里得的板子,于是这个值我们是可以在
而我们如何去二分到这个值呢?
不难想到直接在 S-B 树 上面走,每次枚举走到左子树还是走到右子树,但是在
我们考虑把往一个方向走的连续段缩起来,也就是只枚举在哪些地方拐(eg:原先一直往左子树走,现在变成往右子树走),设
那么这样我们只会拐
总时间复杂度
ABC333G Nearest Fraction
给定一个小于
的正实数 和一个正整数 。 要求在满足 和 的前提下,找到使 最小的二元组 。 如果存在多个这样的二元组 ,输出 值最小的那个。 数据范围:
, 且最多包含 18 位有效数字。
首先
对于当前的区间
这样我们就可以不断减小这个区间,从而不断逼近
每一次搜索都用两个答案
用上一道题类似的倍增优化可以做到
双倍经验:P1298 最接近的分数
SP26073 DIVCNT1 - Counting Divisors
表示 的约数个数。 求
多测,
。
比较难的应用。
首先我们可以对式子进行转化
那么问题就转化成求
我们尝试把这个反比例函数画出图来,那么我们要求的其实就是下图中红色部分中的点数(from x义x):
不妨将其称作 R 区域。
我们希望用一些向量去不断拟合这个反比例函数,于是就想到了这个有理数系:S-B 树。
于是我们需要一个单调栈,里面维护的向量斜率单调递减,进行这样的操作:
- 每次从栈顶取出一个元素,不断加到当前点
上面,直到它差一步就走进 R 区域。
由于我们保证了斜率递减,且向量的斜率表示出来已经是最简分数,所以每一步都是对答案有贡献的,也就是图中棕色区域中的点(from x义x):
容易发现这样的点数一定不会多算,可以做到不重不漏。(在斜线下的点一定在反比例函数中,并且不会多)
如何计算点数?(我们假设每次上都要计算最下面的一行,从而不能计算上面的一行)
画出图来就是这个样子:
把它分成蓝色和橙色的两个部分,而橙色部分一定是上下对称的,所以假设当前向量为
- 不断弹出栈顶,使得攻栈顶刚好走不进 R 区域,称它为
。设最近一次弹的是 (有可能是上一步的栈顶向量)。
现在就是,加上
- 开始二分,用 S-B 树 得到中间向量
,我们称之为 。
-
如果
走一步不会进 R,那么 进栈并把 继续二分。 -
否则
- 如果
,那么继续二分肯定都走不出 R 了,直接停止二分,就是这个样子(from x义x):
- 如果
可以发现如果
- 否则 $l=mid$ 继续二分。
当
这样的时间复杂度是
注意到其他的一些东西:
-
最开始的点我们选择
,因为这样才能做到不重不漏(细节会少很多)。 -
初始栈中的向量是
和 。
这样我们就做完了。代码。
5.3 素数筛法
最重要的知识点,几乎所有题目都需要筛出
首先我们存在一个非常显然的做法——埃氏筛法。
埃氏筛用所有已经筛出的素数倍数标记合数:从
埃氏筛的精髓在于复杂度的证明:不超过
也就是
这说明埃氏筛的时间复杂度是
这里是来自 djq 的证明:
因为每个数只会被其素因子筛到,所以
根据
根据
两边同时取对数
因此
埃氏筛已经非常接近与线性了,但是它也还不是最优的,我们存在一个 线性筛素数 的做法。
它的思想和埃氏筛类似,埃氏筛的优秀之处在于仅用质数的倍数筛去合数,但合数会被多个质因子筛到。让每个合数仅被晒一次就能做到线性。
考虑用每一个合数 最小质因数 筛它:从
综上,有如下步骤:
-
从小到大遍历当前筛出的所有素数
,将 标记为合数。 -
若
,退出循环,因为 ,所以 的最小质因子为 不是 ,再筛就重复了。
时间复杂度线性,于是就可以通过 P3383 【模板】线性筛素数。
for(int i=2;i<N;i++){
if(!vis[i]) prm[++cnt]=i;
for(int j=1;j<=cnt&&prm[j]*i<N;j++){
vis[prm[j]*i]=true;
if(i%prm[j]==0) break;
}
}
5.4 素数判定(Miller Rabin 素性测试)
如何判断一个数是否是素数?
小素数的判定是容易的:直接试除法即可,时间复杂度
至于大素数判定……
费马曾提出一种方法:费马素性测试。
它基于费马小定理:设
为了测试
-
如果
不成立,那么 肯定不是素数。 -
如果
成立,那么 很大概率是素数,尝试的 越多, 是素数的概率就越大。
但是从第
对于某个
特别地,有一些合数,不管选什么
不过这种数很少,
由于上面 费马素性测试 的不完美,我们就引入了 Miller-Rabin 素性测试。
把费马素性测试改进一下,它是已知最快的随机素数测试算法。
原理也非常简单,就是费马素性测试排除 Carmicheael 数。
因为并且 Carmicheael 数在进行次数较多的二次探测定理逆否命题判定后,都会被判定出来。
如果
根据费马测试,如果
而由于
把
于是有:
所以计算出
容易发现,对于多个随机的基值
当
考虑做了
ll mul(ll a,ll b,ll m){ll z=(long double)a/m*b,res=(a*b-z*m)%m;return (res+m)%m;}
ll add(ll a,ll b,ll n){return (a+b>=n)?a+b-n:a+b;}
ll qpow(ll a,ll b,ll m){
ll res=1ll;
while(b){
if(b&1) res=mul(res,a,m);
a=mul(a,a,m);
b>>=1;
}
return res;
}
const ll prm[12]={2,3,5,7,11,13,17,19,23,29,31};
bool miller_rabin(ll n){
if(n<2) return false;
for(int i=0;i<=10;i++){
if(n==prm[i]) return true;
if(n%prm[i]==0) return false;
}
ll u=n-1,t=0;
while(!(u&1)) u>>=1,t++;
for(int i=0;i<=10;i++){
ll a=prm[i],x1=qpow(a,u,n),x2;
for(int j=1;j<=t;j++){
x2=mul(x1,x1,n);
if(x2==1&&x1!=1&&x1!=n-1) return false;//二次探测定理
x1=x2;
}
if(x1!=1) return false;
}
return true;
}
5.5 分解质因数(Pollard-rho)
判断完了素数,我们来尝试分解质因数。
众所周知,可以用试除法做到
所以我们考虑一种更高消的方法——Pollard-rho 启发式方法 分解质因数。
事实上,我们还有一种奇妙的方法来寻找一个合数的因子(from TQX):
scanf("%d",&n);
vector<int> fac;
for(int i=1;i<=???;++i){
int a=rand()%(n-1)+1;
if(n%a==0) fac.push_back(a);
}
这个代码试图通过随即一些数并判断它们是否是
这个算法非常的劣质,但我们的 Pollard_rho 正是基于这个算法而来的,它的核心思想就是随机。
而我们知道 生日悖论:
如果现在有
而当
这个悖论告诉我们:组合随机采样比单个个体满足条件要高得多,这样可以提高正确率。
于是,Pollard 就提出了一个随机数的生成方法:
其中
这个函数的随机性非常好,但是进行了数次生成之后就会出现一个环,在坐标系中表示出来就变成了这个样子:
发现这个图像很像
那么考虑如何判环?
可以用 Floyd 的方法:
假设我们有
于是我们可以用这样随机出来的数相邻两个的差与
但是 gcd 的时间复杂度时
在计算时,我们将每次产生的
但是,如果某一时刻累积下来的样本的乘积为
每次计算的阈值可以倍增,并且加上一个上限,使用
至此,你就可以通过 P4718 【模板】Pollard-Rho 了。代码。
ll mul(ll a,ll b,ll m){ll z=(long double)a/m*b,res=(a*b-z*m)%m;return (res+m)%m;}
ll add(ll a,ll b,ll n){return (a+b>=n)?a+b-n:a+b;}
ll qpow(ll a,ll b,ll m){
ll res=1ll;
while(b){
if(b&1) res=mul(res,a,m);
a=mul(a,a,m);
b>>=1;
}
return res;
}
const ll prm[12]={2,3,5,7,11,13,17,19,23,29,31};
bool miller_rabin(ll n){
if(n<2) return false;
for(int i=0;i<=10;i++){
if(n==prm[i]) return true;
if(n%prm[i]==0) return false;
}
ll u=n-1,t=0;
while(!(u&1)) u>>=1,t++;
for(int i=0;i<=10;i++){
ll a=prm[i],x1=qpow(a,u,n),x2;
for(int j=1;j<=t;j++){
x2=mul(x1,x1,n);
if(x2==1&&x1!=1&&x1!=n-1) return false;
x1=x2;
}
if(x1!=1) return false;
}
return true;
}
ll gcd(ll a,ll b){return (b==0)?a:gcd(b,a%b);}
ll f(ll a,ll c,ll n){return add(mul(a,a,n),c,n);}
ll rad(ll n){return 1ll*rand()*rand()%n*rand()%n+1ll;}
ll pollard_rho(ll n){
if(n==4) return 2;
ll x=rad(n-1),c=rad(n-1),y=x;
x=f(x,c,n);y=f(x,c,n);
for(int lim=1;x!=y;lim=min(128,lim<<1)){
ll cur=1ll,nxt=0;
for(int i=1;i<lim;i++,cur=nxt){
nxt=mul(cur,abs(add(x,n-y,n)),n);
if(!nxt) break;//x=y
x=f(x,c,n);y=f(f(y,c,n),c,n);
}
ll d=gcd(cur,n);
if(d!=1) return d;
}
return n;
}
void sol(ll n){
ll d=pollard_rho(n);
while(d==n) d=pollard_rho(n);
if(miller_rabin(d)) ans=max(d,ans);
else sol(d);
if(miller_rabin(n/d)) ans=max(ans,n/d);
else sol(n/d);
}
void wrk(){
n=read();
if(miller_rabin(n)){puts("Prime");return;}
ans=0;sol(n);print(ans);
}
P4358 [CQOI2016] 密钥破解
一种非对称加密算法的密钥生成过程如下:
1.任选两个不同的质数
2.计算
, 3.选取小于
,且与 互质的整数 4.计算整数
,使得 5.二元组
称为公钥,二元组 称为私钥。 当需要加密消息
时,(假设 是一个小于 的整数,因为任何格式的消息都可转为整数表示),使用公钥 ,按照 运算,可得到密文
对密文
解密时,用私钥 ,按照 运算,可得到原文
。算法正确性证明省略。 由于用公钥加密的密文仅能用对应的私钥解密,而不能用公钥解密,因此称为非对称加密算法。通常情况下,公钥由消息的接收方公开,而私钥由消息的接收方自己持有。这样任何发送消息的人都可以用公钥对消息加密,而只有消息的接收方自己能够解密消息。
现在,你的任务是寻找一种可行的方法来破解这种加密算法,即根据公钥破解出私钥,并据此解密密文。
。
题意写得非常长,但是并不难。
发现只要我们找到
P4714 「数学」约数个数和
给你一个正整数
,请计算 的 所有约数的 约数个数和。 答案可能很大,请输出对
取模的结果。
。
首先,对于
而当
这相当于一个狄利克雷卷积,而根据狄利克雷卷积卷积的性质,我们可以知道
于是现在问题就变成了如何求
发现这是简单的,等价于组合数学中在
根据对称公式,相当于
而
至于分解质因数,交给 Pollard_Rho 就可以了。代码。
6 离散对数问题
离散对数问题就是解决模
这样的
当
6.1 大步小步算法(BSGS)
大步小步算法全称 Baby Step Giant Step,简称 BSGS,适用于
整体思路用到了 根号平衡 的思想,设块长为
于是这个方程就可以被表示成
只要我们直接暴力枚举
根据欧拉定理(稍后会讲),有解就一定在
那么设块长为
于是你就可以通过 P3846 [TJOI2007] 可爱的质数/【模板】BSGS 啦!代码。
int BSGS(int a,int b,int H){
int A=1,B=sqrt(H)+1;a%=H,b%=H;
if(H==1||b==1) return 0;
if(!a) return !b?1:-1;
gp_hash_table<int,int> mp;
for(int i=1;i<=B;i++) mp[1ll*(A=1ll*A*a%H)*b%H]=i;
for(int i=1,AA=A;i<=B;i++,AA=1ll*AA*A%H)
if(mp.find(AA)!=mp.end()) return i*B-mp[AA];
return -1;
}
容易发现你从小到大枚举时就可以得到最小的非负整数解,而关于
在这里还有一些板题:
-
P2485 [SDOI2011] 计算器 :综合一点的板题。代码。
P3306 [SDOI2013] 随机数生成器
最近小 W 准备读一本新书,这本书一共有
页,页码范围为 。 小 W 很忙,所以每天只能读一页书。为了使事情有趣一些,他打算使用 NOI2012 上学习的线性同余法生成一个序列,来决定每天具体读哪一页。
我们用
来表示通过这种方法生成出来的第 个数,也即小 W 第 天会读哪一页。这个方法需要设置 个参数 ,满足 ,且 都是整数。按照下面的公式生成出来一系列的整数: 其中
表示取余操作。 但是这种方法可能导致某两天读的页码一样。
小 W 要读这本书的第
页,所以他想知道最早在哪一天能读到第 页,或者指出他永远不会读到第 页。
是质数
考虑先全部左移一维,那么得到的数其实就是
于是我们就可以将
这就可以用 BSGS 直接完成啦——注意特判
总时间复杂度
P4884 多少个 1?
给定整数
和质数 ,求最小的正整数 ,使得 ( 个 ) 。 说人话:就是
。
, ,保证 是质数。
a 如果我们直接把
我们考虑将这个数表示成一个数,也就是它乘
那么此时这个方程就变成了
于是这个东西就可以直接用 BSGS 完成了。
注意可能在乘的过程中爆 long long。代码。
CF1106F Lunar New Year and a Recursive Sequence
有一串
个数的数列,给你 。当 时: 已知
,问 可能是多少。
, 。
会用到原根相关知识,若不会原根建议先看后面原根部分。
有关于
而这种需要让我们想到了 原根,众所周知,
把指数提下来,同时模数变成了
那么只要我们能求得
发现这个东西刚好是矩阵乘法的形式,也就是我们可以构造
从而得到
于是当我们知道
而由于前面都是
其中
由于
现在问题就转化成如何求出
于是这道题就做完了,实现时注意时模
6.2 扩展大步小步算法(exBSGS)
接下来我们来讨论更为普遍的问题
从已知到未知: 由于
我们运用同余方程的运算法则,我们把两边同时除以
由于
然而这时的
而注意到每除依次,
容易发现这样的时间复杂度是
这样就可以通过 P4195 【模板】扩展 BSGS/exBSGS 和 SP3105 MOD - Power Modulo Inverted 了。代码。
int gcd(int a,int b){return b==0?a:gcd(b,a%b);}
void exgcd(int a,int b,int &x,int &y){
if(b==0) return x=1,y=0,void();
exgcd(b,a%b,y,x),y-=x*(a/b);
}
int inv(int a,int p){exgcd(a,p,x,y);return (x%p+p)%p;}
int BSGS(int a,int b,int H){
int A=1,B=sqrt(H)+1;a%=H,b%=H;
if(H==1||b==1) return 0;
if(!a) return !b?1:-1;
gp_hash_table<int,int> mp;
for(int i=1;i<=B;i++) mp[1ll*(A=1ll*A*a%H)*b%H]=i;
for(int i=1,AA=A;i<=B;i++,AA=1ll*AA*A%H)
if(mp.find(AA)!=mp.end()) return i*B-mp[AA];
return -1;
}
int exBSGS(int a,int b,int H){
int d=gcd(a,H),cnt=0;a%=H,b%=H;
while(d>1){
if(b==1) return cnt;
if(b%d) return -1;
H/=d,b=1ll*(b/d)*inv(a/d,H)%H;
d=gcd(H,a%=H),++cnt;
}
int res=BSGS(a,b,H);
return ~res?res+cnt:-1;
}
P5345 【XR-1】快乐肥宅
给定
组 ,求关于 的方程组 (定义
)在 范围内的最小整数解,或判断在这个范围内无解。
。
难绷。
容易发现,我们的解题过程分成两部分:求解每一个方程的解 和 求线性同余方程组的解。
首先考虑第一部分,如何求解每一个方程的解呢?
根据 exBSGS 的求解过程可以知道,解的情况一定分成了三种:无解、唯一解和无穷解。
画出图来,解一定是一个
其中
而进入环上就变成了无穷解,而这个解是形如线性同余方程;
反之,没有出现过的就是无解。
尾巴的长度是可以直接求
反之对于只有一个解的情况,我们直接记它为答案,最后特判一下就可以了。
接下来考虑第二部分,求
由于不保证模数互质,所以这个需要用 exCRT 解决,而
考虑 exCRT 的过程,中间每一次我们都是 ans+=x*M
,那么只要这次的
所以这个时候我们只需要判断当前的解对于后面的方程是否都成立即可。
这样就做完了,时间复杂度为
7 数论函数基础
数论函数是数论中相当重要的一环,我们先来将一些基本的函数。
7.1 相关定义
-
数论函数:定义域为正整数的函数称为 数论函数。因其在所有正整数处均有定义,故可视作数列。OI 中常见的数论函数的陪域(即可能的取值范围)为整数。
-
加性函数:若对于任意
且 均有 ,则称 为 加性函数。注意区分代数中的加性函数。 -
积性函数:若对于任意
且 均有 ,则称 为 积性函数。容易发现 是必要条件。 -
完全积性函数:若对于任意
均有 ,则称 为 完全积性函数。完全积性函数一定是积性函数。 -
数论函数的 加法:对于数论函数
, 表示 和 对应位置相加,即 。 -
数论函数的 数乘:对于数
和数论函数 , 表示 的各个位置乘 ,即 。一般简记作 。 -
数论函数的 点乘:对于数论函数
, 表示 和 对应位置相乘,即 。为与狄利克雷卷积符号 作区分,点乘符号通常不省略。
7.2 常见函数
设
-
单位函数:
。完全积性函数。 -
常数函数:
。完全积性函数。 -
恒等函数:
。 记作 。完全积性函数。 -
除数函数:
。 表示 的约数个数,记作 或 。 表示 的约约数和,记作 。 有计算式
根据乘法分配律,
等比数列求和后即可。
-
欧拉函数:
,表示 以内与 互质的数的个数。后面我们会有所介绍。 -
本质不同质因子个数函数:
,表示 的本质不同质因子个数。 -
莫比乌斯函数:
以上所有函数除
7.3 线性筛积性函数
在前面的素数部分,我们分析了线性筛素数,它提供了在线性时间内筛出具有特殊性质的积性函数在
只要积性函数
根据积性函数的性质,只要预先求出
下面放一个框架,from Alex_Wei:
for(int i = 2; i < N; i++) {
if(!vis[i]) pr[++pcnt] = i, f[i] = ..., low[i] = i; // 单独算 f(p)
for(int j = 1; j <= pcnt && i * pr[j] < N; j++) {
vis[i * pr[j]] = 1;
if(i % pr[j] == 0) { // i 与 p 不互质
low[i * pr[j]] = low[i] * pr[j];
if(i == low[i]) f[i * pr[j]] = ...; // i = p ^ k,单独算 f(p ^ {k + 1})
else f[i * pr[j]] = f[i / low[i]] * f[low[i * pr[j]]];
break;
}
low[i * pr[j]] = pr[j];
f[i * pr[j]] = f[i] * f[pr[j]]; // i 与 p 互质,f(ip) = f(i)f(p)
}
}
8 狄利克雷卷积
钢筋混凝土。
8.1 定义与性质
两个数论函数
上面简记成
狄利克雷卷积还有一些有趣的性质。
交换律,分配律,结合律。
这三个东西都可以用定义式稍微推一下就证明出来了。
单位元:
。
逆元:若函数
满足 ,则称 互为逆元。
积性函数总是有且仅有一个逆元。
可逆当且仅当 。
下面写出证明和逆元的计算式:
设
这同时也说明了逆元的 唯一性。
的充要条件是 ,其中 。
直接乘上
积性函数的狄利克雷卷积也是积性函数。
设有两个积性函数
第二行中,由于
积性函数的逆元也是积性函数。
设有积性函数
使用归纳法。当
综合后面的两个性质,两个积性函数的积与商都是积性函数,但和与差不是。
8.2 线性筛狄利克雷卷积
给出
我们直接用定义式计算,把枚举因数变成枚举倍数,根据调和级数,时间复杂度为
如果
对于第一项和第三项,我们都可以在线性筛中总代价
如果
当
因此,
综上,使用线性筛求出两个在质数幂处取值已知的积性函数的狄利克雷卷积在
这样我们就得到了积性函数线性筛更弱的条件:可以
P6222 「P6156 简单题」加强版
给出
,多次给出 ,求
。
会用到稍后要讲的莫比乌斯反演。
枚举 gcd,再进行莫比乌斯反演
我们需要线性筛预处理出
这就需要用到我们上面讨论的线性筛狄利克雷卷积。代码。
简单版本:P6156 简单题。
8.3 狄利克雷前缀和
对于任意数论函数
对每个
首先容易想到一种暴力做法,直接枚举每一个数的倍数,这样就可以做到
但是并不是最优的,于是我们有了下面的算法。
将每个数
根据高位前缀和的求法,枚举每一维并将所有下标关于该维做前缀和,可以得到狄利克雷前缀和的实现方法:
初始令
最终得到的
于是你就可以通过 P5495 【模板】Dirichlet 前缀和 了。代码。
#include <bits/stdc++.h>
using namespace std;
#define uint unsigned int
const int N=2e7+5;
uint seed,ans,a[N];
int n,p[N>>3],cnt=0;
bool vis[N];
uint rd(){seed^=seed<<13,seed^=seed>>17,seed^=seed<<5;return seed;}
void init(){
for(int i=2;i<N;i++){
if(!vis[i]) p[++cnt]=i;
for(int j=1;j<=cnt&&i*p[j]<N;j++){
vis[i*p[j]]=true;
if(i%p[j]==0) break;
}
}
}
int main(){
cin>>n>>seed;init();
for(int i=1;i<=n;i++) a[i]=rd();
for(int i=1;i<=cnt;i++) for(int j=1;p[i]*j<=n;j++) a[p[i]*j]+=a[j];
for(int i=1;i<=n;i++) ans^=a[i];
cout<<ans;return 0;
}
9 整除分块
整除分块又称为数论分块,因其解决的问题与整除密切相关而得名。数论分块是用于求解形如
的合适,前提为
9.1 介绍
如果
,不同的 至多只有 个。
证明是简单的:
-
当
时, 只有 个。 -
当
时, ,也只有 个。
根据这条结论,我们枚举
只要我们能对
那么问题就转化成求使得
转化一下,发现
和
两边把
那么我们不重不漏地枚举所有整除值时,可以对于每一个
具体地,我们令当前枚举到的
将
这样每一次的整除值都可以被算一次,时间复杂度
-
当
的上界不为 时,要特殊处理:若 ,每次要把 和 取最小值;反之要特殊处理 的情况。 -
写的时候不要写成了
了!(有些人因为这个常常调很久)
9.2 拓展
拓展也是较为简单的。
我们尝试将 向下取整 转换成 向上取值。
对于左边界
设
最后一步是因为
注意到
另外一个拓展就是我们尝试进行 高维数论分块。
当和式中出现若干个下取整,形如
时,我们只需要每次令
即可。最后要对
时间复杂度
因为我们把存在
9.3 例题
几道简单题。
CF1603C Extreme Extension
对于一个数列,定义一次操作为选择数列中任何一个元素
,将这个位置替换成两个正整数 满足 定义一个数列的极端值为将这个数列变成 不降 数列的最小操作次数
给定一个长度为
的数列 ( ),让你求它的所有非空子段的极端值之和,对 取模 一个测试点中有多组数据,保证所有数据
的规模总和不超过
P2260 [清华集训2012] 模积和
求
mod 19940417 的值
。
P3579 [POI2014] PAN-Solar Panels
对于
组询问,每组询问给定四个整数 ,从区间 和 中任意选取两个整数 和 ,求 的最大值是多少。
, , 。
10 欧拉函数
欧拉函数是一类非常重要的数论函数。
10.1 定义与性质
欧拉函数的定义为在
显然,当
若
为质数,则 。
在
是 积性函数。即若 。
设与
因为
因此,对于每一个
根据算术基本定理,设
被唯一分解为 ,则 。
根据前面两个性质
得证。
当然,我们还可以用 容斥原理 来证明。
考虑去掉所有被至少一个
其中
这样我们也可以反推到积性函数上面。设
于是有
若
,则 。
根据前面
于是我们也容易得到这样的结论:
若
欧拉反演:
考虑
我们把法里级数中非零的部分提出来,例如当
约分得到
我们可以根据它们的分母将这些分数分组,就可以得到
对此,我们如何就是呢?
发现
也就是
这就是 欧拉反演。
同时它可以被表示成狄利克雷卷积的形式,也就是
。
假设
容易发现这种情况显然是满足条件的,即
使得
且 的 的个数为 。
使得这个条件成立的条件可以被表示成
如果我们把条件转移到
10.2 线性筛欧拉函数
根据上面欧拉函数的计算式,我们可以得到它的线性筛法,若
于是代码就是简单的
phi[1]=1;
for(int i=2;i<N;i++){
if(!vis[i]) p[++cnt]=i,phi[i]=i-1;
for(int j=1;j<=cnt&&i*p[j]<N;j++){
vis[i*p[j]]=true;
if(i%p[j]==0){
phi[i*p[j]]=phi[i]*p[j];
break;
}
else phi[i*p[j]]=phi[i]*(p[j]-1);
}
}
10.3 欧拉定理
欧拉发现费马的定理可以用下面的方式推广到非素数的模:
这和费马小定理是相当类似的,于是我们可以尝试用类似于费马小定理的方法来证明它。
将它们相乘并用
而当
10.4 扩展欧拉定理
再把欧拉定理进行一个扩展,可以得到
感性理解:
当
令这个
根据欧拉定理,
比较抽象,非常感性。
严谨证明 摘自 OI-wiki:
-
的从 次, 次到 次幂模 构成的序列中,存在 ,使得前 个数(即从 到 )互不相同,从 开始,每 个数就循环依次。在 3.5 中,我们对这个遗留问题进行了讨论。
-
为素数的情况,该是成立。证明:
-
若模
不能被 整除,那么有 成立,根据欧拉定理容易证明。 -
若模
能被 整除,那么存在 和 使得 ,且 成立。所以根据欧拉定理有 。
又由于
,所以根据欧拉函数的计算式,容易得到 ,即我们有 。所以
,即 ,两边同时乘 ,得 ,因为 。所以
中素因子 的次数 满足: 。我们可以简单变换形式,得到 推论:又由于
,所以 。所以因为
,故有:所以
即
。 -
-
为素数幂的情况,该式成立。证明:
不妨令
,是否依然有 ?答案是肯定的,由于第一点可知存在
使得 ,所以 ,所以令 时,我们能有 。此时有关系
且 ,且 ,由 与 的关系,依然可以得到 。 -
为合数的情况,该是成立。证明:
只证
拆成两个素数的幂的情况,大于两个的用数学归纳法可证。设
,其中 ,而 的循环长度为 ;则
,由于 ,那么 ,所以 ;由
与 的关系,依然可以得到 。
证毕。乐。
于是你就可以通过 P5091 【模板】扩展欧拉定理。代码。
10.5 例题
煎蛋题。
P2158 [SDOI2008] 仪仗队
11 莫比乌斯函数
尝试按照《具体数学》的顺序引入莫比乌斯反演。
11.1 引入
莫比乌斯函数
来定义。
这个式子实际上是一个递归式,左边的
在 1857 年,Richard Dedekind 和 Joseph Liouville 注意到了如下重要的 “反演原理”。
这就是我们熟知的 莫比乌斯反演。
发现这个东西的证明是容易的:如果
反之,如果
于是我们在这里就证明了 莫比乌斯反演。
而莫比乌斯反演还可以被表示成
现在回到最开始的那个式子,发现其实就是
当
由此推出
于是,我们就推出了
如果我们把之前的欧拉反演视作
也就是
接下来我们尝试用容斥的角度来理解莫比乌斯函数。
设
就是对
11.2 线性筛莫比乌斯函数
根据定义是,线性筛莫比乌斯函数是非常容易的。
mu[1]=1;
for(int i=2;i<=maxn;i++){
if(!vis[i]) p[++cnt]=i,mu[i]=-1;
for(int j=1;j<=cnt&&i*p[j]<=maxn;j++){
vis[i*p[j]]=true;
if(i%p[j]==0) break;
mu[i*p[j]]=-mu[i];
}
}
11.3 常见技巧
对于 莫比乌斯反演,我们常常用式子
它把
于是我们就可以这样解决下面的式子:
相当于对最大公约数为
而这个式子是非常容易用整除分块完成。
11.4 例题
接下来,我们来将一些例题,由于比较多,所以单独提出来了。部分例题会用到上面的常见技巧。
前面的一些例题中有些写得比较简略,有些可能是👻乐团中的子问题,如果看不懂可以在那道题题解中找一下,或者直接问我。
P2522 [HAOI2011] Problem b
对于给出的
个询问,每次求有多少个数对 ,满足 , ,且 , 函数为 和 的最大公约数。
, , 。
首先可以对式子做一个二维差分,转换成求下界为
我们把
这就是一个经典的莫反式子,直接作就可以得到
于是就可以用整除分块在
子问题:P3455 [POI2007] ZAP-Queries。
P2257 YY的GCD
给定
,求 , 且 为质数的 有多少对。
, 。
问题就是求
我们用莫反相关知识对其进行推导可以得到
于是我们令
前面可以直接用整除分块完成,而后面的
但细想发现这是可以和莫比乌斯函数一起在线性筛中解决的,所以我们就可以直接预处理出其前缀和,从而做到时间复杂度
双倍经验:P2568 GCD。
SP5971 LCMSUM - LCM Sum
次询问,每次询问给出 ,求
。
首先我们考虑把
我们设
只要我们求出这个,就可以计算上面的式子了。
首先我们可以暴力用莫反计算,也就是
我们发现后面可以被拆成
还有另外一种方法就是对与
而显然
现在原式就变成了
而这里的
P4318 完全平方数
次询问,每次询问给出一个 ,查询第 个不是完全平方数的正整数倍数的数。
。
设
怎么求
发现可以先去掉
相当于做容斥,那么我们就有
直接计算,时间复杂度
P1829 [国家集训队] Crash的数字表格 / JZPTAB
给定
,求
。
后面的式子就和 SP5971 一样了。直接做就是线性的。代码。
P3704 [SDOI2017] 数字表格
Doris 刚刚学习了 fibonacci 数列。用
表示数列的第 项,那么 Doris 用老师的超级计算机生成了一个
的表格, 第
行第 列的格子中的数是 ,其中 表示 的最大公约数。 Doris 的表格中共有
个数,她想知道这些数的乘积是多少。 答案对
取模。
, 。
稍微推一下式子就是
预处理一下括号内的东西,整除分块之后可以做到单次
P5221 Product
给定
,求
,时限 。
有些人喜欢在写完👻之后来写这道题。
话不多说,直接暴力推一下式子就是了。
前面的直接阶乘处理一下,后面的我们枚举 gcd 就可以得到
于是两次整除分块就可以线性完成了,代码相当好些。代码。
P3911 最小公倍数之和
对于
,求
。
给定具体的数是不好处理的,所以我们考虑先令
那么
我们令
于是只要我们令
其中
这样就做完了,查询时可以直接线性完成。代码。(暴力枚举计算
双倍经验:AGC038C LCMs。
CF1285F Classical?
有
个整数 ,求
。
求完了和,我们现在来求
考虑加上一个
容易发现这是可以用单调栈直接维护的,那么我们如何计算出是否有与
通过上面有关
只要这个值不为
时间复杂度
而对于
我们可以枚举每一个
发现如果
P3327 [SDOI2015] 约数个数和
给定
,求
。
在常见技巧中,我们列出了整个式子的转化。
于是式子就是
于是直接对这个式子莫反就可以得到
令
于是式子就是
直接整除分块可以做到时间复杂度
P4619 [SDOI2018] 旧试题
给定
求
。
更为厉害的约数个数和。
多了一个变量,但是按照之前的转化方式是类似的
那么我们就可以对这个东西进行莫反了
我们可以令
然后你发现这还是
真的没有用吗?这个式子在哪些地方是有值的呢?
容易发现我们首先需要保证
那么也就是对于任意的两个点,它们可以成为一组就意味着
这个做法看起来根本不可行,并且直接这样建边是
但是细想一下,这个条件是相当严的,要求
每一次我们枚举了质数的子集,而
下面给出这一部分的代码实现。
for(int g=1;g<=mx;g++)
for(int i=1;i*g<=mx;i++)
if(mu[i*g])
for(int j=i+1;1ll*i*j*g<=mx;j++)
if(mu[j*g]&&gcd(i,j)==1){
int u=i*g,v=j*g,w=i*j*g;
/*...*/
}
进而,我们发现这样的点对并不多,实测发现只有 760741 个!赢麻了!
那么现在成了怎么找三元环,可以尝试阅读 这篇博客。
也就是我们将无向边定向,每一条边从度数大的点连向度数小的,然后这样去枚举的时间复杂度就是
于是我们就做完了!直接枚举就可以了。
for(int x=1;x<=mx;x++){
for(auto i:G[x]) vis[i.fi]=i.se;
for(auto i:G[x]){
int y=i.fi,wxy=i.se;
for(auto j:G[y]){
if(!vis[j.fi]) continue;
int z=j.fi,wyz=j.se,wxz=vis[z],val=mu[x]*mu[y]*mu[z];
/*...*/
}
}
for(auto i:G[x]) vis[i.fi]=0;
}
这样就可以做到
存在把上面式子化成
三倍经验:CF236B Easy Number Challenge 和 CF235E Number Challenge。
P5518 [MtOI2019] 👻乐团 / 莫比乌斯反演基础练习题
东风谷 早苗(Kochiya Sanae)非常喜欢幽灵乐团的演奏,她想对她们的演奏评分。
因为幽灵乐团有
个人,所以我们可以用 个正整数 来表示出乐团演奏的分数,她们的演奏分数可以表示为 因为音乐在不同的部分会有不同的听觉感受,所以
会在 中发生变化,其中: 因为乐团的歌实在太好听了,导致分数特别高,所以她们的分数要对给定的正整数
取模。 因为有很多歌曲要演奏,所以早苗给出了
组询问。
再来一个三重循环的题。
众所周知,cdqz 机房有且仅有三个👻,它们是:
-
:有名字 的👻,也是👻乐团的团长!
-
:没名字 的👻。
-
:直接消失的👻,这样的👻常常从你身边飘过!
这三只👻就是 cdqz 机房的真神!/bx
题目名字都告诉你了,还是赶紧来推式子吧。
首先对于 Type = 0/1 是简单的,直接顺着思路往下面推就可以了,也就是这个样子:
对于 Type = 0,我们首先把
前面两项是可以直接用阶乘完成的,也就是说我们只需要求
这个时候你可以把上面的指数部分单独提出来用一个函数完成,这样就可以做到线性,是可以过的,也就是后面代码1给出的(把三个 Type 分开了写的)。
但是你推到后面会发现这其实也是 Type = 2 的子问题,而这个式子是可以
而中间那一块是可以直接用
恭喜你,获得了
现在我们来研究 Type = 1,用同样的思路,写出来的式子也是简单的
把前面两项提出来就是形如
预处理一下
而对于后面两项,也就是求
这个东西推导方式就较为类似了,笔者并没有把它推导
于是可以把指数部分提出来,就变成了求
然后两次整除分块就可以完成。
很明显和上面 Type = 0 的优化类似,这时可以做到
预处理一下
最后来搞 Type = 2 的情况,按照相同的思路先展开一下
对于前两项,其实就是(后面的
像往常一样,我们直接用
把
容易发现这两项都可以用整除分块做出来的,稍微预处理一下
而对于后面两项,也就是求
按照之前的经验,我们尝试去枚举
然后就有
你可以自己尝试推一下,就会发现只能用
那怎么办呢?
式子中还存在一个
上面的内容和之前推第一部分类似,而现在我们发现其实中间的
其实就是 Type = 0 中我们推的
惊奇地发现这两部分前面的那一块其实是一样的,于是它们就抵消啦!
所以我们需要求的就是
和
于是你就可以用两次整除分块线性求出来啦!
这样你就推出了本题的全部式子,写的时候注意细节,在指数上面的东西注意模的是
本题存在对于 Type = 2 的情况求
P3700 [CQOI2017] 小 Q 的表格
为了完成任务,小 Q 需要列一个表格,表格有无穷多行,无穷多列,行和列都从
开始标号。为了完成任务,表格里面每个格子都填了一个整数,为了方便描述,小 Q 把第 行第 列的整数记为 。为了完成任务,这个表格要满足一些条件:
- 对任意的正整数
,都要满足 ; - 对任意的正整数
,都要满足 。 为了完成任务,一开始表格里面的数很有规律,第
行第 列的数恰好等于 ,显然一开始是满足上述两个条件的。为了完成任务,小 Q 需要不断的修改表格里面的数,每当修改了一个格子的数之后,为了让表格继续满足上述两个条件,小 Q 还需要把这次修改能够波及到的全部格子里都改为恰当的数。由于某种神奇的力量驱使,已经确保了每一轮修改之后所有格子里的数仍然都是整数。为了完成任务,小 Q 还需要随时获取前 行前 列这个有限区域内所有数的和是多少,答案可能比较大,只需要算出答案 之后的结果。
, , 。
建议盯着式子发呆。
首先,由于
而你把第二个式子换一种写法就是
多写几个就变成了
那么换一种写法就是
我们发现这个东西和辗转相除求 gcd 是一样的,那么我们就可以得到这个很好的性质
那么
于是我们要求的答案其实就是
那么现在就是希望求
我们需要预处理出所有的
但是按照上一道题我们推式子的方法,发现单次是
这时我们考虑在 SP5971 中求的东西
带进去就是
于是这个东西就可以线性递推出来了,只要我们能够动态维护
但是我们发现上面所说的
具体每个块每个元素维护块内的前缀和,在每一个块末尾维护块间的前缀和,修改时暴力修改即可。
总时间复杂度
12 奇妙筛法
已开工。
12.1 杜教筛
12.2 Min-25 筛
Min_25 筛可以在
说形象点就是
要求:
是关于 的低阶多项式,或者可以快速求出。 可以快速求出。
以下
同时我们钦定
首先,我们尝试对其进行推导,分成质数和合数两个部分
前面的部分是枚举素数,而后面的部分就是对于每一个合数枚举了其最小的质因子以及其次数。
考虑第一块的质数部分如何求值?
由于我们在前面要求了
那么
现在问题就变成了求
我们只对一个
我们尝试用一个神秘的 dp,设
(至于它有什么用,你先别急)
我们考虑
不难想到用
由于
而仔细一想,发现我们把质数多算了一次,它一共出现了
其实后面的
这样一来,我们就有了递推式
发现这是可以依次 dp 下去的,只要我们知道
就可以用质数处理出来。(但是时间复杂度根本过不去啊?你先别急)
回到引入
我们发现如果
所以如果通过 dp 推出
假设我们现在已经可以快速算出
发现后面的那一块也是不好处理的,所以我们采用类似的 dp 方式,设
其中
而后面的
这样的递归是可以直接求的,时间复杂度为
现在我们回到对
发现直接处理出
这不就是类似于整除分块的东西么?
我们在第一章证明过
所以在整个过程中只会用到
在第九章我们知道个数是
再来讲一讲具体的实现细节,我们用
同时,由于
于是就可以通过 P5325 【模板】Min_25 筛 了。代码。
const ll H=1e9+7;
int cnt=0,tot=0;
ll p[N],g1[N],g2[N],s1[N],s2[N],w[N],id1[N],id2[N],B,n,i6,i2;
bool vis[N];
ll qpow(ll a,ll b){
ll res=1;
while(b){if(b&1) res=res*a%H;a=a*a%H,b>>=1;}
return res;
}
ll S1(ll x){x%=H;return x*(x+1)%H*i2%H;}
ll S2(ll x){x%=H;return x*(x+1)%H*(2*x+1)%H*i6%H;}
ll pos(ll x){return x<=B?id1[x]:id2[n/x];}
ll S(ll x,int y){
if(x<=p[y]) return 0;
int k=pos(x),res=((g2[k]-g1[k]-s2[y]+s1[y])%H+H)%H;
for(int i=y+1;i<=cnt&&p[i]*p[i]<=x;i++)
for(ll j=1,nw=p[i];nw<=x;++j,nw*=p[i]){
ll z=nw%H;
res=(res+(z-1)*z%H*(S(x/nw,i)+(j>1))%H)%H;
}
return res;
}
int main(){
ios::sync_with_stdio(false);cin.tie(0);cout.tie(0);
cin>>n,B=sqrt(n);
i2=qpow(2,H-2),i6=qpow(6,H-2);
for(int i=2;i<=B;i++){
if(!vis[i]){
p[++cnt]=i;
s1[cnt]=(s1[cnt-1]+i)%H;
s2[cnt]=(s2[cnt-1]+1ll*i*i%H)%H;
}
for(int j=1;j<=cnt&&i*p[j]<=B;j++){
vis[i*p[j]]=true;
if(i%p[j]==0) break;
}
}
for(ll l=1,r;l<=n;l=r+1){
r=min(n/(n/l),n);
w[++tot]=n/l;
g1[tot]=S1(n/l)-1,g2[tot]=S2(n/l)-1;
if(n/l<=B) id1[n/l]=tot;
else id2[n/(n/l)]=tot;
}
for(int i=1;i<=cnt;i++){//g 数组的 dp,用了滚动数组
for(int j=1;p[i]*p[i]<=w[j]&&j<=tot;j++){
int cur=pos(w[j]/p[i]);
g1[j]=(g1[j]-p[i]*(g1[cur]-s1[i-1])%H+H)%H;
g2[j]=(g2[j]-p[i]*p[i]%H*(g2[cur]-s2[i-1])%H+H)%H;
}
}
cout<<(S(n,0)+1)%H;
return 0;
}
现在我们换几个函数来筛,也就是我们最熟悉的
它们在质数处的取值也是相当容易的,
直接把两个分开写出来就是这样的
const int N=1e6+5;
int cnt=0,tot=0;
ll p[N],w[N],id1[N],id2[N],n,B;
bool vis[N>>2];
int pos(ll x){return x<=B?id1[x]:id2[n/x];}
void Pr(int n){
memset(vis,0,sizeof(vis));
for(int i=2;i<=n;i++){
if(!vis[i]) p[++cnt]=i;
for(int j=1;j<=cnt&&i*p[j]<=n;j++){
vis[i*p[j]]=true;
if(i%p[j]==0) break;
}
}
}
namespace Mu{
ll g[N];
ll S(ll x,int y){
if(x<=p[y]) return 0;
ll res=-g[pos(x)]+y;
for(int i=y+1;i<=cnt&&p[i]*p[i]<=x;i++)
res-=S(x/p[i],i);
return res;
}
void SOL(){
for(int i=1;i<=tot;i++) g[i]=w[i]-1;
for(int i=1;i<=cnt;i++){
for(int j=1;p[i]*p[i]<=w[j]&&j<=tot;j++){
int cur=pos(w[j]/p[i]);
g[j]+=i-g[cur]-1;
}
}
cout<<(S(n,0)+1)<<'\n';
}
}
namespace Phi{
ll g1[N],g2[N],s[N];
ll S(ll x,int y){
if(x<=p[y]) return 0;
int k=pos(x);
ll res=g1[k]-g2[k]-(s[y]-y);
for(int i=y+1;i<=cnt&&p[i]*p[i]<=x;i++){
for(ll j=1,nw=1;nw*p[i]<=x;++j,nw*=p[i])
res+=nw*(p[i]-1)*(S(x/(nw*p[i]),i)+(j>1));
}
return res;
}
void SOL(){
for(int i=1;i<=cnt;i++) s[i]=s[i-1]+p[i];
for(int i=1;i<=tot;i++) g1[i]=1ll*w[i]*(w[i]+1)/2-1,g2[i]=w[i]-1;
for(int i=1;i<=cnt;i++)
for(int j=1;p[i]*p[i]<=w[j]&&j<=tot;j++){
int cur=pos(w[j]/p[i]);
g1[j]-=p[i]*(g1[cur]-s[i-1]);
g2[j]+=i-1-g2[cur];
}
cout<<(S(n,0)+1)<<' ';
}
}
void SOLVE(){
cin>>n,B=sqrt(n),cnt=tot=0;
memset(vis,0,sizeof(vis));
Pr(B);
for(ll l=1,r;l<=n;l=r+1){
r=min(n,n/(n/l));
w[++tot]=n/l;
if(n/l<=B) id1[n/l]=tot;
else id2[n/(n/l)]=tot;
}
Phi::SOL(),Mu::SOL();
}
int main(){
/*2024.4.26 H_W_Y P4213 【模板】杜教筛 Min_25 筛*/
ios::sync_with_stdio(false);cin.tie(0);cout.tie(0);
int _;cin>>_;
while(_--) SOLVE();
return 0;
}
但结果并不理想:(前两个是 Min_25,最下面的是杜教筛)
于是我们开始了常数优化。
优化也是简单的,你把预处理放在一起,就可以跑到 2.95s:
进而,把两个函数一起做(再把不需要 long long 的地方换成 int),于是就可以做到 1.44s!!!
代码。
const int N=1e6+5;
int cnt=0,tot=0,p[N],n,B,w[N],id1[N],id2[N],s[N];
ll g1[N],g2[N];
bool vis[N>>2];
inline int pos(int x){return x<=B?id1[x]:id2[n/x];}
inline void Pr(int n){
for(int i=2;i<=n;i++){
if(!vis[i]) p[++cnt]=i,s[cnt]=s[cnt-1]+i;
for(int j=1;j<=cnt&&i*p[j]<=n;j++){
vis[i*p[j]]=true;
if(i%p[j]==0) break;
}
}
}
inline pll S(int x,int y){
if(x<=p[y]) return {0,0};
int k=pos(x);
pll res={g1[k]-g2[k]-(s[y]-y),-g2[k]+y};
for(int i=y+1;i<=cnt&&p[i]*p[i]<=x;i++)
for(int j=1,nw=1;1ll*nw*p[i]<=x;++j,nw*=p[i]){
pll cur=S(x/(nw*p[i]),i);
if(nw==1) res.se-=cur.se;
res.fi+=1ll*nw*(p[i]-1)*(cur.fi+(j>1));
}
return res;
}
inline void SOLVE(){
cin>>n,B=sqrt(n),tot=0;
for(ll l=1,r;l<=n;l=r+1){
r=n/(n/l);
w[++tot]=n/l;
g1[tot]=1ll*(n/l)*(n/l+1)/2-1,g2[tot]=n/l-1;
if(n/l<=B) id1[n/l]=tot;
else id2[n/(n/l)]=tot;
}
for(int i=1;i<=cnt&&p[i]<=B;i++)
for(int j=1;p[i]*p[i]<=w[j]&&j<=tot;j++){
int cur=pos(w[j]/p[i]);
g1[j]-=1ll*p[i]*(g1[cur]-s[i-1]);
g2[j]+=i-1-g2[cur];
}
pll res=S(n,0);
cout<<(res.fi+1)<<' '<<(res.se+1)<<'\n';
}
int main(){
/*2024.4.26 H_W_Y P4213 【模板】杜教筛 Min_25 筛*/
ios::sync_with_stdio(false);cin.tie(0);cout.tie(0);
int _;cin>>_;Pr(46340);
while(_--) SOLVE();
return 0;
}
12.3 Powerful Number 筛
12.4 扩展埃氏筛
13 阶与原根
模
的简化剩余系在模 乘法意义下封闭且构成群。容易证明其满足封闭性和结合律,存在逆元和单位元。将群相关定义应用在其上,得到阶和原根(循环群和生成元)的概念。
——Alex_Wei
13.1 阶
由欧拉定理可知,对于
因此满足同余式
阶有一些性质。
模 两两不同余。
这和我们在 3.5 里面探讨的遗留问题较为相似。
考虑反证,假设存在
但是显然有
若
,则 。
设
若
与阶的定义矛盾(最小性)。得证。
设
,则 的充分必要条件是
-
必要性:因为
和 可知由于上面的性质,我们知道
所以我们有
即
。 -
充分性:因为
那么
所以
又由于
,同理可得
所以
而我们有
故
所以
得证。
设
,则
很重要的一个性质。
注意到
而由于
所以
故
得证。
容易发现其实上面的证明方法都挺类似的。/kk
那如何求
-
直接使用 BSGS 求
的最小正整数解。时间复杂度 。 -
根据阶的性质,对于
,必然有 。因此,我们考虑先让
,然后对于 的每个质因子 ,用 不断试除 知道 无法整除 或者 。时间复杂度
,因为需要分解质因数和求解欧拉函数,用 Pollard_rho 可以优化为 。若
不大,则可以 线性处理每个数最小质因子,单次查询即可 。
13.2 原根
原根——数论的神!
设
我们需要原根的原因是它可以 生成 模
原根判定定理:设
这和我们上面讲的求阶的方法非常像。
-
必要性:显然。
-
充分性:因为
的所有因子取遍了 除了其本身的所有因子,所以若存在 ,使得 ,则必然存在 的质因子 使得 ,这就说明存在 ,与假设矛盾。得证
原根个数:若一个数
若
所以如果
而满足
原根存在定理:一个数存在原根当且仅当
可以当一个结论背住,具体证明比较复杂,可见 OI-wiki。
最小原根的范围估计:
在 1959 年,中国数学家王元证明了最小原根的级别为
关于原根,还有一个有趣且重要的性质:
若
这是因为对于任何一个
同时,
这也是为什么对于特定模数,我们可以使用原根代替
于是就得到了 NTT。
13.3 例题
放两道简单题。
P5605 小 A 与两位神仙
这个游戏的规则是这样的:两个神仙先选定一个正整数
,保证 是一个奇质数 的正整数次幂。然后进行 轮游戏,每轮中造梦者选定一个正整数 ,杰瑞米选定一个正整数 ,保证 ,即 与 互质, 与 互质,接下来询问小 A 是否存在非负整数 使得 。 神仙们说小 A 只有在每一轮游戏中都回答正确才能回到正常的生活中,不得已之下他只好求助于聪明的你。
, , 。
注意到模数
我们尝试将
14 高次剩余问题
15 卢卡斯定理
后记
参考文献:
Pollard_Rho 算法 - cjTQX - 博客园 (cnblogs.com)
《具体数学》第四章
本文作者:H_W_Y
本文链接:https://www.cnblogs.com/H-W-Y/p/18140410/NumberTheory
版权声明:本作品采用知识共享署名-非商业性使用-禁止演绎 2.5 中国大陆许可协议进行许可。
【推荐】编程新体验,更懂你的AI,立即体验豆包MarsCode编程助手
【推荐】凌霞软件回馈社区,博客园 & 1Panel & Halo 联合会员上线
【推荐】抖音旗下AI助手豆包,你的智能百科全书,全免费不限次数
【推荐】博客园社区专享云产品让利特惠,阿里云新客6.5折上折
【推荐】轻量又高性能的 SSH 工具 IShell:AI 加持,快人一步