算法竞赛如何训练数论这一块？
算法竞赛中的数论训练，就像在迷宫里寻找宝藏，每一步都需要扎实的理论基础和巧妙的计算技巧。要想在这片领域游刃有余，就得循序渐进，把基础打牢，然后逐步攻克难题。
第一步：夯实基础，筑牢地基
核心概念的理解：
整除与模运算： 这是数论的基石，要深刻理解 $a pmod{m}$ 的含义，以及整除的性质。比如，如果 $a|b$ 且 $b|c$，那么 $a|c$。这些看似简单的性质，在很多复杂的证明和算法设计中起着至关重要的作用。
同余方程与模运算性质： 掌握同余的传递性、加法、减法、乘法性质。比如，若 $a equiv b pmod{m}$ 且 $c equiv d pmod{m}$，则 $a+c equiv b+d pmod{m}$ 和 $ac equiv bd pmod{m}$。这是后面学习线性同余方程等内容的基础。
素数与合数： 什么是素数？什么是合数？了解素数的定义以及素数定理（虽然不需要严格证明，但知道素数分布的大致规律很有帮助）。
最大公约数 (GCD) 与最小公倍数 (LCM)： 理解 GCD 和 LCM 的定义，以及它们之间 $GCD(a,b) imes LCM(a,b) = |a imes b|$ 的关系。
欧几里得算法 (Euclidean Algorithm)： 这是计算 GCD 的经典算法，熟练掌握其原理和实现。很多数论问题都离不开高效的 GCD 计算。
经典定理的学习：
费马小定理 (Fermat's Little Theorem)： 如果 $p$ 是一个素数，那么对于任意整数 $a$，$a^p equiv a pmod{p}$。特别地，如果 $a$ 不是 $p$ 的倍数，则 $a^{p1} equiv 1 pmod{p}$。这是后面学习欧拉定理和模幂运算的关键。
欧拉定理 (Euler's Totient Theorem)： 对于任意互质整数 $a$ 和 $n$，有 $a^{phi(n)} equiv 1 pmod{n}$，其中 $phi(n)$ 是欧拉函数，表示小于等于 $n$ 的正整数中与 $n$ 互质的数的个数。这是费马小定理的推广。
中国剩余定理 (Chinese Remainder Theorem, CRT)： 当我们有若干个线性同余方程组，且模数两两互质时，CRT 可以帮我们找到一个解。这是解决一些涉及多个模数的复杂问题的利器。
练习与实现：
GCD 的各种应用： 比如约分分数，判断是否互质等。
模幂运算 (Modular Exponentiation)： 高效计算 $a^b pmod{m}$，这在密码学和组合数学中非常常见。可以使用快速幂算法（Binary Exponentiation）。
欧拉函数的计算： 了解欧拉函数的性质，比如 $phi(p^k) = p^k p^{k1}$ 和 $phi(mn) = phi(m)phi(n)$（当 $m, n$ 互质时）。学习如何通过素因数分解计算任意数的欧拉函数值。
简单同余方程的求解： 比如 $ax equiv b pmod{m}$。
第二步：进阶技巧，深化理解
线性同余方程组与中国剩余定理的深入应用：
掌握如何求解形如 $a_1 x equiv b_1 pmod{m_1}$，$a_2 x equiv b_2 pmod{m_2}$，... 的线性同余方程组。
了解当模数不互质时，如何处理 CRT。这通常需要先将方程化简，或者通过其他技巧合并。
通过大量练习，体会 CRT 在解决“还原问题”和“组合问题”时的强大威力。
素性测试 (Primality Testing)：
了解 MillerRabin 素性测试，这是目前算法竞赛中常用的高效概率性素性测试算法。理解其背后的数学原理（费马小定理和二次探测定理）。
虽然不需要自己从头实现 MillerRabin，但要理解其思想，能在需要时调用相关库函数或理解现有代码。
模逆元 (Modular Multiplicative Inverse)：
理解模逆元存在的条件（即 $a$ 与 $m$ 互质）。
学习使用扩展欧几里得算法 (Extended Euclidean Algorithm) 计算模逆元。
了解费马小定理和欧拉定理也可以用来计算模逆元（当模数是素数时，费马小定理更方便）。
数论函数 (Arithmetic Functions)：
了解一些常见的数论函数，如欧拉函数 $phi(n)$、莫比乌斯函数 $mu(n)$、约数和函数 $sigma(n)$、单位函数 $epsilon(n)$ (或称狄利克雷卷积的单位元)。
理解狄利克雷卷积 (Dirichlet Convolution) 的概念：$(f g)(n) = sum_{d|n} f(d)g(n/d)$。这是理解和推导许多数论函数性质的关键。
学习莫比乌斯反演 (Möbius Inversion Formula)：如果 $g(n) = sum_{d|n} f(d)$，那么 $f(n) = sum_{d|n} mu(n/d) g(d)$。这是解决许多数论问题，特别是涉及 GCD 求和的问题的利器。
生成函数 (Generating Functions)：
虽然生成函数更偏向组合数学，但在数论中也有其身影，尤其是在处理一些数列的性质时。了解其基本思想和应用，例如利用生成函数推导数论函数之间的关系。
第三步：实战演练，攻克难关
从简单题入手： 在牛客网、洛谷、Codeforces 等平台上，搜索“数论”标签的题目，从普及组或提高组的简单题开始练起。
类型题专项训练：
GCD/LCM 相关： 包含约数个数、约数和、GCD 求和等题目。
模运算/同余方程： 涉及模幂、同余方程求解、模逆元、线性同余方程组等。
中国剩余定理： 各种需要 CRT 解题的题目。
数论函数/莫比乌斯反演： 很多 ACM/ICPC 风格的题目会用到这些。
素数相关： 质因数分解、素数筛法（埃氏筛法、线性筛法）、素性测试等。
代码实现能力：
熟练掌握各种数论算法的 C++ 实现，包括快速幂、扩展欧几里得、中国剩余定理模板、埃氏筛、线性筛等。
注意代码的常数优化和边界处理，特别是在处理大数和模运算时。
理解题意与模型转化：
很多数论题目的形式可能比较抽象，需要将其转化为熟悉的数论模型。例如，很多组合计数问题可能背后隐藏着数论的规律。
学会从题目中提取关键信息，比如“求和”、“最大值”、“个数”、“存在性”等，并思考哪种数论工具适用于解决这些问题。
阅读他人题解与学习优秀代码：
遇到难题时，不要卡太久。参考优秀的题解和代码，理解别人是怎么思考和解决问题的。但要避免盲目复制代码，而是要消化吸收其思想。
分析代码实现细节，学习别人的代码风格和技巧。
参与模拟赛：
模拟赛是检验学习成果的最佳方式。在模拟赛中，强迫自己按照考试的节奏去思考和解决问题，包括时间分配、题目选择、题目的难度评估等。
赛后认真复盘，找出在数论题目上的薄弱环节，并针对性地进行补充训练。
进阶提升
数论在其他领域的应用： 了解数论在密码学（RSA、ElGamal 等）、编码理论、图论等领域的应用，这有助于开阔思路。
更深入的数论概念： 如果有精力，可以了解一些更深入的数论概念，例如二次剩余、平方数、丢番图方程等，这些在一些高级竞赛中可能会出现。
训练心态
数论的学习曲线可能比较陡峭，需要耐心和坚持。不要因为一开始遇到困难就气馁。把每一次的练习都看作是一次学习和进步的机会，不断积累经验，你会发现自己在这片奇妙的数学世界里越走越远。
祝你在算法竞赛的数论之路上取得成功！
第一步：夯实基础，筑牢地基
核心概念的理解：
整除与模运算： 这是数论的基石，要深刻理解 $a pmod{m}$ 的含义，以及整除的性质。比如，如果 $a|b$ 且 $b|c$，那么 $a|c$。这些看似简单的性质，在很多复杂的证明和算法设计中起着至关重要的作用。
同余方程与模运算性质： 掌握同余的传递性、加法、减法、乘法性质。比如，若 $a equiv b pmod{m}$ 且 $c equiv d pmod{m}$，则 $a+c equiv b+d pmod{m}$ 和 $ac equiv bd pmod{m}$。这是后面学习线性同余方程等内容的基础。
素数与合数： 什么是素数？什么是合数？了解素数的定义以及素数定理（虽然不需要严格证明，但知道素数分布的大致规律很有帮助）。
最大公约数 (GCD) 与最小公倍数 (LCM)： 理解 GCD 和 LCM 的定义，以及它们之间 $GCD(a,b) imes LCM(a,b) = |a imes b|$ 的关系。
欧几里得算法 (Euclidean Algorithm)： 这是计算 GCD 的经典算法，熟练掌握其原理和实现。很多数论问题都离不开高效的 GCD 计算。
经典定理的学习：
费马小定理 (Fermat's Little Theorem)： 如果 $p$ 是一个素数，那么对于任意整数 $a$，$a^p equiv a pmod{p}$。特别地，如果 $a$ 不是 $p$ 的倍数，则 $a^{p1} equiv 1 pmod{p}$。这是后面学习欧拉定理和模幂运算的关键。
欧拉定理 (Euler's Totient Theorem)： 对于任意互质整数 $a$ 和 $n$，有 $a^{phi(n)} equiv 1 pmod{n}$，其中 $phi(n)$ 是欧拉函数，表示小于等于 $n$ 的正整数中与 $n$ 互质的数的个数。这是费马小定理的推广。
中国剩余定理 (Chinese Remainder Theorem, CRT)： 当我们有若干个线性同余方程组，且模数两两互质时，CRT 可以帮我们找到一个解。这是解决一些涉及多个模数的复杂问题的利器。
练习与实现：
GCD 的各种应用： 比如约分分数，判断是否互质等。
模幂运算 (Modular Exponentiation)： 高效计算 $a^b pmod{m}$，这在密码学和组合数学中非常常见。可以使用快速幂算法（Binary Exponentiation）。
欧拉函数的计算： 了解欧拉函数的性质，比如 $phi(p^k) = p^k p^{k1}$ 和 $phi(mn) = phi(m)phi(n)$（当 $m, n$ 互质时）。学习如何通过素因数分解计算任意数的欧拉函数值。
简单同余方程的求解： 比如 $ax equiv b pmod{m}$。
第二步：进阶技巧，深化理解
线性同余方程组与中国剩余定理的深入应用：
掌握如何求解形如 $a_1 x equiv b_1 pmod{m_1}$，$a_2 x equiv b_2 pmod{m_2}$，... 的线性同余方程组。
了解当模数不互质时，如何处理 CRT。这通常需要先将方程化简，或者通过其他技巧合并。
通过大量练习，体会 CRT 在解决“还原问题”和“组合问题”时的强大威力。
素性测试 (Primality Testing)：
了解 MillerRabin 素性测试，这是目前算法竞赛中常用的高效概率性素性测试算法。理解其背后的数学原理（费马小定理和二次探测定理）。
虽然不需要自己从头实现 MillerRabin，但要理解其思想，能在需要时调用相关库函数或理解现有代码。
模逆元 (Modular Multiplicative Inverse)：
理解模逆元存在的条件（即 $a$ 与 $m$ 互质）。
学习使用扩展欧几里得算法 (Extended Euclidean Algorithm) 计算模逆元。
了解费马小定理和欧拉定理也可以用来计算模逆元（当模数是素数时，费马小定理更方便）。
数论函数 (Arithmetic Functions)：
了解一些常见的数论函数，如欧拉函数 $phi(n)$、莫比乌斯函数 $mu(n)$、约数和函数 $sigma(n)$、单位函数 $epsilon(n)$ (或称狄利克雷卷积的单位元)。
理解狄利克雷卷积 (Dirichlet Convolution) 的概念：$(f g)(n) = sum_{d|n} f(d)g(n/d)$。这是理解和推导许多数论函数性质的关键。
学习莫比乌斯反演 (Möbius Inversion Formula)：如果 $g(n) = sum_{d|n} f(d)$，那么 $f(n) = sum_{d|n} mu(n/d) g(d)$。这是解决许多数论问题，特别是涉及 GCD 求和的问题的利器。
生成函数 (Generating Functions)：
虽然生成函数更偏向组合数学，但在数论中也有其身影，尤其是在处理一些数列的性质时。了解其基本思想和应用，例如利用生成函数推导数论函数之间的关系。
第三步：实战演练，攻克难关
从简单题入手： 在牛客网、洛谷、Codeforces 等平台上，搜索“数论”标签的题目，从普及组或提高组的简单题开始练起。
类型题专项训练：
GCD/LCM 相关： 包含约数个数、约数和、GCD 求和等题目。
模运算/同余方程： 涉及模幂、同余方程求解、模逆元、线性同余方程组等。
中国剩余定理： 各种需要 CRT 解题的题目。
数论函数/莫比乌斯反演： 很多 ACM/ICPC 风格的题目会用到这些。
素数相关： 质因数分解、素数筛法（埃氏筛法、线性筛法）、素性测试等。
代码实现能力：
熟练掌握各种数论算法的 C++ 实现，包括快速幂、扩展欧几里得、中国剩余定理模板、埃氏筛、线性筛等。
注意代码的常数优化和边界处理，特别是在处理大数和模运算时。
理解题意与模型转化：
很多数论题目的形式可能比较抽象，需要将其转化为熟悉的数论模型。例如，很多组合计数问题可能背后隐藏着数论的规律。
学会从题目中提取关键信息，比如“求和”、“最大值”、“个数”、“存在性”等，并思考哪种数论工具适用于解决这些问题。
阅读他人题解与学习优秀代码：
遇到难题时，不要卡太久。参考优秀的题解和代码，理解别人是怎么思考和解决问题的。但要避免盲目复制代码，而是要消化吸收其思想。
分析代码实现细节，学习别人的代码风格和技巧。
参与模拟赛：
模拟赛是检验学习成果的最佳方式。在模拟赛中，强迫自己按照考试的节奏去思考和解决问题，包括时间分配、题目选择、题目的难度评估等。
赛后认真复盘，找出在数论题目上的薄弱环节，并针对性地进行补充训练。
进阶提升
数论在其他领域的应用： 了解数论在密码学（RSA、ElGamal 等）、编码理论、图论等领域的应用，这有助于开阔思路。
更深入的数论概念： 如果有精力，可以了解一些更深入的数论概念，例如二次剩余、平方数、丢番图方程等，这些在一些高级竞赛中可能会出现。
训练心态
数论的学习曲线可能比较陡峭，需要耐心和坚持。不要因为一开始遇到困难就气馁。把每一次的练习都看作是一次学习和进步的机会，不断积累经验，你会发现自己在这片奇妙的数学世界里越走越远。
祝你在算法竞赛的数论之路上取得成功！
网友意见
题主你好，我是浙江师范大学2016级数学专业的ACM集训队队员，曾在计蒜客实习，对ACM与OI竞赛都有了解，希望我的回答对你有所帮助。
如果你是大学生，可以考虑了解一下抽象代数，可以从一个奇妙的角度来理解数论。
下面是我整理的数论算法，全文非常长，在 markdown 格式下约 6w 字符，建议放到收藏夹吃灰。
文中某些代码仅为了帮助读者理解算法而写，建议理解后自己写模板。
文中不包含积性函数前缀和筛，例如min25筛，因为笔者认为这样的知识点太偏太难，考起来只是大佬们的自嗨。（划掉）（更新）听评论区说想要这些杂技内容，那我就把模板贴这了，没有原理解析，也不会有的。不过多项式求逆是没有的，因为我没遇到过，模板也就没写过。
另外，厚颜无耻地推广自己的bilibili账号，里面有我校自用的ACM课程，我要退役了，你们加油！
数论算法详解
1. 整数的取余运算
定义 带余除法 设 是整数，且 ，使得
且 ，称 为商， 为余数
显然带余除法中的商和余数都是唯一的，在下文中将商记为 ，将余数记为 ，“/”与“%”的运算优先级与乘除法相同，当然，在C语言中二者分别对应 a/b 与 a%b
1.1 整数的加减乘在 中的运算
1.1.1 定理 整数的加、减、乘对于取余的右分配律
注意，在计算减法时，通常需要加 ，防止变成负数
在计算乘法时，如果 较大（但不超过64位整数范围），可以使用快速乘法进行计算，原理与快速幂运算类似，复杂度为
ll fastMul(ll a,ll b,ll p){ a%=p; ll ans=0; while(b>0){ if(b&1)ans=(ans+a)%p; b>>=1; a=(a+a)%p; } return ans; } 如果需要更快的快速乘法，可以使用 long double 数据类型进行计算，复杂度为
ll modMul(ll a,ll b,ll p){ if (p<=1000000000) return a*b%p; else if (p<=1000000000000ll) return (((a*(b>>20)%p)<<20)+(a*(b&((1<<20)-1))))%p; else { ll d=(ll)floor(a*(long double)b/p+0.5); ll ret=(a*b-d*p)%p; if (ret<0) ret+=p; return ret; } } 1.2 整数的除法在 中的运算( 为素数)
虽然取余运算对于“+”、“-”、“×”不难，但通常情况下
如何计算 ？我们有时可以找一个乘法逆元 ，使得 ，那么就有
如果 是素数，有下面的定理
1.2.1 定理 费马小定理 设 是一个整数， 是一个素数，二者互质，那么
将上式改写一下，得到
因此取 即可，一般需要用快速幂计算，但要注意，与除数不能为0类似，要保证
ll inv(ll a,ll p){ return fpow(a,p-2); } 1.3 整数的除法在 中的运算
上面的方法给出了模数为素数的解决方案，如果模数不是素数，可以用下面的方法
1.3.1 定理 若 ，则
这样可以不使用逆元来求 中的除法
其他逆元相关内容将在下文中介绍
（如果你像我一样无聊，可以把整环 中的数写成这样一个类）
class mint{ static const ll mo=1e9+7; inline ll fpow(ll a,ll b){ll c=1;while(b){if(b&1)c=c*a%mo;b>>=1;a=a*a%mo;}return c;} inline ll inv(ll a){return fpow(a,mo-2);} inline ll norm(ll a){return a<0?(a%mo+mo):a%mo;} public: ll v; mint(){} mint(ll x):v(x){} mint &operator =(const ll b){v=norm(b);return *this;} mint &operator +=(const mint &b){v+=b.v;if(v>=mo)v-=mo;return *this;} mint &operator -=(const mint &b){v-=b.v;if(v<0)v+=mo;return *this;} mint &operator *=(const mint &b){v=v*b.v%mo;return *this;} mint &operator /=(const mint &b){v=v*inv(b.v)%mo;return *this;} mint operator + (const mint &b)const{mint a{*this};return a+=b;} mint operator - (const mint &b)const{mint a{*this};return a-=b;} mint operator * (const mint &b)const{mint a{*this};return a*=b;} mint operator / (const mint &b)const{mint a{*this};return a/=b;} mint operator ^ (const mint &b){return mint(fpow(v,b.v));} mint &operator++(){v++;if(v==mo)v-=mo;return *this;} mint &operator--(){v--;if(v==-1)v+=mo;return *this;} mint operator -()const{return mint(v?(mo-v):0);} bool operator ==(const mint &b)const{return v==b.v;} bool operator <(const mint &b)const{return v<b.v;} bool operator >(const mint &b)const{return v<b.v;} bool operator !=(const mint &b)const{return v!=b.v;} explicit operator ll()const{return v;} friend ostream &operator << (ostream &o,const mint & b){o<<b.v;return o;} friend istream &operator >> (istream &in,mint &a){in>>a.v;a.v%=mo;return in;} }; 2. 最大公因数与最小公倍数
2.1 最大公因数
顾名思义，最大公因数就是公因数中最大的那个，我们记 的最大公因数为 ，有如下性质
2.1.1 性质
2.1.2 性质
2.1.3 性质
2.1.4 性质
2.1.5 性质
2.1.1 性质是显然的，2.1.2 性质是辗转相减法的原理，2.1.3 性质可以视为2.1.2 性质的“一步到位”版本，2.1.4 性质指出多个数的最大公因数可以递推地进行求解，2.1.5 性质说明 对乘法有分配律
2.2 辗转相除法
根据2.1.3 性质，得到辗转相除法的参考代码模板
typedef long long ll; ll gcd(ll a, ll b){ return b?gcd(b,a%b):a; } 注意当 时，返回值为 gcd(b,a%b) 而不是 gcd(a%b,b) ，否则会不断递归导致栈溢出
2.3 最小公倍数
最小公倍数就是公倍数中最小的那个，我们记 的最小公倍数为 ，有如下性质
2.3.1 性质
下面是最小公倍数的参考代码模板
ll lcm(ll a, ll b){ return a/gcd(a,b)*b; } 注意是先除后乘，避免在中间过程中数据超出64位整数的范围
2.4 扩展欧几里得算法
考虑二元一次不定方程
其中 是已知的正整数，如何求出方程的解呢？
2.4.1 定理 上述方程有解的充要条件是 （ 是 的倍数）
可以理解为， 是 可以表示出的最小正整数
2.4.2 定理 方程 的所有解为
其中 是一组特解
2.4.3 定理 方程 的所有解为
其中 是方程 的一组特解
下面是参考代码模板：
ll ext_gcd(ll a,ll b,ll& x,ll& y){ ll d = a; if (!b){ x = 1;y = 0; }else{ d = ext_gcd(b,a%b,y,x); y -= a/b*x; } return d; } 在求逆元时，要找到 使得 ，实质上是求解方程 中的 ，因此可以用扩展欧几里得算法来求逆元，当然只有 时才有解，否则逆元不存在
2.4.4 推论 逆元的存在性 存在 的充要条件是
3. 同余方程组
有方程就有方程组
其中 是整数， 是正整数
3.1.1 中国剩余定理（孙子定理） 设上述方程组中 两两互质，则方程组的通解为
其中
下面是参考代码模板，需要调用前面的扩展欧几里得算法模板
ll Sunzi(ll *m,ll *a,int len){ ll lcm = 1; ll ans = 0; for (int i=0;i<len;i++){ ll k0,ki; ll d = ext_gcd(lcm,m[i],k0,ki); if ((a[i]-ans)%d!=0) return -1; else { ll t = m[i]/d; k0 = ( k0*(a[i]-ans)/d%t + t)%t; ans = k0*lcm + ans; lcm = lcm/d*m[i]; } } return ans; } 4. 素数
素数是只有1和本身两个因数的数，1不是素数
4.1 素数的判断
bool isPrime(ll n){ if(n==1)return false; for(ll i=2;i*i<=n;i++) if(n%i==0)return false; return true; } 这是最简单的素数的判断的参考代码模板，复杂度为
原理其实很简单，对于一个大于1的整数，如果 是它的一个大于 的因子，那么 是它的小于 的因子
在大多数情况下，这种判断方式的复杂度已经足够小了，如果要追求更高的效率，可以考虑 法
一个大于1的整数如果不是素数，那么一定有素因子，因此在枚举因子时只需要考虑可能为素数的因子即可。 法即枚举形如 的数，例如取 ，那么 都不可能为素数，只需要枚举形如 的数即可，这样复杂度降低了 。
下面的模板是 法 的版本
bool isPrime(ll n){ if(n==2||n==3||n==5)return 1; if(n%2==0||n%3==0||n%5==0||n==1)return 0; ll c=7,a[8]={4,2,4,2,4,6,2,6}; while(c*c<=n)for(auto i:a){if(n%c==0)return 0;c+=i;} return 1; } 如果 极大，可以使用素数测试算法，素数测试算法可以通过控制迭代次数来间接控制正确率，常用的有下面的Miller-Rabin方法
ll Rand(){ static ll x=(srand((int)time(0)),rand()); x+=1000003; if(x>1000000007)x-=1000000007; return x; } bool Witness(ll a,ll n){ ll t=0,u=n-1; while(!(u&1))u>>=1,t++; ll x=fpow(a,u,n),y; while(t--){ y=x*x%n; if(y==1 && x!=1 && x!=n-1)return true; x=y; } return x!=1; } bool MillerRabin(ll n,ll s){ if(n==2||n==3||n==5)return 1; if(n%2==0||n%3==0||n%5==0||n==1)return 0; while(s--){ if(Witness(Rand()%(n-1)+1,n))return false; } return true; } 当然， Rand() 怎么写可以自由发挥，这会影响其性能
4.2 素数筛
如果要求出不超过 的所有素数，素数筛是最好的选择，下面是一种朴素的筛法
void getPrime(bool p[],int n){ for(int i=1;i<=n;i++)p[i]=true; p[1]=false; for(int i=2;i<=n;i++){ if(p[i]){ for(int j=i+i;j<=n;j+=i)p[j]=false; } } } 这种方法的原理是从小到大将素数的倍数筛掉，复杂度为 ，注意到每个合数如果有多个素因子，那么就会被重复筛掉，造成复杂度的浪费，因此，用下面的方法可以保证每个合数只被它最小的素因子筛掉一遍，以 的复杂度解决上述问题
ll getPrime(ll n,bool vis[],ll prime[]){ ll tot=0; for(ll i=1;i<=n;i++)vis[i]=0; for(ll i=2;i<=n;i++){ if(!vis[i])prime[tot++]=i; for(ll j=0;j<tot;j++){ if(prime[j]*i>n)break; vis[prime[j]*i]=1; if(i%prime[j]==0)break; } } return tot; } 5. 组合数取余
5.1 组合数
5.1.1 定义 组合数 在 个不同元素中选取 个元素，不同的取法记为
5.2 杨辉三角
组合数与杨辉三角中的数字是一一对应的
杨辉三角的自然数形式
杨辉三角的组合数形式
按照上面的写法，杨辉三角的第 行第 列即为
注意到上图中每个数等于其左边的数与上边的数（如果有的话）之和，这就是杨辉恒等式
5.3 组合数取余的求法
在ACM竞赛中，我们常常需要计算 ，可以参考下面几种方法
- 如果 很小（不超过50），可以用C++的库函数
double tgamma(double x)，这是一个欧拉积分
在整数点处的取值满足
因此代码可以这么写
ll C(ll n,ll m){ return (ll)round(tgamma(n+1)/tgamma(m+1)/tgamma(n-m+1)); } 效率并不高，但是对于追求手速来说足够了
- 如果 不大，可以开 的空间，可以利用杨辉恒等式来预处理组合数表
const ll mo=1e9+7; ll C[1005][1005]; void getC(int n){ for(int i=0;i<=n;i++){ for(int j=0;j<=i;j++){ if(j==0 || j==i) C[i][j]=1; else C[i][j]=(C[i-1][j-1]+C[i-1][j])%mo; } } } - 如果 比较大，可以开 的空间，可以利用前文所述的逆元来求解，当然，要保证 是素数
const ll mo=1e9+7; ll C(ll n,ll m){ static ll M=0,inv[N],mul[N],invMul[N]; while(M<=n){ if(M){ inv[M]=M==1?1:(mo-mo/M)*inv[mo%M]%mo; mul[M]=mul[M-1]*M%mo; invMul[M]=invMul[M-1]*inv[M]%mo; } else mul[M]=1,invMul[M]=1; M++; } return mul[n]*invMul[m]%mo*invMul[n-m]%mo; } 上面的代码中用 的复杂度处理了 的逆元，处理 次 的询问的总复杂度为
- 如果 更大， 是素数，可以用Lucas定理来求解
5.3.1 Lucas定理
若 是素数，则
其中
即将 表示成 进制形式
5.3.2 推论
ll Lucas(ll n,ll m,ll p){ ll ans=1; while(n|m)ans=ans*C(n%P,m%P)%P,n/=P,m/=P; return ans; } 6. 与 的结构
6.1 的结构
6.1.1 算数基本定理 任何一个大于1的整数 ，都可以唯一地表示成素数乘积的形式
其中 是素数
对于一个较大的数，有用来分解其因数的 Pollard Rho 算法
typedef long long ll; typedef pair<int,int> pii; typedef pair<ll,ll> pll; namespace pollard_rho{ const int C=2307; const int S=10; typedef pair<ll,int> pli; mt19937 rd(time(0)); vector<ll>ve; ll gcd(ll a,ll b){return b?gcd(b,a%b):a;} ll mul(ll a,ll b,ll mod){return (__int128)a*b%mod;} ll power(ll a,ll b,ll mod){ ll res=1;a%=mod; while(b){ if(b&1)res=mul(res,a,mod); a=mul(a,a,mod); b>>=1; } return res; } bool check(ll a,ll n){ ll m=n-1,x,y; int j=0; while(!(m&1))m>>=1,j++; x=power(a,m,n); for(int i=1;i<=j;x=y,i++){ y=mul(x,x,n); if(y==1&&x!=1&&x!=n-1)return 1; } return y!=1; } bool miller_rabin(ll n){ ll a; if(n==1)return 0; if(n==2)return 1; if(!(n&1))return 0; for(int i=0;i<S;i++)if(check(rd()%(n-1)+1,n))return 0; return 1; } ll pollard_rho(ll n,int c){ ll i=1,k=2,x=rd()%n,y=x,d; while(1){ i++;x=(mul(x,x,n)+c)%n,d=gcd(y-x,n); if(d>1&&d<n)return d; if(y==x)return n; if(i==k)y=x,k<<=1; } } void findfac(ll n,int c){ if(n==1)return ; if(miller_rabin(n)){ ve.push_back(n); return ; } ll m=n; while(m==n)m=pollard_rho(n,c--); findfac(m,c); findfac(n/m,c); } vector<pli> solve(ll n){ vector<pli>res; ve.clear(); findfac(n,C); sort(ve.begin(),ve.end()); for(auto x:ve){ if(res.empty()||res.back().fi!=x)res.push_back({x,1}); else res.back().se++; } return res; } } 6.1.2 推论 若
则
6.1.3 定理 的阶乘中 的幂次为
6.1.4 定理
6.2
的结构
定理6.2.1 是循环群，即存在 ，使得
这样的 称为 的原根
素数一定有原根，原根不唯一，部分合数也有原根
原根一般不大，暴力枚举即可
1000000007的原根为5，998244353的原根为3
7. 离散对数与BSGS算法
考虑求解方程
这样的 称为离散对数，可以写为
Baby step giant step 算法
设 （ 为某常正整数），则原方程可以写成
将 存入表 (table，C++中可以用unordered_map) 中，然后枚举 ，在表中查找 即可，复杂度为 ，取 ，那么复杂度为
ll bsgs(ll a,ll b,ll p){ static unordered_map<ll,ll> tab; tab.clear(); ll u=(ll)sqrt(p)+1; ll now=1,step; rep(i,0,u-1){ ll tmp=b*inv(now,p)%p; if(!tab.count(tmp))tab[tmp]=i; (now*=a)%=p; } step=now; now=1; for(ll i=0;i<p;i+=u){ if(tab.count(now))return i+tab[now]; (now*=step)%=p; } return -1; } 8. 高次同余方程
考虑求解方程
先求 的原根 ，设 ， ，用BSGS算法求出 ，方程可写成
进而有
用扩展欧几里得算法求出 ，也就求出了
ll ModEquationSolve(ll a,ll y,ll p){ a%=p-1; ll g=primitiveRoot(p),t=bsgs(g,y,p),z,z_,d=ext_gcd(a,p-1,z,z_); if(t%d!=0)return -1; ll tmp=(p-1)/d; z=(t/d*z%tmp+tmp)%tmp; return fpow(g,z,p); } 9. 积性函数
9.1 积性函数
9.1.1 定义 积性函数 映射 如果满足：任意 且 ，有
那么称 为积性函数
最简单的积性函数 与
其中 为示性函数，当条件 成立时取1，否则取0
9.1.2 定理 积性函数必满足
9.2 欧拉函数
9.2.1 定义 欧拉函数
欧拉函数是积性函数，且有
vector<ll> factors(ll x){ vector<ll> fac; ll y=x; for(ll i=2;i*i<=x;i++){ if(y%i==0){ fac.push_back(i); while(y%i==0)y/=i; if(y==1)return fac; } } if(y!=1)fac.push_back(y); return fac; } ll Euler(ll n){ vector<ll> fac=factors(n); ll ans=n; for(auto p:fac)ans=ans/p*(p-1); return ans; } 9.2.2 欧拉定理 若 是大于1的整数， ，则
这又是一个求逆元的方法，其结果与扩展欧几里得算法相同
9.3 莫比乌斯函数
9.3.1 定义 莫比乌斯函数
ll Mobius(ll n){ vector<ll> fac=factors(n); for(auto p:fac)n/=p; return n>1?0:(fac.size()&1)?-1:1; } 9.3.2 定理 莫比乌斯反演
9.3.2' 定理 莫比乌斯反演的另一种形式
9.4 积性函数筛
积性函数可以通过前文中的线性素数筛在线性时间复杂度内求出，某些取值依赖于自变量的素因子的数论函数也能通过其求出
#define rep(i,a,b) for(ll i=a;i<=b;i++) const ll N=10005; //积性数论函数筛 void Number_Theory_Function_Solve(vector<ll> &a,ll n,string function_name){ static bool vis[N]; static ll prime[N]; a.resize(n+1); ll tot=0; rep(i,1,n)vis[i]=0; a[1]=1; if(function_name=="IsPrime"){ /* { 1, if x is a prime f(x) = { { 0, otherwise 素数判断函数不是积性函数，但可以用积性函数筛求解 */// rep(i,2,n){ if(!vis[i])prime[tot++]=i,a[i]=1; rep(j,0,tot-1){ if(prime[j]*i>n)break; vis[prime[j]*i]=1; if(i%prime[j]==0){ a[i*prime[j]]=0; break; } a[i*prime[j]]=0; } } a[1]=0; }else if(function_name=="E"){ /* { 1, x=1 f(x) = I(x=1) = { { 0, otherwise 积性函数群中的单位元 */// rep(i,2,n)a[i]=0; }else if(function_name=="Euler"){ /* n f(x) = sum I(gcd(i,x)=1) = x product (1-1/p) i=1 p|x, p is Prime 欧拉函数 */// rep(i,2,n){ if(!vis[i])prime[tot++]=i,a[i]=i-1; rep(j,0,tot-1){ if(prime[j]*i>n)break; vis[prime[j]*i]=1; if(i%prime[j]==0){ a[i*prime[j]]=a[i]*prime[j]; break; } a[i*prime[j]]=a[i]*a[prime[j]]; } } }else if(function_name=="Mobius"){ /* { 1, x=1 f(x) = { (-1)^k, n=p1 p2 ... pk, pi is prime { 0, otherwise 莫比乌斯函数 */// rep(i,2,n){ if(!vis[i])prime[tot++]=i,a[i]=-1; rep(j,0,tot-1){ if(prime[j]*i>n)break; vis[prime[j]*i]=1; if(i%prime[j]==0){ a[i*prime[j]]=0; break; } a[i*prime[j]]=a[i]*a[prime[j]]; } } }else if(function_name=="FactorCounter"){ /* n f(x) = sum I(i|x) i=1 因子个数 */// static ll t[N]; rep(i,2,n){ if(!vis[i])prime[tot++]=i,a[i]=2,t[i]=1; rep(j,0,tot-1){ if(prime[j]*i>n)break; vis[prime[j]*i]=1; if(i%prime[j]==0){ a[i*prime[j]]=a[i]/(t[i]+1)*(t[i]+2); t[i*prime[j]]=t[i]+1; break; } a[i*prime[j]]=a[i]*a[prime[j]]; t[i*prime[j]]=1; } } }else if(function_name=="FactorSum"){ /* n f(x) = sum I(i|x)*i i=1 因子和 */// static ll pw[N]; rep(i,2,n){ if(!vis[i])prime[tot++]=i,a[i]=i+1,pw[i]=i; rep(j,0,tot-1){ if(prime[j]*i>n)break; vis[prime[j]*i]=1; if(i%prime[j]==0){ a[i*prime[j]]=(pw[i]*prime[j]*prime[j]-1)/(prime[j]-1)*a[i/pw[i]]; pw[i*prime[j]]=pw[i]*prime[j]; break; } a[i*prime[j]]=a[i]*a[prime[j]]; pw[i*prime[j]]=prime[j]; } } }else if(function_name=="MinFactor"){ rep(i,2,n){ if(!vis[i])prime[tot++]=i,a[i]=i; rep(j,0,tot-1){ if(prime[j]*i>n)break; vis[prime[j]*i]=1; if(i%prime[j]==0){ a[i*prime[j]]=prime[j]; break; } a[i*prime[j]]=prime[j]; } } }else{ assert(false); } } 9.5 数论函数的狄利克雷乘积
9.5.1 定义 数论函数的狄利克雷乘积
9.5.2 定理 所有积性数论函数构成一个群，该群有单位元
任意积性数论函数 ，都存在唯一的积性数论函数 ，使得
且
称 为 的狄利克雷逆函数
9.5.3 推论 数论函数间的关系 设 为欧拉函数， 为莫比乌斯函数， 为单位数论函数， 为常数函数（ ）， 为恒等函数（ ），那么
bool mulFunCheck(const vector<ll> &val){ ll n=val.size()-1; static bool vis[N]; static ll prime[N],pw[N],t[N],tot=0; if(val[1]!=1)return false; rep(i,1,n)vis[i]=0; rep(i,2,n){ if(!vis[i])prime[tot++]=i,pw[i]=i,t[i]=1; rep(j,0,tot-1){ if(prime[j]*i>n)break; vis[prime[j]*i]=1; if(i%prime[j]==0){ if(val[i*prime[j]]!=val[i/pw[i]]*val[pw[i]*prime[j]])return false; pw[i*prime[j]]=pw[i]*prime[j],t[i*prime[j]]=t[i]+1; break; } if(val[i*prime[j]]!=val[i]*val[prime[j]])return false; pw[i*prime[j]]=prime[j],t[i*prime[j]]=1; } } return true; } //狄利克雷乘积 h(n)=sum_{d|n}f(d)g(n/d) vector<ll> DirichletProduct(const vector<ll> &f,const vector<ll> &g){ ll n=f.size()-1; vector<ll> h(n+1,0); rep(i,1,n)h[i]=0; rep(i,1,n)rep(j,1,n/i)h[i*j]+=f[i]*g[j]; return h; } //狄利克雷逆函数 f*g=g*f=[n==1] vector<ll> DirichletInversion(const vector<ll> &f){ ll n=f.size()-1; vector<ll> g(n+1,0); g[1]=1; rep(i,1,n)rep(j,2,n/i)g[i*j]-=f[j]*g[i]; return g; } //Mobius反演 f(n)=sum(g(d),d|n) g(n)=sum(f(d)mu(n/d),d|n) vector<ll> MobiusInversion(const vector<ll> &f){ ll n=f.size()-1; static vector<ll> mu; Number_Theory_Function_Solve(mu,n,"Mobius"); return DirichletProduct(mu,f); } 10. 多项式乘积算法
试想这样一个问题，求两个多项式
的乘积
使用传统的方法至少需要 的复杂度，如何加速上述过程？
10.1 快速傅里叶变换
首先考虑如何用其他方式表示多项式 .
任取n个不同的数(可以是整数、实数，甚至是复数)
将其代入 中，就得到一个线性方程组
只要 足够大，就能够唯一地确定一个多项式，换言之，上述方程组可以表示一个多项式，将这两种多项式的表示方法分别称为系数表示和点值表示.
利用快速傅里叶变换来求多项式乘积的总体思路是
- 选取合适的 个不同的数
- 将多项式 与 转化为点值表示(称为离散傅里叶变换，简称 )
- 计算 的点值表示
- 将 转化为系数表示(称为逆离散傅里叶变换，简称 )
1. 选取合适的 个不同的数
我们选取复数域上 的 个不同的值(或称 个 次单位复根)作为 的值，即
至于指数形式的复数 ，用大家所熟知的欧拉公式即可求得其代数形式
经过简单计算可知
当 为偶数时
其中 为 除以 的余数，上述两等式将在后文中使用.
我们为什么要费尽周折选取如此复杂的点呢?是为了使用快速傅里叶变换.
2. 将多项式 与 转化为点值表示
考虑多项式 ，当 时(当不满足该条件时，向 补充系数为0的高次项来扩大 使其满足该条件)，将其化为两个多项式
则有
进而
也就是说，要求 在 个不同点处的值，只需要求 与 在 个不同点处的值，由于 ，可对 与 重复进行上述过程，最终经过 步后得到 个函数
之后回推得到 的点值表示，上述过程就是快速傅里叶变换的过程，复杂度为 即 .
当然，还需要对 进行同样的变换.
3. 计算 的点值表示
点值表示的优点是可以快速地求出两个选取了相同点值的多项式的乘积，例如多项式
与多项式
的乘积
只需要 的复杂度即可求得.
4. 将 转化为系数表示
下面以 为例，讲解如何将多项式从点值表示转化为系数表示，此过程又称多项式的插值.
将 的点值表示写成矩阵形式 即
此处矩阵 中的1可视为 .
现在我们已知的是 和 ，要求的是 ， 是一范德蒙德矩阵，可求得其逆矩阵
因此 即
也就是说，只需将 与 对换，将 换成 ，再乘上系数 ，进行类似步骤2的变换，即可进行逆快速傅里叶变换，算法复杂度同样为 .
按照上述方法将 转化为系数表示，本题得解.
#include<bits/stdc++.h> using namespace std; const int maxn=1<<18; //typedef complex<double> C; struct C { double a,b; C(){} C(double a,double b):a(a),b(b){} C operator = (int a){*this=C(a*1.0,0);return *this;} C operator + (const C &t){return C(a+t.a,b+t.b);} C operator - (const C &t){return C(a-t.a,b-t.b);} C operator * (const C &t){return C(a*t.a-b*t.b,a*t.b+b*t.a);} }; C wn(int n,int f) { return C(cos(acos(-1.0)/n),f*sin((acos(-1.0))/n)); } C inv(int n) { return C(1.0/n,0); } C a[maxn],b[maxn],c[maxn]; int g[maxn]; void FFT(C *a,int n,int f) { for(int i=0;i<n;i++)if(i>g[i])swap(a[i],a[g[i]]); for(int i=1;i<n;i<<=1) { C w=wn(i,f),x,y; for(int j=0;j<n;j+=i+i) { C e;e=1; for(int k=0;k<i;e=e*w,k++) { x=a[j+k]; y=a[j+k+i]*e; a[j+k]=x+y; a[j+k+i]=x-y; } } } if(f==-1) { C Inv=inv(n); for(int i=0;i<n;i++)a[i]=a[i]*Inv; } } void conv(C *a,int n,C *b,int m,C *c) { int k=0,s=2; while((1<<k)<max(n,m)+1)k++,s<<=1; for(int i=1;i<s;i++)g[i]=(g[i/2]/2)|((i&1)<<k); FFT(a,s,1); FFT(b,s,1); for(int i=0;i<s;i++)c[i]=a[i]*b[i]; FFT(c,s,-1); } int main() { int n,m; scanf("%d%d",&n,&m); for(int i=0;i<=n;i++)scanf("%lf",&a[i]); for(int i=0;i<=m;i++)scanf("%lf",&b[i]); conv(a,n,b,m,c); for(int i=0;i<=n+m;i++)printf("%d ",(int)(c[i].a+0.5)); return 0; } 10.2 快速数论变换
快速傅里叶变换可以用于计算两个实系数多项式的乘积，但毕竟有精度问题，对于整系数多项式和模 意义下的乘法，是否有无损精度的算法呢？
FFT中精度损失的关键出在第一步—— 的选取，算法的实现依赖一个基本的等式
这个等式导致， 有 个不同的取值，但它们的平方只有 个不同的取值，这 个值的平方只有 个不同的取值，以此类推
若 满足某些条件时，是可以找到整数 的，例如当 ， 时，取
此时恰好满足类似的性质，后续的算法步骤类似，但这种变换就不是离散傅里叶变换（DFT）了，而是数论变换（NTT）
const ll mo=998244353; ll fpow(ll a, ll b){ ll ans=1; while(b>0){if(b&1)ans=ans*a%mo;b>>=1;a=a*a%mo;} return ans; } ll D(ll x) { ((x>=mo) && (x-=mo)) || ((x<0) && (x+=mo)); return x; } void NTT(ll a[],ll n,ll op) { for(ll i=1,j=n>>1;i<n-1;++i) { if(i<j) swap(a[i],a[j]); ll k=n>>1; while(k<=j) { j-=k; k>>=1; } j+=k; } for(ll len=2;len<=n;len<<=1) { ll rt=fpow(3,(mo-1)/len); for(ll i=0;i<n;i+=len) { ll w=1; for(ll j=i;j<i+len/2;++j) { ll u=a[j],t=1LL*a[j+len/2]*w%mo; a[j]=D(u+t),a[j+len/2]=D(u-t); w=1LL*w*rt%mo; } } } if(op==-1) { reverse(a+1,a+n); ll in=fpow(n,mo-2); for(ll i=0;i<n;++i) a[i]=1LL*a[i]*in%mo; } } vector<ll> Conv(vector<ll> const &A,vector<ll> const &B,ll N) { static ll a[2000005],b[2000005]; auto Make2=[](ll x)->ll { return 1<<((32-__builtin_clz(x))+((x&(-x))!=x)); }; ll n=Make2(A.size()+B.size()-1); for(ll i=0;i<n;++i) { a[i]=i<A.size()?A[i]:0; b[i]=i<B.size()?B[i]:0; } NTT(a,n,1);NTT(b,n,1); for(ll i=0;i<n;++i) a[i]=1LL*a[i]*b[i]%mo; NTT(a,n,-1); vector<ll> C(N); for (ll i=0;i<N;i++) C[i]=a[i]; return C; } 杂技
1. 积性函数前缀和
1.1 min25筛
用于求积性函数 的前 项和，要求 在素数处为一多项式，在素数的整数幂处可以直接求值。
#define rep(i,a,b) for(ll i=a;i<=b;i++) typedef long long ll; class MultiplicativeFunction{ private: static const ll SQRTN=1100000; ll prime[SQRTN],pcnt,cnt; bool vis[SQRTN]; ll d,n; ll val[SQRTN*2],id1[SQRTN],id2[SQRTN]; ll G[SQRTN*2],fpval[SQRTN*2]; ll g[SQRTN*2]; function<ll(ll,ll)> f_prime_pow; ll sum_pow_n_k(ll n,ll k){ if(k==0)return n; if(k==1)return n*(n+1)/2; if(k==2)return n*(n+1)*(2*n+1)/6; return -1; } void getPrime(ll n){ pcnt=0; fill(vis+1,vis+1+d,0); rep(i,2,n){ if(!vis[i])prime[++pcnt]=i; rep(j,1,pcnt){ if(prime[j]*i>n)break; vis[prime[j]*i]=1; if(i%prime[j]==0)break; } } } ll powll(ll x,ll k){ ll ans=1; while(k--)ans*=x; return ans; } ll &index(ll x){return x<=d?id1[x]:id2[n/x];} void sieve(ll x,ll k){ rep(i,1,pcnt)fpval[i]=powll(prime[i],k)+fpval[i-1]; rep(i,1,cnt)g[i]=sum_pow_n_k(val[i],k)-1; rep(i,1,pcnt){ if(prime[i]*prime[i]>n)break; rep(j,1,cnt){ if(prime[i]*prime[i]>val[j])break; g[j]=g[j]-powll(prime[i],k)*(g[index(val[j]/prime[i])]-fpval[i-1]); } } } ll S(ll x,ll now){ if(x<=1 || x<prime[now])return 0; ll ans=G[index(x)]-fpval[now-1]; rep(i,now,pcnt){ if(prime[i]*prime[i]>x)break; ll prod=prime[i]; for(ll k=1;prod*prime[i]<=x;k++,prod*=prime[i]){ ans+=f_prime_pow(prime[i],k)*S(x/prod,i+1)+f_prime_pow(prime[i],k+1); } } return ans; } public: ll solve(ll x,const vector<ll> &c,function<ll(ll,ll)> fp){ f_prime_pow=fp; n=x;d=sqrt(n);cnt=0; getPrime(d); for(ll i=1;i<=n;i=n/(n/i)+1){ val[++cnt]=n/i; index(val[cnt])=cnt; } fill(G+1,G+1+cnt,0); for(int k=0;k<(int)c.size();k++){ if(c[k]!=0){ sieve(x,k); rep(i,1,cnt)G[i]+=c[k]*g[i]; } } rep(i,1,pcnt){ fpval[i]=0; for(ll k=(ll)c.size()-1;k>=0;k--){ fpval[i]=fpval[i]*prime[i]+c[k]; } fpval[i]+=fpval[i-1]; } return S(n,1)+1; } }; 使用说明：调用 solve 函数，c 是函数 在素数处的解析式（多项式）的系数，fp 是 在素数的整数幂处的取值计算函数。
例：
MultiplicativeFunction f; //Euler函数 cout<<f.solve(n,{-1,1},[](ll p,ll k){ll a=1;while(k--)a*=p;return a/p*(p-1);})<<endl; //Mobius函数 cout<<f.solve(n,{-1},[](ll p,ll k){return k==1?-1:0;})<<endl; 1.2 另一种筛
用于求积性函数 的前 项和，要求找一个函数 ，使得 与 的狄利克雷乘积、 可以快速求前 项和。
class sieve{ private: static const ll SQRTN=100005; ll n,d,cnt; ll val[SQRTN*2],id1[SQRTN],id2[SQRTN],s[SQRTN*2]; vector<ll> sum; bool vis[SQRTN*2]; ll &index(ll x){return x<=d?id1[x]:id2[n/x];} function<ll(ll)> sum_g,sum_f_g; ll S(ll x){ if(x<sum.size())return sum[x]; if(x==1)return 1; ll p=index(x); if(vis[p])return val[p]; val[p]=sum_f_g(x); for(ll i=2;i<=x;i=x/(x/i)+1){ val[p]-=(sum_g(x/(x/i))-sum_g(i-1))*S(x/i); } vis[p]=1; return val[p]; } public: void init(function<ll(ll)> sg,function<ll(ll)> sfg,const vector<ll> &v={}){sum_g=sg;sum_f_g=sfg;sum=v;} ll operator () (ll x){ n=x; d=sqrt(x); cnt=0; for(ll i=1;i<=x;i=x/(x/i)+1){ val[++cnt]=x/i; vis[cnt]=0; index(val[cnt])=cnt; } return S(n); } }; 使用说明：调用 init 函数进行初始化，sg 为 的前 项和，sfg 为 的前 项和，如果需要，可以预处理前若干项的和放到 v 中进行剪枝；直接使用 () 进行求和。
例：
sieve f; //Mobius函数 mu(x) g=1 f*g=E f.init([](ll x){return x;},[](ll x){return x>=1;}); cout<<f(n)<<endl; //Euler函数 phi(x) g=1 f*g=id f.init([](ll x){return x;},[](ll x){return x*(x+1)/2;}); cout<<f(n)<<endl; //x^k phi(x) g=id^k f*g=id^(k+1) f.init([](ll x){return x*(x+1)/2;},[](ll x){return x*(x+1)*(2*x+1)/6;}); cout<<f(n)<<endl; //x^k mu(x) g=id^k f*g=E f.init([](ll x){return x*(x+1)/2;},[](ll x){return x>=1;}); cout<<f(n)<<endl; 一个简化版本，使用方式相同：
class sieve2{ private: vector<ll> sum; unordered_map<ll,ll> mp; function<ll(ll)> sum_g,sum_f_g; ll S(ll x){ if(x<sum.size())return sum[x]; if(mp.find(x)!=mp.end())return mp[x]; ll ans=sum_f_g(x); for(ll i=2;i<=x;i=x/(x/i)+1){ ans-=(sum_g(x/(x/i))-sum_g(i-1))*S(x/i); } mp[x]=ans; return ans; } public: void init(function<ll(ll)> sg,function<ll(ll)> sfg,const vector<ll> &v={0,1}){ mp.clear(); sum_g=sg;sum_f_g=sfg;sum=v; } ll operator () (ll x){ return S(x); } }; 2. 二次剩余
Shanks 算法，求解二次同余方程 ， 是素数。
ll fpow(ll a,ll b,ll c){a%=c;ll ans=1;while(b>0){if(b&1)ans=ans*a%c;b>>=1;a=a*a%c;}return ans;} bool IsQuadraticResidue(ll a,ll p){//判断a是否为p的二次剩余 p是素数 return fpow(a,(p-1)>>1,p)==1; } ll Shanks(ll a,ll p){//求解二次同余方程x^2=a(mod p) p是素数 if(a==0)return 0; ll q=p-1,e=0; while(!(q&1))q>>=1,e++; static mt19937_64 rd(time(0)); ll n=rd()%(p-1)+1;//随机选取p的一个非二次剩余，若p为定值，n也可为定值 while(IsQuadraticResidue(n,p))n=rd()%(p-1)+1; ll z=fpow(n,q,p),y=z,r=e,x=fpow(a,(q-1)>>1,p),b=a*x%p*x%p; x=a*x%p; while(b!=1){ ll temp=b*b%p,m=1; while(temp!=1)(temp*=temp)%=p,m++; if(m==r)return -1; ll t=y; rep(i,1,r-m-1)(t*=t)%=p; y=t*t%p,r=m%p,x=x*t%p,b=b*y%p; } return x; } 3. 组合数
在 Lucas 定理的基础上修改，可以求 除以 的余数， 可以非素数。
class Combine{ private: ll fpow(ll a,ll b,ll c){a%=c;ll ans=1;while(b>0){if(b&1)ans=ans*a%c;b>>=1;a=a*a%c;}return ans;} void factorization(ll x,vector<ll> &p,vector<ll> &w,vector<ll> &pw){ p.clear(),w.clear(),pw.clear(); ll y=x; for(ll i=2;i*i<=x;i++)if(y%i==0){ p.push_back(i),w.push_back(0),pw.push_back(1); while(y%i==0)y/=i,w.back()++,pw.back()*=i; } if(y!=1)p.push_back(y),w.push_back(1),pw.push_back(y); } ll powll(ll a,ll b){ll ans=1;while(b)ans*=a,b--;return ans;} ll ext_gcd(ll a,ll b,ll& x,ll& y){ ll d=a; if(!b)x=1,y=0; else d=ext_gcd(b,a%b,y,x),y-=a/b*x; return d; } ll inv(ll a,ll b){ ll x,y,d=ext_gcd(a,b,x,y); assert(d==1); return (x%b+b)%b; } ll crt(const vector<ll> &m,const vector<ll> &b){ ll ans=0,M=1; for(auto i:m)M*=i; for(int i=0;i<(int)m.size();i++){ (ans+=M/m[i]*inv(M/m[i],m[i])%M*b[i])%=M; } return ans; } public: ll mod; vector<ll> fac[105],p,w,pw,res; void init(ll mo){ mod=mo; factorization(mod,p,w,pw); res.resize(p.size()); for(int i=0;i<(int)p.size();i++){ fac[i].resize(pw[i]); fac[i][0]=1; for(int j=1;j<pw[i];j++){ fac[i][j]=(j%p[i]==0)?fac[i][j-1]:(fac[i][j-1]*j%pw[i]); } } } Combine(){} Combine(ll mod):mod(mod){init(mod);} ll factorial(ll n,ll i){ if(n<p[i])return fac[i][n]; else return fpow(fac[i][pw[i]-1],n/pw[i],pw[i])*fac[i][n%pw[i]]%pw[i]*factorial(n/p[i],i)%pw[i]; } ll C_pw(ll n,ll m,ll i){ ll c=0; for(ll j=p[i];j<=n;j*=p[i])c+=n/j; for(ll j=p[i];j<=m;j*=p[i])c-=m/j; for(ll j=p[i];j<=n-m;j*=p[i])c-=(n-m)/j; return factorial(n,i)*inv(factorial(m,i),pw[i])%pw[i]*inv(factorial(n-m,i),pw[i])%pw[i]*fpow(p[i],c,pw[i])%pw[i]; } ll operator () (ll n,ll m){ assert(m<=n && m>=0); for(int i=0;i<(int)p.size();i++)res[i]=C_pw(n,m,i); return crt(pw,res); } }C; 4. 一个求和技巧（常被称为类欧几里得算法）
求 ， 表示取整。
ll sum_pow(ll n,ll k){ if(k==0)return n; else if(k==1)return n*(n+1)/2; else if(k==2)return n*(n+1)*(2*n+1)/6; else if(k==3)return n*n*(n+1)*(n+1)/4; else if(k==4)return n*(2*n+1)*(n+1)*(3*n*n+3*n-1)/30; else assert(false); } ll EuclidLike1(ll a,ll b,ll c,ll n){ if(a==0)return b/c*(n+1); else if(a>=c || b>=c)return (a/c)*sum_pow(n,1)+(b/c)*(n+1)+EuclidLike1(a%c,b%c,c,n); else return (a*n+b)/c*n-EuclidLike1(c,c-b-1,a,(a*n+b)/c-1); } 5. 模数任意的 FFT
模数任意（可取 甚至非素数）的 FFT，你可以忘掉 NTT 了。
模板需要 define 相应的关键字（FFT_conv,FFT_conv_mod,FFT_conv_mod_fast）来启用相关功能。
namespace FFT{ // FFT_conv:普通多项式乘积 // FFT_conv_mod:多项式乘积取模 // FFT_conv_mod_fast:优化的多项式乘积取模 #define FFT_conv_mod_fast const int maxn=1<<20; //typedef complex<double> cp; struct cp{ long double a,b; cp(){} cp(long double a,long double b):a(a),b(b){} cp(ll a):a(a*1.0),b(0){} cp operator + (const cp &t)const{return cp(a+t.a,b+t.b);} cp operator - (const cp &t)const{return cp(a-t.a,b-t.b);} cp operator * (const cp &t)const{return cp(a*t.a-b*t.b,a*t.b+b*t.a);} cp conj()const{return cp(a,-b);} }; cp wn(int n,int f){ const long double pi=acos(-1.0); return cp(cos(pi/n),f*sin(pi/n)); } int g[maxn]; void DFT(cp *a,int n,int f){ for(int i=0;i<n;i++)if(i>g[i])swap(a[i],a[g[i]]); for(int i=1;i<n;i<<=1){ cp w=wn(i,f),x,y; for(int j=0;j<n;j+=i+i){ cp e(1,0); for(int k=0;k<i;e=e*w,k++){ x=a[j+k]; y=a[j+k+i]*e; a[j+k]=x+y; a[j+k+i]=x-y; } } } if(f==-1){ cp Inv(1.0/n,0); for(int i=0;i<n;i++)a[i]=a[i]*Inv; } } #ifdef FFT_conv cp a[maxn],b[maxn],c[maxn]; vector<ll> conv(const vector<ll> &u,const vector<ll> &v,ll sz){ int k=0,s=1,n=(int)u.size()-1,m=(int)v.size()-1; while(s<n+m+1)k++,s<<=1; for(int i=1;i<s;i++)g[i]=(g[i/2]/2)|((i&1)<<(k-1)); for(int i=0;i<s;i++)a[i]=i<=n?u[i]:0; for(int i=0;i<s;i++)b[i]=i<=m?v[i]:0; DFT(a,s,1); DFT(b,s,1); for(int i=0;i<s;i++)c[i]=a[i]*b[i]; DFT(c,s,-1); vector<ll> ans(sz+1,0); for(int i=0;i<=sz;i++)ans[i]=i<s?ll(round(c[i].a)):0; return ans; } #endif #ifdef FFT_conv_mod cp a1[maxn],a2[maxn],b1[maxn],b2[maxn],c[maxn]; vector<ll> conv_mod(const vector<ll> &u,const vector<ll> &v,int sz,ll mod){ int k=0,s=1,n=(int)u.size()-1,m=(int)v.size()-1,M=sqrt(mod)+1; while(s<n+m+1)k++,s<<=1; for(int i=1;i<s;i++)g[i]=(g[i/2]/2)|((i&1)<<(k-1)); for(int i=0;i<s;i++){ if(i<=n)a1[i]=u[i]%mod/M,b1[i]=u[i]%mod%M; else a1[i]=b1[i]=0; if(i<=m)a2[i]=v[i]%mod/M,b2[i]=v[i]%mod%M; else a2[i]=b2[i]=0; } DFT(a1,s,1); DFT(b1,s,1); DFT(a2,s,1); DFT(b2,s,1); vector<ll> ans(sz+1,0); for(int i=0;i<s;i++)c[i]=a1[i]*a2[i]; DFT(c,s,-1); for(int i=0;i<min(sz+1,s);i++)(ans[i]+=ll(round(c[i].a))%mod*M%mod*M)%=mod; for(int i=0;i<s;i++)c[i]=a1[i]*b2[i]+a2[i]*b1[i]; DFT(c,s,-1); for(int i=0;i<min(sz+1,s);i++)(ans[i]+=ll(round(c[i].a))%mod*M)%=mod; for(int i=0;i<s;i++)c[i]=b1[i]*b2[i]; DFT(c,s,-1); for(int i=0;i<min(sz+1,s);i++)(ans[i]+=ll(round(c[i].a))%mod)%=mod; return ans; } #endif #ifdef FFT_conv_mod_fast cp a[maxn],b[maxn],Aa[maxn],Ab[maxn],Ba[maxn],Bb[maxn]; vector<ll> conv_mod_fast(const vector<ll> &u,const vector<ll> &v,int sz,ll mod){ int k=0,s=1,n=(int)u.size()-1,m=(int)v.size()-1,M=sqrt(mod)+1; while(s<n+m+1)k++,s<<=1; for(int i=1;i<s;i++)g[i]=(g[i/2]/2)|((i&1)<<(k-1)); for(int i=0;i<s;i++){ if(i<=n)a[i].a=u[i]%mod%M,a[i].b=u[i]%mod/M; else a[i]=0; if(i<=m)b[i].a=v[i]%mod%M,b[i].b=v[i]%mod/M; else b[i]=0; } DFT(a,s,1); DFT(b,s,1); for(int i=0;i<s;i++){ int j=(s-i)%s; cp t1=(a[i]+a[j].conj())*cp(0.5,0); cp t2=(a[i]-a[j].conj())*cp(0,-0.5); cp t3=(b[i]+b[j].conj())*cp(0.5,0); cp t4=(b[i]-b[j].conj())*cp(0,-0.5); Aa[i]=t1*t3,Ab[i]=t1*t4,Ba[i]=t2*t3,Bb[i]=t2*t4; } for(int i=0;i<s;i++){ a[i]=Aa[i]+Ab[i]*cp(0,1); b[i]=Ba[i]+Bb[i]*cp(0,1); } DFT(a,s,-1); DFT(b,s,-1); vector<ll> ans(sz+1,0); for(int i=0;i<min(s,sz+1);i++){ ll t1=(ll)round(a[i].a)%mod; ll t2=(ll)round(a[i].b)%mod; ll t3=(ll)round(b[i].a)%mod; ll t4=(ll)round(b[i].b)%mod; (ans[i]+=t1+(t2+t3)*M%mod+t4*M*M)%=mod; } return ans; } #endif } 6. 多项式插值
6.1 Lagrange 插值
#define rep(i,a,b) for(ll i=a;i<=b;i++) const ll mo=1e9+7; ll fpow(ll a,ll b){ll ans=1;while(b>0){if(b&1)ans=ans*a%mo;b>>=1;a=a*a%mo;}return ans;} ll lagrange(vector<ll> x,vector<ll> y,ll X){ auto p=y.begin(); ll ans=0; for(auto k:x){ ll a=*p++%mo,b=1; for(auto j:x)if(j!=k)a=(X-j)%mo*a%mo,b=(k-j)%mo*b%mo; ans=(ans+mo+a*fpow(b,mo-2)%mo)%mo; } return ans; } 插值点连续的情况下， ，时间复杂度优化到 ， 是插值点的个数。
#define rep(i,a,b) for(ll i=a;i<=b;i++) ll Lagrange(const vector<ll> &y,ll x,ll mod){ ll n=y.size()-1; static const ll N=100005; static ll lx[N],rx[N],inv[N],mul[N]; rep(i,1,n)inv[i]=(i==1)?1:((mod-mod/i)*inv[mod%i]%mod); rep(i,0,n)lx[i]=rx[i]=(x-1-i)%mod; rep(i,1,n)(lx[i]*=lx[i-1])%=mod; for(ll i=n-1;i>=0;i--)(rx[i]*=rx[i+1])%=mod; rep(i,0,n)mul[i]=(i==0)?1:(mul[i-1]*inv[i])%mod; ll ans=0; rep(i,0,n){ ll now=1; if(i>0)(now*=lx[i-1]*mul[i]%mod)%=mod; if(i<n)(now*=rx[i+1]*mul[n-i]%mod)%=mod; if((n-i)&1)now=-now; (ans+=y[i]*now)%=mod; } return ans+(ans<0)*mod; } 二维的 Lagrange 插值
#define rep(i,a,b) for(ll i=a;i<=b;i++) type lagrange2(vector<type> x,vector<type> y,vector<vector<type> > z,type X,type Y){ int M=x.size()-1,N=y.size()-1; type ans=0; rep(m,0,M)rep(n,0,N){ type t=z[m][n]; rep(i,0,M)if(i!=m)t*=(X-x[i])/(x[m]-x[i]); rep(i,0,N)if(i!=n)t*=(Y-y[i])/(y[n]-y[i]); ans+=t; } return ans; } 6.2 牛顿插值
#define rep(i,a,b) for(ll i=a;i<=b;i++) ll fpow(ll a,ll b){ll ans=1;while(b>0){if(b&1)ans=ans*a%mo;b>>=1;a=a*a%mo;}return ans;} class NewtonPoly{ public: ll f[105],d[105],x[105],n=0; void add(ll X,ll Y){ x[n]=X,f[n]=Y%mo; rep(i,1,n)f[n-i]=(f[n-i+1]-f[n-i])%mo*fpow((x[n]-x[n-i])%mo,mo-2)%mo; d[n++]=f[0]; } ll cal(ll X){ ll ans=0,t=1; rep(i,0,n-1)ans=(ans+d[i]*t)%mo,t=(X-x[i])%mo*t%mo; return ans+mo*(ans<0); } }; 6.3 利用插值求 前 项和
#define rep(i,a,b) for(ll i=a;i<=b;i++) const ll mo=1e9+7; ll sum_pow_n_k(ll n,ll k){ static ll inv[1005],d[1005][1005],init_flag=1; if(init_flag){ init_flag=0; const ll all=1002; inv[1]=1; rep(i,2,all)inv[i]=(mo-mo/i)*inv[mo%i]%mo; d[1][1]=1; rep(i,2,all)rep(j,1,i)d[i][j]=(j*j*d[i-1][j]+(j-1)*d[i-1][j-1])%mo*inv[j]%mo; } n%=mo; ll ans=0,now=n; rep(i,1,k+1){ (ans+=d[k+1][i]*now)%=mo; (now*=(n-i))%=mo; } return ans+(ans<0)*mo; } 7. 线性递推数列
7.1 BM 算法
根据数列前若干项求递推公式。
例如斐波那契数列，传入 {1,1,2,3,5,8,13}，传出 {1,-1,-1}（ ）。
#define rep(i,a,b) for(ll i=a;i<=b;i++) const ll N=25; const ll mo=1e9+7; ll fpow(ll a,ll b){ll ans=1;while(b>0){if(b&1)ans=ans*a%mo;b>>=1;a=a*a%mo;}return ans;} //BM算法 求线性递推数列的递推公式 vector<ll> BM(const vector<ll> &s){ vector<ll> C={1},B={1},T; ll L=0,m=1,b=1; rep(n,0,s.size()-1){ ll d=0; rep(i,0,L)d=(d+s[n-i]%mo*C[i])%mo; if(d==0)m++; else{ T=C; ll t=mo-fpow(b,mo-2)*d%mo; while(C.size()<B.size()+m)C.push_back(0); rep(i,0,B.size()-1)C[i+m]=(C[i+m]+t*B[i])%mo; if(2*L>n)m++; else L=n+1-L,B=T,b=d,m=1; } } return C; } 7.2 线性递推数列求单项值
利用多项式求 ，复杂度为 ， 为递推公式长度。
s 为前若干项，C 为递推公式系数， 从 开始。
#define rep(i,a,b) for(ll i=a;i<=b;i++) ll polySolve(const vector<ll> &s,const vector<ll> &C,ll n){ if(n<s.size())return s[n]; static ll g[N],f[N],d[N]; ll k=(ll)C.size()-1,w=1; rep(i,0,k)f[i]=i==1,d[i]=i==k?1:C[k-i]; while((w<<1)<=n)w<<=1; while(w>>=1){ rep(i,0,k+k-2)g[i]=0; rep(i,0,k-1)if(f[i])rep(j,0,k-1)(g[i+j]+=f[i]*f[j])%=mo; for(ll i=k+k-2;i>=k;i--)if(g[i])rep(j,1,k)(g[i-j]-=g[i]*d[k-j])%=mo; rep(i,0,k-1)f[i]=g[i]; if(w&n)for(ll i=k;i>=0;i--)f[i]=i==k?f[i-1]:(i==0?-f[k]*d[i]:(f[i-1]-f[k]*d[i]))%mo; } ll ans=0; rep(i,0,k-1)(ans+=f[i]*s[i])%=mo; return ans+(ans<0)*mo; } 7.3 RS 算法
当模数非素数时，使用 RS 算法替代 BM 算法。
模板是网上搜刮来的，它终结了线性递推题目在ACM中的历史，模板作者留了一个 FFT 的优化接口，供下一个毒瘤出题人使用。
struct LinearSequence { typedef vector<long long> vec; static void extand(vec &a, size_t d, ll value = 0) { if (d <= a.size()) return; a.resize(d, value); } static vec BerlekampMassey(const vec &s, ll mod) { function<ll(ll)> inverse = [&](ll a) { return a == 1 ? 1 : (ll) (mod - mod / a) * inverse(mod % a) % mod; }; vec A = {1}, B = {1}; ll b = s[0]; for (size_t i = 1, m = 1; i < s.size(); ++i, m++) { ll d = 0; for (size_t j = 0; j < A.size(); ++j) { d += A[j] * s[i - j] % mod; } if (!(d %= mod)) continue; if (2 * (A.size() - 1) <= i) { auto temp = A; extand(A, B.size() + m); ll coef = d * inverse(b) % mod; for (size_t j = 0; j < B.size(); ++j) { A[j + m] -= coef * B[j] % mod; if (A[j + m] < 0) A[j + m] += mod; } B = temp, b = d, m = 0; } else { extand(A, B.size() + m); ll coef = d * inverse(b) % mod; for (size_t j = 0; j < B.size(); ++j) { A[j + m] -= coef * B[j] % mod; if (A[j + m] < 0) A[j + m] += mod; } } } return A; } static void exgcd(ll a, ll b, ll &g, ll &x, ll &y) { if (!b) x = 1, y = 0, g = a; else { exgcd(b, a % b, g, y, x); y -= x * (a / b); } } static ll crt(const vec &c, const vec &m) { int n = c.size(); ll M = 1, ans = 0; for (int i = 0; i < n; ++i) M *= m[i]; for (int i = 0; i < n; ++i) { ll x, y, g, tm = M / m[i]; exgcd(tm, m[i], g, x, y); ans = (ans + tm * x * c[i] % M) % M; } return (ans + M) % M; } static vec ReedsSloane(const vec &s, ll mod) { auto inverse = [](ll a, ll m) { ll d, x, y; exgcd(a, m, d, x, y); return d == 1 ? (x % m + m) % m : -1; }; auto L = [](const vec & a, const vec & b) { int da = (a.size() > 1 || (a.size() == 1 && a[0])) ? a.size() - 1 : -1000; int db = (b.size() > 1 || (b.size() == 1 && b[0])) ? b.size() - 1 : -1000; return max(da, db + 1); }; auto prime_power = [&](const vec & s, ll mod, ll p, ll e) { // linear feedback shift register mod p^e, p is prime vector<vec> a(e), b(e), an(e), bn(e), ao(e), bo(e); vec t(e), u(e), r(e), to(e, 1), uo(e), pw(e + 1);; pw[0] = 1; for (int i = pw[0] = 1; i <= e; ++i) pw[i] = pw[i - 1] * p; for (ll i = 0; i < e; ++i) { a[i] = {pw[i]}, an[i] = {pw[i]}; b[i] = {0}, bn[i] = {s[0] * pw[i] % mod}; t[i] = s[0] * pw[i] % mod; if (t[i] == 0) { t[i] = 1, u[i] = e; } else { for (u[i] = 0; t[i] % p == 0; t[i] /= p, ++u[i]); } } for (size_t k = 1; k < s.size(); ++k) { for (int g = 0; g < e; ++g) { if (L(an[g], bn[g]) > L(a[g], b[g])) { ao[g] = a[e - 1 - u[g]]; bo[g] = b[e - 1 - u[g]]; to[g] = t[e - 1 - u[g]]; uo[g] = u[e - 1 - u[g]]; r[g] = k - 1; } } a = an, b = bn; for (int o = 0; o < e; ++o) { ll d = 0; for (size_t i = 0; i < a[o].size() && i <= k; ++i) { d = (d + a[o][i] * s[k - i]) % mod; } if (d == 0) { t[o] = 1, u[o] = e; } else { for (u[o] = 0, t[o] = d; t[o] % p == 0; t[o] /= p, ++u[o]); int g = e - 1 - u[o]; if (L(a[g], b[g]) == 0) { extand(bn[o], k + 1); bn[o][k] = (bn[o][k] + d) % mod; } else { ll coef = t[o] * inverse(to[g], mod) % mod * pw[u[o] - uo[g]] % mod; int m = k - r[g]; extand(an[o], ao[g].size() + m); extand(bn[o], bo[g].size() + m); for (size_t i = 0; i < ao[g].size(); ++i) { an[o][i + m] -= coef * ao[g][i] % mod; if (an[o][i + m] < 0) an[o][i + m] += mod; } while (an[o].size() && an[o].back() == 0) an[o].pop_back(); for (size_t i = 0; i < bo[g].size(); ++i) { bn[o][i + m] -= coef * bo[g][i] % mod; if (bn[o][i + m] < 0) bn[o][i + m] -= mod; } while (bn[o].size() && bn[o].back() == 0) bn[o].pop_back(); } } } } return make_pair(an[0], bn[0]); }; vector<tuple<ll, ll, int>> fac; for (ll i = 2; i * i <= mod; ++i) if (mod % i == 0) { ll cnt = 0, pw = 1; while (mod % i == 0) mod /= i, ++cnt, pw *= i; fac.emplace_back(pw, i, cnt); } if (mod > 1) fac.emplace_back(mod, mod, 1); vector<vec> as; size_t n = 0; for (auto && x : fac) { ll mod, p, e; vec a, b; tie(mod, p, e) = x; auto ss = s; for (auto && x : ss) x %= mod; tie(a, b) = prime_power(ss, mod, p, e); as.emplace_back(a); n = max(n, a.size()); } vec a(n), c(as.size()), m(as.size()); for (size_t i = 0; i < n; ++i) { for (size_t j = 0; j < as.size(); ++j) { m[j] = get<0>(fac[j]); c[j] = i < as[j].size() ? as[j][i] : 0; } a[i] = crt(c, m); } return a; } LinearSequence(const vec &s, const vec &c, ll mod) : init(s), trans(c), mod(mod), m(s.size()) {} LinearSequence(const vec &s, ll mod, bool is_prime = true) : mod(mod) { vec A; if (is_prime) A = BerlekampMassey(s, mod); else A = ReedsSloane(s, mod); if (A.empty()) A = {0}; m = A.size() - 1; trans.resize(m); for (int i = 0; i < m; ++i) { trans[i] = (mod - A[i + 1]) % mod; } reverse(trans.begin(), trans.end()); init = {s.begin(), s.begin() + m}; } ll calc(ll n) { if (mod == 1) return 0; if (n < m) return init[n]; vec v(m), u(m << 1); int msk = !!n; for (ll m = n; m > 1; m >>= 1) msk <<= 1; v[0] = 1 % mod; for (int x = 0; msk; msk >>= 1, x <<= 1) { fill_n(u.begin(), m * 2, 0); x |= !!(n & msk); if (x < m) u[x] = 1 % mod; else {// can be optimized by fft/ntt for (int i = 0; i < m; ++i) { for (int j = 0, t = i + (x & 1); j < m; ++j, ++t) { u[t] = (u[t] + v[i] * v[j]) % mod; } } for (int i = m * 2 - 1; i >= m; --i) { for (int j = 0, t = i - m; j < m; ++j, ++t) { u[t] = (u[t] + trans[j] * u[i]) % mod; } } } v = {u.begin(), u.begin() + m}; } ll ret = 0; for (int i = 0; i < m; ++i) { ret = (ret + v[i] * init[i]) % mod; } return ret; } vec init, trans; ll mod; int m; }; 8. 数值积分
8.1 复化 Simpson 公式
#define rep(i,a,b) for(ll i=a;i<=b;i++) //复化Simpson公式 O(n) template<typename type,typename fun> type Simpson(fun f,type a,type b,int n){ type h=(b-a)/n,s=0; rep(i,0,n-1)s+=h/6.0*(f(a+h*i)+4.0*f(a+h*(i+0.5))+f(a+h*(i+1.0))); return s; } 8.2 Romberg 算法
#define rep(i,a,b) for(ll i=a;i<=b;i++) //Romberg算法 O(2^n) template<typename type,typename fun> type Romberg(fun f,type a,type b,int n){ assert(n>=3); type h=b-a,T[n+2],S[n+2],C[n+2],R[n+2]; T[1]=h/2.0*(f(a)+f(b)); rep(t,0,n-1){ type s=0; rep(i,0,(1<<t)-1)s+=f(a+h*(i+0.5)); T[t+2]=T[t+1]/2.0+h*s/2.0; h/=2.0; } rep(i,1,n)S[i]=T[i+1]*4.0/3.0-T[i]/3.0; rep(i,1,n-1)C[i]=S[i+1]*16.0/15.0-S[i]/15.0; rep(i,1,n-2)R[i]=C[i+1]*64.0/63.0-C[i]/63.0; return R[n-2]; } 类似的话题
-
算法竞赛中的数论训练，就像在迷宫里寻找宝藏，每一步都需要扎实的理论基础和巧妙的计算技巧。要想在这片领域游刃有余，就得循序渐进，把基础打牢，然后逐步攻克难题。第一步：夯实基础，筑牢地基 核心概念的理解： 整除与模运算： 这是数论的基石，要深刻理解 $a pmod{m}$ 的含义，以及整............. -
算法竞赛中的焦虑，这几乎是所有参赛者都绕不过去的坎。说实话，我刚开始参加比赛的时候，也像个无头苍蝇一样，脑子里全是“万一做不出来怎么办？”“别人是不是已经秒出了？”“这次又要垫底了吧？”…… 这种感觉，就像是手里攥着一把沙子，越是想握紧，它就越是从指缝溜走，只留下空虚和恐慌。后来慢慢摸索，也吃了些苦.............
-
DeepMind 在 2021 年 12 月公布的 AlphaCode，在算法竞赛领域引起了巨大的轰动，也成为了人工智能领域的一项重要里程碑。评价 AlphaCode 需要从多个维度进行深入分析，包括其技术实现、性能表现、潜在影响以及局限性。 AlphaCode 评价：一次深入的剖析 1. 技术实现.............
-
2022 年秋招算法岗的竞争有多激烈？这事儿，相信但凡投过简历、面过试的同行们，心里都有数。简直就是一场“千军万马过独木桥”的硬仗，而且这桥还建在悬崖边上，下面是无尽的深渊——找不到工作的焦虑。咱们从几个方面掰开了聊聊，为啥今年算法岗的战场会这么硝烟弥漫：一、 “算法”这个词自带的光环与误读： .............
-
“我国电竞人才缺口达50万”——这个数字一出来，确实让人挺惊讶的，也引发了不少讨论。那么，咱们就好好聊聊这个事儿。首先，怎么看待这个“50万”的缺口？这个数字，背后反映的是电竞产业蓬勃发展但人才供给跟不上的现实。你可以理解成，现在国内做电竞这行的人才，离行业实际发展的需求，还差了差不多50万个“对口.............
-
将并行计算纳入算法竞赛，这是一个非常有趣且具有深远意义的问题。答案是总体上合适的，甚至可以说是必然的趋势。然而，是否合适以及如何合适地纳入，需要我们详细地探讨其原因、挑战、潜在的好处以及实现方式。 一、为什么将并行计算纳入算法竞赛是合适的？ 1. 符合现实世界的计算需求和趋势： 硬件发展驱动： .............
-
“Tourist” 是 Ilya Zakharyaschev 的游戏名，他被广泛认为是历史上最强的算法竞赛选手之一。然而，要断言他是“史上最强”是一个复杂的问题，因为这涉及到对“强”的定义、不同时代比赛的难度以及选手个人特点的评估。以下是支持 Tourist 可能是史上最强的一些论据以及需要考虑的方.............
-
信息学竞赛，尤其是像IOI（国际信息学奥林匹克竞赛）这种级别的，说它是“边缘竞赛”，那绝对是有点冤枉了。但如果要说它“大众化”到像高考语文、数学那样人人皆知，那确实也还差点火候。所以，我们可以说信息学竞赛，在某些圈子里是顶流，但在更广阔的社会认知里，还处在一个“逐渐被看见”的阶段。咱们得把“边缘”这.............
-
今年校招的情况确实非常卷，尤其是机器学习算法岗，应届生数量简直可以用“洪流”来形容。身边不少朋友都深有体会，投出去简历，要么石沉大海，要么收到的面试通知也是寥寥无几，即便拿到offer，薪资待遇也可能不如预期。这让很多人开始担忧，是不是现在涌入机器学习这个赛道的人太多了，再过几年，这个领域会不会饱和.............
-
国际石油市场，就其本质而言，更接近于 寡头市场，并且带有 垄断竞争 和 一定程度的卡特尔（Cartel） 特征。简单地说，它不是一个整齐划一的、符合任何一种经典市场模型的产品。让我们一层一层地剥开来看，为什么会是这样。为什么不是完全竞争？完全竞争市场有几个铁打不动的特征： 无数的买家和卖家： 没.............
-
美团回应取消支付宝支付，并甩出“饿了么 App 也不支持微信支付”这样一句回击，这事儿在很多人看来，确实够味儿，也确实引发了不少讨论。要说这是不是正常的行业竞争，得掰开了揉碎了好好说道说道。首先，我们得明白这事的背景。美团和饿了么，这是国内即时配送领域绝对的两大巨头，可以说是“二龙戏珠”，用户粘性非.............
-
这是一个非常深刻且重要的问题，触及了现代社会最核心的议题之一：数据主权。算法确实离不开大数据，而大数据的基础正是我们每一个人在数字世界留下的痕迹。因此，“我们是不是应该拥有主导数据的权利？” 这个问题的答案，在道德、法律和技术层面都值得深入探讨。一、 问题的根源：我们是谁，以及数据如何与我们关联首先.............
-
选择一本优秀的算法书是踏上算法学习之旅的关键一步。一本好的算法书不仅能帮助你理解枯燥的概念，更能激发你的学习兴趣，培养解决问题的能力。下面我将从多个维度为你详细讲解如何选择算法书：一、明确你的学习目标和基础：这是选择算法书最重要的一步。在开始挑选之前，问自己几个问题： 你的编程语言偏好？ （C+.............
-
作为算法工程师，尤其是在负责算法策略制定和优化时，我们常常会面临一个绕不开的挑战：不确定性。这就像在浓雾中航行，你知道目标在那里，但前进的路充满了未知。尤其是在策略上线后，如果效果不如预期，甚至出现负面影响，那么随之而来的绩效考核和职业发展问题，更是让我们感到压力山大。那么，面对这种“没效果”的风险.............
-
在我看来，大多数情况下我们确实会优先选择时间复杂度更低的算法，也就是在您提到的例子里是算法A（O(n²)）。毕竟，更快的执行速度通常意味着更好的用户体验和更少的资源消耗。但是，任何事情都有例外，而且在某些特定的场景下，一个看似“慢”的算法（比如算法B，O(n³)）反而可能比一个“快”的算法（算法A，.............
-
算法岗位的招聘门槛，尤其是在大厂和研究机构，确实让不少求职者感到焦虑。很多人觉得“没有顶会就别想进”，这话虽然有点绝对，但背后折射出的确实是当下算法岗位对求职者硬实力，特别是科研能力的极高要求。顶会，为什么这么重要？首先，我们需要理解为什么顶会论文会被如此看重。 技术实力的直接证明： 顶会（如 .............
-
算法工程师的“落地能力”，这可不是一句空洞的“会调参”就能概括的。它是一个工程师真正能把脑袋里的模型、想法，变成实实在在能运行、能产生价值的东西的能力的总和。用大白话说，就是能把算法从纸面、从notebook里，真正变成一部能稳定运转、能解决实际问题的机器。这能力涵盖的维度非常广，绝不是只懂模型架构.............
-
好，我们来深入探讨一下有向图强连通分量的求解算法（Kosaraju算法）以及你提出的那个修改方案的错误之处。Kosaraju算法是求解有向图强连通分量的经典算法之一，它的核心思想是利用两次深度优先搜索（DFS）以及图的转置（反向图）。Kosaraju算法的步骤详解：1. 第一次DFS（在原图上）：.............
-
作为一个算法工程师，你是否应该持续阅读论文？这是一个非常值得深思的问题，而且答案绝对不是简单的“是”或“否”。我认为，持续阅读论文是算法工程师职业生涯中一个极其重要但又需要策略性对待的环节，它关乎你的技术深度、广度以及未来的发展潜力。我们先来抛开AI的痕迹，像一个真实的、有经验的同行一样，掰开了揉碎.............
-
这年头，算法岗的风口好像真的有点变了。曾经风光无限，现在听到的更多是“诸神黄昏”和“内卷”。对于很多想入行或者已经在算法领域摸爬滚打的初级从业者来说，这无疑是个让人焦虑的问题：我的路在哪里？我该怎么选？别急，我们先冷静下来，看看现在这个“诸神黄昏”和“内卷”到底是怎么回事，再聊聊我们该怎么在这个新局.............
本站所有内容均为互联网搜索引擎提供的公开搜索信息，本站不存储任何数据与内容，任何内容与数据均与本站无关，如有需要请联系相关搜索引擎包括但不限于百度，google,bing,sogou 等
© 2025 tinynews.org All Rights Reserved. 百科问答小站 版权所有