数论学习笔记
本文中出现的代数式、数,如无特殊说明,均为整数。
对于 OI 中的数论题,答案不对请先尝试:
#define int __int128
typedef long long ll;
#define ll __int128
数论函数指的是定义域为正整数的函数,也可以视为一个序列。
基础部分
\(\perp\) 在数论中表示“互质”,例如 \(2\perp3\)。
对于一个序列 \(f_n=af_{n-1}+bf_{n-2}\),若 \(\gcd(a,b)=1\),则有:
质数
Fermat 素性测试
由费马小定理,对于质数 \(p\),若 \(\gcd(a,p)=1\),满足 \(a^{p-1}\equiv 1\pmod p\),\(a^{p-2}\equiv a^{-1}\pmod p\)。
但是 \(a^{p-1}\equiv1\pmod p\) 并不能推出 \(p\) 为质数。
因此可以随机选择来测试。
Miller–Rabin 素性测试
Miller-Rabin 素性测试,取质数集合 \(A=\lbrace{2,3,5,7,11,13,17,19,23,29,31,37}\rbrace\),则可以通过随机化确定 long long 范围(\([0,2^{64})\))内的任意整数 \(n\) 是否为质数。
从 \(A\) 中取一底数 \(a\),若:
- \(n=a\),\(n\) 为质数。
- \(n\bmod a=0\land n>a\),\(n\) 不为质数。
在都不成立的情况下,进行 Miller-Rabin 素性测试。
将 \(n-1\) 转化:
其中,\(d,r\in\N^*\),也就是说,\(d\) 是 \(n-1\) 在二进制位上的一个前缀,满足前缀的后缀不为 \(0\)。
二次探测定理
使用平方差公式可证。
Miller-Rabin 素性测试的实现
基于 \(a\) 执行 \(r\) 轮测试,通过二次探测定理判断。
参考代码
constexpr const ll A[12]={2,3,5,7,11,13,17,19,23,29,31,37};
bool MillerRabin(ll n){
if(n<=1){
return false;
}
for(ll a:A){
if(n==a){
return true;
}
if(n%a==0){
return false;
}
bool no=true;
ll d=n-1,r=0;
while(!(d&1)){
r++;
d>>=1;
}
ll pl=qpow(a,d,n);
if(pl==1){
no=false;
}
for(int i=0;no&&i<r;i++){
if(pl==n-1){
no=false;
}
pl=(__int128)pl*pl%n;
}
if(no){
return false;
}
}
return true;
}
最大公约数
裴蜀定理
对于 \(\forall a,b\in \Z\):
-
\(\exist x,y\in \Z\),使 \(ax+by=\gcd(a,b)\)。
-
\(\forall x,y\in\Z\),\(\gcd(a,b)\mid(ax+by)\)。
逆定理
\(\forall a,b\in \Z\),若 \(d>0\) 且 \(d\) 为 \(a,b\) 的公因数,且存在 \(x,y\),使得 \(ax+by=d\),则:
乘法逆元
若存在 \(ax\equiv 1\pmod p\),则称 \(x\) 为 \(a\) 关于 \(p\) 的逆元,记 \(x=a^{-1}\)。
当 \(\gcd(a,p)=1\) 的时候(\(p\) 为质数)逆元一定存在。
费马小定理
对于质数 \(p\),若 \(\gcd(a,p)=1\),满足 \(a^{p-1}\equiv 1\pmod p\),\(a^{p-2}\equiv a^{-1}\pmod p\)。
常用于分数取模,即:
详见费马小定理。
欧拉定理
欧拉函数
记 $\varphi(n)$ 为 $1\sim n$ 中与 $n$ 互质的数的个数。
令 $n=p_1^{c_1}p_2^{c_2}p_3^{c_3}\cdots p_k^{c_k}$,满足 $p_i$ 为质数,$c_i>0$。
则:
$$ \varphi(n)=n\left(1-\frac 1{p_1}\right)\left(1-\frac 1{p_2}\right)\left(1-\frac 1{p_3}\right)\cdots\left(1-\frac 1{p_k}\right) $$
推导见容斥原理求解欧拉函数。
若 \(\gcd(a,p)=1\),则 \(a^{\varphi(p)}\equiv1\pmod p\)。
显然,当 \(p\) 为质数时,有 \(\varphi(p)=p-1\),即费马小定理。因此,费马小定理为欧拉定理的特殊形式。
扩展欧拉定理
对于任意的 \(a,p,b\in \N^*\),若满足 \(\gcd(a,p)>1,b\geq \varphi(p)\),有:
当 \(\gcd(a,p)=1\) 时即欧拉定理。
可以使用初等证明或群论证明。
常用于降幂。
扩展欧几里得算法 exgcd
注意不是“类欧几里得算法”。
用于求关于 \(x,y\) 的不定方程 \(ax+by=\gcd(a,b)\) 的一组整数特解。
不妨令 \(a>b\)。
考虑到 \(a=\left\lfloor\dfrac{a}{b}\right\rfloor b+a\bmod b\),代入可得:
因此:
令 \(a'=b,x'=\left\lfloor\dfrac ab\right\rfloor x+y,b'=a\bmod b,y'=x\),有:
显然,可以递归求解。
那么直到 \(b=0\) 时,可以直接解得一组整数特解 \(\begin{cases}x=1\\y=0\end{cases}\)。
设最终递归出来的解为 \(\begin{cases}x=x_0\\y=y_0\end{cases}\)。
令 \(\Delta b=\left\vert\dfrac{b}{\gcd(a,b)}\right\vert,\Delta a=\left\vert\dfrac a{\gcd(a,b)}\right\vert\),则对于 \(\forall k\in \N^*\),方程存在一组解为 \(\begin{cases}x=x_0+k\Delta b\\y=y_0-k\Delta a\end{cases}\)。
参考代码
void exgcd(int a,int &x,int b,int &y){
if(!b){
x=1;
y=0;
return;
}
int tmp;
exgcd(b,tmp,a%b,x);
y=tmp-a/b*x;
}
exgcd 求乘法逆元
如果需要讨论 \(a\) 模 \(p\) 意义下的乘法逆元 \(a^{-1}\),显然有 \(\gcd(a,p)=1\)。
那么可列方程 \(ax+py=1\),显然 \(ax\equiv1\pmod p\)。
则通过扩展欧几里得算法求出 \(x\),即求出了 \(a\) 的逆元。
求出 \(x\) 后,为了使 \(x\) 为正整数,只需要让 \(x\leftarrow(x\bmod p+p)\bmod p\) 即可。
参考代码
//#include<bits/stdc++.h>
#include<algorithm>
#include<iostream>
#include<cstring>
#include<iomanip>
#include<cstdio>
#include<string>
#include<vector>
#include<cmath>
#include<ctime>
#include<deque>
#include<queue>
#include<stack>
#include<list>
using namespace std;
void exgcd(int a,int &x,int b,int &y){
if(!b){
x=1;
y=0;
return;
}
int tmp;
exgcd(b,tmp,a%b,x);
y=tmp-a/b*x;
}
int main(){
/*freopen("test.in","r",stdin);
freopen("test.out","w",stdout);*/
ios::sync_with_stdio(false);
cin.tie(0);cout.tie(0);
int a,b;
cin>>a>>b;
int x0,y0;
exgcd(a,x0,b,y0);
x0%=b;
if(x0<0){
x0+=b;
}
cout<<x0<<'\n';
cout.flush();
/*fclose(stdin);
fclose(stdout);*/
return 0;
}
例题:正整数最大/小解
给定不定方程 \(ax+by=c\),\(a,b,c\in \N^*\),判断其是否有解。
若有正整数解,分别求 \(x,y\) 的最大/小值。
否则,求 \(x,y\) 得最小正整数解。
对于无解,即 \(\gcd(a,b)\nmid c\),很好判断。
先通过 exgcd 求出一组特解 \((x_0,y_0)\),则 \(\begin{cases}x=x_0+k\Delta b\\y=y_0-k\Delta a\end{cases}\) 也是原方程的解。\(\Delta b\) 为 \(b\) 的单次变化量,最小为 \(\Delta b=\left\vert\dfrac{b}{\gcd(a,b)}\right\vert\),\(\Delta a=\left\vert\dfrac a{\gcd(a,b)}\right\vert\)。
先求 \(x\) 的最小正整数解 \(x_{\min}\),即 \(x_0+k\Delta b\) 最小,显然 \(x_{\min}=x_0\bmod \Delta b\)。
但是为了保证 \(x_{\min}\) 为正整数,因此当 \(x_0\bmod \Delta b\leq0\) 时,\(x_{\min}=x_0\bmod \Delta b+\Delta b\)。
\(a,b>0\),因此 \(x\) 最小时 \(y\) 最大,因此可以确定 \(y_{\max}=\dfrac{c-ax_{\min}}{b}\)。
同理,可求出 \(x_{\max},y_{\min}\)。
当 \(x_{\min}\) 为正整数,此时若 \(y_{\max}>0\),则存在正整数解 \(\begin{cases}x=x_{\min}\\y=y_{\max}\end{cases}\),因此可以判断正整数解(判断 \(x_{\max}>0\) 也行)。
正整数解的个数即:
因为从 \(x_{\min}\) 开始,\(x_{\min}+k\Delta b\) 对应一组 \(x,y\)。
参考代码
//#include<bits/stdc++.h>
#include<algorithm>
#include<iostream>
#include<cstring>
#include<iomanip>
#include<cstdio>
#include<string>
#include<vector>
#include<cmath>
#include<ctime>
#include<deque>
#include<queue>
#include<stack>
#include<list>
using namespace std;
#define int ll
typedef long long ll;
int gcd(int a,int b){
while(b){
int tmp=a;
a=b;
b=tmp%b;
}
return a;
}
void exgcd(int a,int &x,int b,int &y){
if(!b){
x=1;
y=0;
return;
}
int tmp;
exgcd(b,tmp,a%b,x);
y=tmp-a/b*x;
}
main(){
/*freopen("test.in","r",stdin);
freopen("test.out","w",stdout);*/
ios::sync_with_stdio(false);
cin.tie(0);cout.tie(0);
int T;
cin>>T;
while(T--){
int a,b,c;
cin>>a>>b>>c;
int gcdAB=gcd(a,b);
if(c%gcdAB){
cout<<"-1\n";
}else{
int w=c/gcdAB;
int x0,y0;
exgcd(a,x0,b,y0);
x0*=w,y0*=w;
// cerr<<"x0="<<x0<<" y0="<<y0<<endl;
int xMin,xMax,yMin,yMax;
int deltaA=a/gcdAB;
int deltaB=b/gcdAB;
xMin=x0%deltaB;
if(xMin<=0){
xMin+=deltaB;
}
yMax=(c-a*xMin)/b;
yMin=y0%deltaA;
if(yMin<=0){
yMin+=deltaA;
}
xMax=(c-b*yMin)/a;
// cerr<<"A="<<a<<" B="<<b<<endl;
// cerr<<"ΔA="<<deltaA<<" ΔB="<<deltaB<<endl;
// cerr<<"xMax="<<xMax<<" xMin="<<xMin<<" yMax="<<yMax<<" yMin="<<yMin<<endl;
if(yMax>0){
cout<<(xMax-xMin)/deltaB+1<<' '<<xMin<<' '<<yMin<<' '<<xMax<<' '<<yMax<<'\n';
}else{
cout<<xMin<<' '<<yMin<<'\n';
}
}
}
cout.flush();
/*fclose(stdin);
fclose(stdout);*/
return 0;
}
/*
1
3 18 6
2 1
*/
exgcd 与同余方程
exgcd 求解的是不定方程:
而同余方程形如:
可以将同余方程转化为:
进而使用 exgcd 求解。
线性求逆元
有一种方法,可以在 \(\mathcal O(\log V+n)\) 的时间内求出任意 \(n\) 个数 \(a_1,a_2,\cdots,a_n\) 关于 \(p\) 的逆元。
记:
那么就可以 \(\mathcal O\left(\log V\right)\) 计算:
\(n-1\sim 1\) 递推,可以推出:
随后再次递推,可以得到 \(a_i\) 关于 \(p\) 的逆元 \(\textit{inv}_i\):
参考代码
pre[0]=1;
for(int i=1;i<=n;i++){
pre[i]=1ll*pre[i-1]*a[i]%P;
}
inv[n]=qpow(pre[n],P-2);
for(int i=n-1;i>=1;i--){
inv[i]=1ll*inv[i+1]*(a[i+1])%P;
}
for(int i=1;i<=n;i++){
inv[i]=1ll*inv[i]*pre[i-1]%P;
}
CRT 中国剩余定理
CRT 中的运算溢出
在使用 CRT 或 exCRT 时,一律建议使用 __int128,防止溢出。
尽管题目大多会保证 $\operatorname{lcm}(p_1,p_2,\cdots,p_n)\leq10^{18}$,但是中间计算时会溢出。
CRT(中国剩余定理)被用于求解线性同余方程组:
其中,\(p_1,p_2,p_3,\cdots,p_n\) 两两互质。
令 \(L=\prod\limits_{i=1}^np_i\),则对于 \(k\in\N\) 有 \(k\cdot L\equiv0\pmod{p_i}\)。
记 \(q_i=\dfrac{L}{p_i}\),设 \(c_i\equiv q_i^{-1}\pmod{p_i}\),即 \(q_i\) 模 \(p_i\) 意义下的逆元。
则,方程组的最小整数解为:
同时,对于 \(\forall k\in\N^*\),\(x+kL\) 均为原方程组的解。
CRT 其实就是一个构造式的做法,易证 \(\forall i\neq j,a_iq_ic_i\equiv0\pmod{p_j}\)。
例题:中国剩余定理(CRT)/ 曹冲养猪
题目中的 \(a_i,b_i\) 即分别为 \(p_i,a_i\)。
直接套 CRT 即可,但是需要注意的是,计算 \(a_iq_ic_i\) 时需要 __int128,会爆 long long。
参考代码
//#include<bits/stdc++.h>
#include<algorithm>
#include<iostream>
#include<cstring>
#include<iomanip>
#include<cstdio>
#include<string>
#include<vector>
#include<cmath>
#include<ctime>
#include<deque>
#include<queue>
#include<stack>
#include<list>
using namespace std;
typedef long long ll;
constexpr const int N=10;
int n,a[N+1],p[N+1];
void exgcd(ll a,ll &x,ll b,ll &y){
if(!b){
x=1;
y=0;
return;
}
ll tmp;
exgcd(b,tmp,a%b,x);
y=tmp-a/b*x;
}
ll inverse(ll a,ll p){
ll x,y;
exgcd(a,x,p,y);
x%=p;
if(x<0){
x+=p;
}
return x;
}
int main(){
/*freopen("test.in","r",stdin);
freopen("test.out","w",stdout);*/
ios::sync_with_stdio(false);
cin.tie(0);cout.tie(0);
cin>>n;
ll L=1;
for(int i=1;i<=n;i++){
cin>>p[i]>>a[i];
L*=p[i];
}
ll x=0;
for(int i=1;i<=n;i++){
ll q=L/p[i],c=inverse(q,p[i]);
x=(x+(__int128)1*a[i]%L*q%L*c)%L;
}
if(x<0){
x+=L;
}
cout<<x<<'\n';
cout.flush();
/*fclose(stdin);
fclose(stdout);*/
return 0;
}
扩展 CRT/exCRT
exCRT 实际上与 CRT 关系不大,如同 exLucas 与 Lucas 一般。
exCRT 解决的是当模数 \(p_1,p_2,p_3,\cdots,p_n\) 不一定两两互质的时候(这时有可能无解)。
先考虑两个同余方程:
设 \(r,s\),将其转化为两个不定方程:
因此可得 \(r\cdot p_1-s\cdot p_2=a_2-a_1\)。
由裴蜀定理,若 \((a_2-a_1)\nmid\gcd(p_1,-p_2)\),则该不定方程无解,则原同余方程组无解。
否则通过 exgcd 求出 \(r,s\),则:
既然 \(x=a_1+r\cdot p_1\),因此有:
这样取 \(\operatorname{lcm}(p_1,p_2)\) 是为了将两个方程合并。
这样,两两顺次合并,最终合并为一个方程即可解。
参考代码
for(int i=2;i<=n;i++){
ll w=(a[i]-a[i-1])/gcd(p[i-1],-p[i]);
ll r,s;
exgcd(p[i-1],r,-p[i],s);
r*=w,s*=w;
a[i]=(p[i-1]*r+a[i-1])%lcm(p[i-1],p[i]);
p[i]=lcm(p[i-1],p[i]);
if(i==n){
if(a[n]<0){
a[n]+=p[n];
}
cout<<a[n]<<'\n';
}
}
例题参考代码
//#include<bits/stdc++.h>
#include<algorithm>
#include<iostream>
#include<cstring>
#include<iomanip>
#include<cstdio>
#include<string>
#include<vector>
#include<cmath>
#include<ctime>
#include<deque>
#include<queue>
#include<stack>
#include<list>
using namespace std;
namespace IO{
inline char getchar(){
static int p1,p2;
static char buf[1<<20];
if(p1==p2)p2=fread(buf,1,1<<20,stdin),p1=0;
return (p1==p2?EOF:buf[p1++]);
}
template<typename T>
inline void scanf(T &x){
x=0;
register int f=1;
register char ch=IO::getchar();
for(;ch<'0'||'9'<ch;ch=IO::getchar());
for(;'0'<=ch&&ch<='9';ch=IO::getchar())x=(x<<3)+(x<<1)+ch-'0';
x*=f;
}
static int p;
static char pbuf[1<<20];
inline void putchar(char ch){
pbuf[p++]=ch;
if(p==(1<<20)){
fwrite(pbuf,1,1<<20,stdout);
p=0;
}
}
template<typename T>
inline void printf(T x){
static char s[101];
int top=0;
do{
s[++top]=x%10+'0';
x/=10;
}while(x);
while(top)IO::putchar(s[top--]);
}
inline void end(){
fwrite(pbuf,1,p,stdout);
}
struct Tool{
~Tool(){
end();
}
}tool;
}
typedef __int128 lll;
#define int lll
#define ll lll
//typedef long long ll;
constexpr const int N=1e5;
int n;
ll a[N+1],p[N+1];
void exgcd(ll a,ll &x,ll b,ll &y){
if(!b){
x=1;
y=0;
return;
}
ll tmp;
exgcd(b,tmp,a%b,x);
y=tmp-a/b*x;
}
ll gcd(ll a,ll b){
while(b){
ll tmp=a;
a=b;
b=tmp%b;
}
return a;
}
ll lcm(ll a,ll b){
return a/gcd(a,b)*b;
}
main(){
/*freopen("test.in","r",stdin);
freopen("test.out","w",stdout);*/
/*ios::sync_with_stdio(false);
cin.tie(0);cout.tie(0);*/
IO::scanf(n);
for(int i=1;i<=n;i++){
IO::scanf(p[i]);
IO::scanf(a[i]);
}
//合并(i-1,i)
for(int i=2;i<=n;i++){
ll w=(a[i]-a[i-1])/gcd(p[i-1],-p[i]);
ll r,s;
exgcd(p[i-1],r,-p[i],s);
r*=w,s*=w;
a[i]=(p[i-1]*r+a[i-1])%lcm(p[i-1],p[i]);
p[i]=lcm(p[i-1],p[i]);
if(i==n){
if(a[n]<0){
a[n]+=p[n];
}
IO::printf(a[n]);
}
}
cout.flush();
/*fclose(stdin);
fclose(stdout);*/
return 0;
}
exexCRT
即形如:
多了一个系数 \(f_i\)。
同样可以用类似 exCRT 的方法两两合并解决,具体见此。
离散对数
所谓离散对数,即模意义下取对数。
形如:给定模数 \(p\),及整数 \(a,b\),求整数 \(x\) 使得 \(a^x\equiv b\pmod p\)。
注意:离散对数可能不存在。
BSGS 算法
在 OI 中,BSGS 算法(Baby-Step Giant-Step,大步小步算法北上广深算法)常用来求解离散对数。
BSGS 算法要求 \(a\perp p\),求 \(x\) 使得 \(a^x\equiv b\pmod p\)。
若有解,则存在 \(x\leq\varphi(p)\) 的解;因为欧拉定理说明,\(a\perp p\) 时,\(a^{\varphi(p)}\equiv1\pmod p\)。
若枚举 \(x\) 是 \(\mathcal O(\varphi(p))\) 的,当 \(p\) 为质数的时候不能接受,因此可以考虑分块。
令 \(B=\left\lceil\sqrt{\varphi(p)}\right\rceil\),\(x=qB+r\),且 \(0\leq q,r\leq B\)。
则有:
\(a\perp p\),则 \(a\) 的乘法逆元 \(a^{-1}\) 一定存在,有:
枚举 \(r\),将其与其对应的 \(b\left(a^{-1}\right)^r\) 一同存储在数据结构(一般是哈希表)中,随后枚举 \(q\),在数据结构中找 \(a^{qB}\) 对应的 \(r\),找到了 \(r\) 便找到了一个解 \(x=qB+r\)。同时,因为 \(a\perp p\),因此 \(a^{-1},a^{-2},a^{-3},\cdots,a^{-B}\) 模 \(p\) 意义下互不相同。
时间复杂度:\(\mathcal O\left(\sqrt{\varphi(p)}\right)\),在 \(p\) 为质数时最劣,\(\mathcal O\left(\sqrt p\right)\)。
参考代码
B=ceil(sqrt(euler(P)));
c=qpow(a,B);
for(int r=0,powR=1,R=qpow(a,P-2);r<B;r++){
m[1ll*b*powR%P]=r;
powR=1ll*powR*R%P;
}
for(int q=0,powC=1;q<=B;q++){
if(m.count(powC)){
cout<<q*B+m[powC]<<endl;
return 0;
}
powC=1ll*powC*c%P;
}
cout<<"no solution\n";
参考代码
//#include<bits/stdc++.h>
#include<algorithm>
#include<iostream>
#include<cstring>
#include<iomanip>
#include<cstdio>
#include<string>
#include<vector>
#include<cmath>
#include<ctime>
#include<deque>
#include<queue>
#include<stack>
#include<list>
#include<unordered_map>
using namespace std;
int a,b,P,B,c;
unordered_map<int,int>m;
int euler(int n){
int ans=n;
for(int i=2;1ll*i*i<=n;i++){
if(n%i==0){
ans=ans/i*(i-1);
while(n%i==0){
n/=i;
}
}
}
if(n>1){
ans=ans/n*(n-1);
}
return ans;
}
int qpow(int base,int n){
int ans=1;
while(n){
if(n&1){
ans=1ll*ans*base%P;
}
base=1ll*base*base%P;
n>>=1;
}
return ans;
}
int main(){
/*freopen("test.in","r",stdin);
freopen("test.out","w",stdout);*/
ios::sync_with_stdio(false);
cin.tie(0);cout.tie(0);
cin>>P>>a>>b;
B=ceil(sqrt(euler(P)));
c=qpow(a,B);
for(int r=0,powR=1,R=qpow(a,P-2);r<B;r++){
m[1ll*b*powR%P]=r;
powR=1ll*powR*R%P;
}
for(int q=0,powC=1;q<=B;q++){
if(m.count(powC)){
cout<<q*B+m[powC]<<endl;
return 0;
}
powC=1ll*powC*c%P;
}
cout<<"no solution\n";
cout.flush();
/*fclose(stdin);
fclose(stdout);*/
return 0;
}
exBSGS 扩展 BSGS 算法
exBSGS 可解决 \(a\perp p\) 不成立的情况。
由扩展欧拉定理可知,若 \(x\) 有解,则 \(x\leq2\varphi(p)\) 范围内也有解。
首先特判掉 \(x=0\) 的情况,即 \(b=1\) 或 \(p=1\)。
考虑分块,令 \(B=\left\lceil\sqrt{2\varphi(p)}\right\rceil,x=qB-r,0\leq q,r\leq B\),则有:
因此可以预处理 \(a^{qB}\bmod p\),随后枚举 \(r\);但是注意到前者是后者的充分条件而不是充要条件,因此找出来的 \(x=qB+r\) 只是“可能的解”,需要检验。
对于多个不同的 \(q\) 产生的模 \(p\) 意义下相同的 \(a^{qB}\),可以只存储 \(q\) 最小的两个。
证明
由扩展欧拉定理,$a^x$ 是在 $k\leq\varphi(p)$ 值后以 $\varphi(p)$ 为周期循环的。
循环周期内的 $a^{qB}$ 显然至多存储 $1$ 个,而非循环部分 $a^{qB}$ 也至多存储 $1$ 个。
即:保留最小的 $2$ 个。
时间复杂度:仍然为 \(\mathcal O\left(\sqrt{\varphi(p)}\right)\)。
参考代码
a%=P;b%=P;
if(b==1||P==1){
cout<<"0\n";
continue;
}
B=ceil(2*sqrt(euler(P)));
c=qpow(a,B);
for(int q=0,powC=1;q<B;q++){
auto &pl=m[powC];
if(pl.size()<2){
pl.push_back(q);
}
powC=1ll*c*powC%P;
}
int x=2147483647;
for(int r=0,powA=1;r<B;r++){
auto &pl=m[1ll*b*powA%P];
for(int q:pl){
int x0=(1ll*q*B-r)%P;
if(x0<0){
continue;
}
if(qpow(a,x0)==b){
x=min(x,x0);
}
}
powA=1ll*powA*a%P;
}
if(x<2147483647){
cout<<x<<'\n';
}else{
cout<<"No Solution\n";
}
参考代码
//#include<bits/stdc++.h>
#include<algorithm>
#include<iostream>
#include<cstring>
#include<iomanip>
#include<cstdio>
#include<string>
#include<vector>
#include<cmath>
#include<ctime>
#include<deque>
#include<queue>
#include<stack>
#include<list>
#include<unordered_map>
#include<ext/pb_ds/assoc_container.hpp>
#include<ext/pb_ds/tree_policy.hpp>
using namespace std;
int a,b,P,B,c;
//pb_ds 哈希表常数较 unordered_map 更小
__gnu_pbds::gp_hash_table<int,vector<int> >m;
//unordered_map<int,vector<int> >m;
//map<int,vector<int> >m;
int euler(int n){
int ans=n;
for(int i=2;1ll*i*i<=n;i++){
if(n%i==0){
ans=ans/i*(i-1);
while(n%i==0){
n/=i;
}
}
}
if(n>1){
ans=ans/n*(n-1);
}
return ans;
}
int qpow(int base,int n){
int ans=1;
while(n){
if(n&1){
ans=1ll*ans*base%P;
}
base=1ll*base*base%P;
n>>=1;
}
return ans;
}
main(){
/*freopen("test.in","r",stdin);
freopen("test.out","w",stdout);*/
ios::sync_with_stdio(false);
cin.tie(0);cout.tie(0);
while(true){
m.clear();
cin>>a>>P>>b;
if(a==P&&P==b&&b==0){
break;
}
a%=P;b%=P;
if(b==1||P==1){
cout<<"0\n";
continue;
}
B=ceil(2*sqrt(euler(P)));
c=qpow(a,B);
for(int q=0,powC=1;q<B;q++){
auto &pl=m[powC];
if(pl.size()<2){
pl.push_back(q);
}
powC=1ll*c*powC%P;
}
int x=2147483647;
for(int r=0,powA=1;r<B;r++){
auto &pl=m[1ll*b*powA%P];
for(int q:pl){
int x0=(1ll*q*B-r)%P;
if(x0<0){
continue;
}
if(qpow(a,x0)==b){
x=min(x,x0);
}
}
powA=1ll*powA*a%P;
}
if(x<2147483647){
cout<<x<<'\n';
}else{
cout<<"No Solution\n";
}
}
cout.flush();
/*fclose(stdin);
fclose(stdout);*/
return 0;
}
Lucas 定理
Lucas 定理
Lucas 定理用于求解大组合数取质数模。(对于模数不为质数的情况,请参考exLucas)。
Lucas 定理的内容很简单:
考虑如何证明。
证明
由二项式定理:
考虑 \(\dbinom pi\) 在模 \(p\) 意义下的取值。
因为:
那么如果化简之后,分子 \(p!\) 中的 \(p\) 项没有被约分掉,则有 \(\dbinom pi\equiv p\cdot\dfrac{(p-1)!}{i!(p-i)!}\equiv0\pmod p\)。
因为 \(p\) 为质数,所以 \(p\) 项能被约分掉当且仅当 \(i!\) 中含有 \(p\) 项或 \((p-i)!\) 中含有 \(p\) 项。
即:\(\dbinom pi\not\equiv0\pmod p\) 当且仅当 \(i\equiv0\pmod p\) 或 \(p-i\equiv 0\pmod p\)。
在二项式定理中,\(i\) 满足 \(0\leq i\leq p\),所以 \(\dbinom pi\not\equiv0\) 当且仅当 \(i=0\) 或 \(i=p\)。
在这两种情况中,都可以计算得到 \(\dbinom pi=1\),即 \(\dbinom pi\equiv[i=0\lor i=p]\)。
重新带回二项式定理,可得:
考虑一个二项式 \((1+x)^n\bmod p\)。
由二项式定理,\((1+x)^n\) 的 \(x^m\) 项系数为 \(\dbinom nm\)。
而想要从 \((1+x^p)^{\lfloor\frac np\rfloor}(1+x)^{n\bmod p}\) 中得到 \(x^m\),即从 \((1+x^p)^{\lfloor\frac np\rfloor}\) 中选取 \(\lfloor\frac mp\rfloor\) 个 \(x^p\),再从 \((1+x)^{n\bmod p}\) 中选取 \(m\bmod p\) 个 \(x\)。
即:
你可能有一个问题:这看起来明明应当是一个等式,但为什么是同余呢?
即,Lucas 定理应当表述为:
但是,你要知道,只有在模 \(p\) 意义下才有 \((1+x)^{p\lfloor\frac np\rfloor}=(1+x^p)^{\lfloor\frac np\rfloor}\)。
因此,只有在模 \(p\) 意义下才有:
即:
应用
当 \(n\) 比较大而无法使用其他方法(例如预处理 \(1\sim n\) 的阶乘再利用乘法逆元)直接求解组合数时,可以使用 Lucas 定理。
Lucas 定理只需要递归使用即可,即递归计算 \(\dbinom{\lfloor\frac np\rfloor}{\lfloor\frac mp\rfloor}\),递归终点即 \(\dbinom{0}{0}\)。
其时间复杂度为:\(\mathcal O\left(f(p)\log m\right)\)。
其中,\(f(p)\) 表示单次计算 \(\dbinom{n\bmod p}{m\bmod p}\) 的复杂度,因写法而异。
可以使用乘法逆元,则时间复杂度为 \(\mathcal O(\log p\log m)\)。
也可以 \(\mathcal O(p\log p)\) 递推,时间复杂度为 \(\mathcal O(p\log p\log m)\)。
推荐使用乘法逆元。
很简单,注意是 \(\dbinom{n+m}{n}\) 即可。
参考代码
//#include<bits/stdc++.h>
#include<algorithm>
#include<iostream>
#include<cstring>
#include<iomanip>
#include<cstdio>
#include<string>
#include<vector>
#include<cmath>
#include<ctime>
#include<deque>
#include<queue>
#include<stack>
#include<list>
using namespace std;
const int N=1e5;
int ksm(int a,int n,int p){
if(n==0)return 1;
int t=ksm(a,n>>1,p);
t=1ll*t*t%p;
if(n&1)t=1ll*t*a%p;
return t;
}
int C(int n,int m,int p){
static int ans[N+1]={1};
for(int i=1;i<=m;i++){
ans[i]=1ll*ans[i-1]*(n-i+1)%p*ksm(i,p-2,p)%p;
}return ans[m];
}
int Lucas(int n,int m,int p){
if(m==0)return 1;
return (1ll*Lucas(n/p,m/p,p)*C(n%p,m%p,p))%p;
}
int main(){
/*freopen("test.in","r",stdin);
freopen("test.out","w",stdout);*/
int T,n,m,p;
scanf("%d",&T);
while(T--){
scanf("%d %d %d",&n,&m,&p);
printf("%d\n",Lucas(n+m,n,p));
}
/*fclose(stdin);
fclose(stdout);*/
return 0;
}
exLucas 算法(exLucas 定理)
exLucas 算法可以求 \(\dbinom{n}{m}\bmod P\),\(P\) 不一定为质数。(但是,exLucas 实际上和 Lucas 没有多大关系……)
由唯一分解定理,可以对 \(P\) 进行质因数分解:
CRT 求解
我们如果能够求出 \(a_1,a_2,\cdots,a_n\) 使得:
那么我们就可以通过 CRT 求解出 \(\dbinom nm\bmod P\)。因为在 CRT 中,恰好也有 \(P=p_1^{c_1}p_2^{c_2}\cdots p_n^{c_n}\)。
求解模质数幂下余数
展开定义式:
因此:
因为 \(p_i^{c_i},m!,(n-m)!\) 不一定互质,因此乘法逆元不一定存在。
考虑先提出分子和分母中所有的 \(p_i\) 次幂,随后便可以使用逆元求解。
记 \(x\) 的质因数分解中质数 \(p\) 的幂次为 \(\nu_p(x)\),剩余积为 \((x)_p\),即:
则有:
那么,现在只需要考虑对于 \(x,p\),如何高效地求出 \(\nu_{p}(x!)\bmod p_i^{c_i},(x!)_p\bmod p_i^{c_i}\)。
考虑:
容易发现,在 \(p,2p,3p,\cdots,\left\lfloor\dfrac xp \right\rfloor p\) 中 \(p\) 的个数为 \(\left\lfloor\dfrac xp\right\rfloor+\nu_p\left(\left\lfloor\dfrac xp \right\rfloor!\right)\)。其中 \(\nu_p\left(\left\lfloor\dfrac xp \right\rfloor!\right)\) 是 \(1,2,3,\cdots,\left\lfloor\dfrac xp\right\rfloor\) 中的 \(p\) 的个数。
\(p\) 为质数,所以在 \(1,2,3,\cdots,p-1,p+1,\cdots\) 中不会含有 \(p\)。
将递推式展开:
因此可以 \(\mathcal O(\log_p x)\) 计算 \(\nu_p(x!)\):
int v(ll n,ll p){
int ans=0;
do{
n/=p;
ans+=n;
}while(n);
return ans;//这里取不取模其实都无所谓,因为幂次不会太多,想取也可以.
}
现在考虑如何计算 \((x!)_p\bmod p^c\);显然不能利用定义式 \((x!)_p=\dfrac{x!}{\nu_p(x!)}\),而需要其他方法(因为无法得知 \(x!\))。
不难进行递推:
由 Wilson 定理的推论(\(m=2,4,p^c,2p^c\) 时 \(k\equiv-1\),否则为 \(1\)):
因此:
可以发现,每次计算 \((x!)_p\) 时,\(p^c\) 是固定的,因此可以预处理:
因此,得到最终推导式:
可以递归或迭代处理,时间复杂度 \(\mathcal O(p^c+\log n x)\)。
int Wilson(int n,int p,int pc){
int ans=1;
vector<int>f(pc)
f[0]=1;
for(int i=1;i<pc;i++){
f[i]=i%p?f[i-1]*i%pc:f[i-1];
}
bool flag=p!=2||pc<=4;
while(n>1){
if((n/pc)&flag){
ans=pc-ans;
}
ans=ans*f[n%pc]%pc;
n/=p;
}
return ans;
}
参考代码
//#include<bits/stdc++.h>
#include<algorithm>
#include<iostream>
#include<cstring>
#include<iomanip>
#include<cstdio>
#include<string>
#include<vector>
#include<cmath>
#include<ctime>
#include<deque>
#include<queue>
#include<stack>
#include<list>
using namespace std;
typedef long long ll;
#define ll __int128
#define int __int128
constexpr const int PMAX=1e6;
void exgcd(ll a,ll &x,ll b,ll &y){
if(!b){
x=1;
y=0;
return;
}
ll tmp;
exgcd(b,tmp,a%b,x);
y=tmp-a/b*x;
}
ll inverse(ll a,ll p){
ll x,y;
exgcd(a,x,p,y);
x%=p;
if(x<0){
x+=p;
}
return x;
}
int qpow(int base,int n){
int ans=1;
while(n){
if(n&1){
ans=1ll*ans*base;
}
base=1ll*base*base;
n>>=1;
}
return ans;
}
int CRT(vector<int>a,vector<int>p){
ll L=1;
for(int &i:p){
L*=i;
}
ll x=0;
for(int i=0;i<a.size();i++){
ll q=L/p[i];
x=(x+1ll*a[i]*q%L*inverse(q,p[i])%L)%L;
}
if(x<0){
x+=L;
}
return x;
}
int v(ll n,ll p){
int ans=0;
do{
n/=p;
ans+=n;
}while(n);
return ans;
}
int Wilson(int n,int p,int pc){
int ans=1;
vector<int>f(pc);
f[0]=1;
for(int i=1;i<pc;i++){
f[i]=i%p?f[i-1]*i%pc:f[i-1];
}
bool flag=p!=2||pc<=4;
while(n>1){
if((n/pc)&flag){
ans=pc-ans;
}
ans=ans*f[n%pc]%pc;
n/=p;
}
return ans;
}
void breakDown(int P,vector<int>&p,vector<int>&pc){
int pp=P;
for(int i=2;i*i<=pp;i++){
if(pp%i==0){
p.push_back(i);
pc.push_back(1);
while(pp%i==0){
pc.back()*=i;
pp/=i;
}
}
}
if(pp>1){
p.push_back(pp);
pc.push_back(pp);
}
}
long long exLucas(ll n,ll m,ll P){
if(n<m){
return 0;
}
vector<int>p,pc;
breakDown(P,p,pc);
vector<int>a(pc.size());
for(int i=0;i<p.size();i++){
int nWilson=Wilson(n,p[i],pc[i]);
int mWilson=Wilson(m,p[i],pc[i]);
int nmWilson=Wilson(n-m,p[i],pc[i]);
a[i]=nWilson*inverse(mWilson*nmWilson%pc[i],pc[i])*qpow(p[i],v(n,p[i])-v(m,p[i])-v(n-m,p[i]));
}
return CRT(a,pc);
}
main(){
/*freopen("test.in","r",stdin);
freopen("test.out","w",stdout);*/
ios::sync_with_stdio(false);
cin.tie(0);cout.tie(0);
long long n,m,P;
cin>>n>>m>>P;
cout<<exLucas(n,m,P)<<'\n';
cout.flush();
/*fclose(stdin);
fclose(stdout);*/
return 0;
}
/*
1000000000000000000 500000000000000000 998243
0
*/
参考代码
//#include<bits/stdc++.h>
#include<algorithm>
#include<iostream>
#include<cstring>
#include<iomanip>
#include<cstdio>
#include<string>
#include<vector>
#include<cmath>
#include<ctime>
#include<deque>
#include<queue>
#include<stack>
#include<list>
using namespace std;
typedef long long ll;
#define ll __int128
#define int __int128
constexpr const int M=5;
long long w[M+1];
constexpr const int PMAX=1e5;
void exgcd(ll a,ll &x,ll b,ll &y){
if(!b){
x=1;
y=0;
return;
}
ll tmp;
exgcd(b,tmp,a%b,x);
y=tmp-a/b*x;
}
ll inverse(ll a,ll p){
ll x,y;
exgcd(a,x,p,y);
x%=p;
if(x<0){
x+=p;
}
return x;
}
int qpow(int base,int n){
int ans=1;
while(n){
if(n&1){
ans=1ll*ans*base;
}
base=1ll*base*base;
n>>=1;
}
return ans;
}
int CRT(vector<int>a,vector<int>p){
ll L=1;
for(int &i:p){
L*=i;
}
ll x=0;
for(int i=0;i<a.size();i++){
ll q=L/p[i];
x=(x+1ll*a[i]*q%L*inverse(q,p[i])%L)%L;
}
if(x<0){
x+=L;
}
return x;
}
int v(ll n,ll p){
int ans=0;
do{
n/=p;
ans+=n;
}while(n);
return ans;
}
int Wilson(int n,int p,int pc){
int ans=1;
vector<int>f(pc);
f[0]=1;
for(int i=1;i<pc;i++){
f[i]=i%p?f[i-1]*i%pc:f[i-1];
}
bool flag=p!=2||pc<=4;
while(n>1){
if((n/pc)&flag){
ans=pc-ans;
}
ans=ans*f[n%pc]%pc;
n/=p;
}
return ans;
}
void breakDown(int P,vector<int>&p,vector<int>&pc){
int pp=P;
for(int i=2;i*i<=pp;i++){
if(pp%i==0){
p.push_back(i);
pc.push_back(1);
while(pp%i==0){
pc.back()*=i;
pp/=i;
}
}
}
if(pp>1){
p.push_back(pp);
pc.push_back(pp);
}
}
long long exLucas(ll n,ll m,ll P){
if(n<m){
return 0;
}
vector<int>p,pc;
breakDown(P,p,pc);
vector<int>a(pc.size());
for(int i=0;i<p.size();i++){
int nWilson=Wilson(n,p[i],pc[i]);
int mWilson=Wilson(m,p[i],pc[i]);
int nmWilson=Wilson(n-m,p[i],pc[i]);
a[i]=nWilson*inverse(mWilson*nmWilson%pc[i],pc[i])*qpow(p[i],v(n,p[i])-v(m,p[i])-v(n-m,p[i]));
}
return CRT(a,pc);
}
main(){
/*freopen("test.in","r",stdin);
freopen("test.out","w",stdout);*/
ios::sync_with_stdio(false);
cin.tie(0);cout.tie(0);
long long n,m,P;
cin>>P>>n>>m;
long long ans=1;
for(int i=1;i<=m;i++){
cin>>w[i];
ans=ans*exLucas(n,w[i],P)%P;
n-=w[i];
if(n<0){
cout<<"Impossible"<<endl;
return 0;
}
}
cout<<ans<<'\n';
cout.flush();
/*fclose(stdin);
fclose(stdout);*/
return 0;
}
数论分块
数论分块(整除分块)用于快速计算:
数论分块的原理即本质不同的 \(\left\lfloor\dfrac ni\right\rfloor\) 只有至多 \(\mathcal O\left(\sqrt n\right)\) 种。
证明
- 当 $i\leq\sqrt n$ 时,$i$ 只有 $\mathcal O\left(\sqrt n\right)$ 种取值。
- 当 $i>\sqrt n$ 时,$\dfrac ni\leq\sqrt n$,$\left\lfloor\dfrac ni\right\rfloor$ 只有 $\mathcal O\left(\sqrt n\right)$ 种取值。
故,本质不同的 $\left\lfloor\dfrac ni\right\rfloor$ 只有至多 $\mathcal O\left(\sqrt n\right)$ 种取值。
因此 \(\left\lfloor\dfrac ni\right\rfloor\) 相同值肯定是连续的,那么就可以找出这 \(\mathcal O\left(\sqrt n\right)\) 区间 \([l_1,r_1],[l_2,r_2],\cdots,[l_k,r_k]\)(\(r_i+1=l_{i+1}\)),求出:
将 \(i\) 个答案相加即可。
若分块后可以 \(\mathcal O(1)\) 求出,则数论分块是 \(\mathcal O\left(\sqrt n\right)\) 的。
寻找区间
对于一个 \(l\) 所对应的 \(\left\lfloor\dfrac nl\right\rfloor\),要寻找一个最大的 \(r\) 使得 \(\left\lfloor\dfrac nl\right\rfloor=\left\lfloor\dfrac nr\right\rfloor\),从而找到区间 \([l,r]\)。
则有:
而对于向上取整的情形,即对于 \(l\) 对应的 \(\left\lceil\dfrac nl\right\rceil\),要寻找一个最大的 \(r\) 使得 \(\left\lceil\dfrac nr\right\rceil\),从而找到区间 \([l,r]\)。
则有:
单一参数
给定 \(n,k\in\N^*\),满足 \(1\leq n,k\leq10^9\),求:
\[G(n,k)=\sum_{i=1}^nk\bmod i \]
显然,\(k\bmod i=k-i\left\lfloor\dfrac ki\right\rfloor\)。
因此有:
考虑快速求出 \(\sum\limits_{i=1}^ni\left\lfloor\dfrac ki\right\rfloor\),可以进行数论分块。
对于数论分块中一组确定的 \(l,r\),有:
这可以 \(\mathcal O(1)\) 计算,因此总计算时间复杂度为 \(\mathcal O\left(\sqrt n\right)\)。
参考代码
//#include<bits/stdc++.h>
#include<algorithm>
#include<iostream>
#include<cstring>
#include<iomanip>
#include<cstdio>
#include<string>
#include<vector>
#include<cmath>
#include<ctime>
#include<deque>
#include<queue>
#include<stack>
#include<list>
using namespace std;
#define int ll
typedef long long ll;
ll G(int n,int k){
ll ans=1ll*n*k;
for(int l=1,r;l<=n;l=r+1){
int t=k/l;
if(!t){
r=n;
}else{
r=min(k/t,n);
}
ans-=1ll*(r-l+1)*t*(l+r)>>1;
}
return ans;
}
main(){
/*freopen("test.in","r",stdin);
freopen("test.out","w",stdout);*/
ios::sync_with_stdio(false);
cin.tie(0);cout.tie(0);
int n,k;
cin>>n>>k;
cout<<G(n,k)<<'\n';
cout.flush();
/*fclose(stdin);
fclose(stdout);*/
return 0;
}
多参数
此时寻找 \(r\),会取 \(\min\),即:
给定 \(1\leq n,m\leq10^9\),求:
\[\left(\sum_{i=1}^n\sum_{j=1}^m[i\neq j](n\bmod i)(m\bmod j)\right)\bmod 19940417 \]
不妨令 \(n\leq m\),否则交换 \(n,m\)。
模运算不好处理,转换一下:
\([i\neq j]\) 不好处理,因此可以先全部求出来,再减去相等的部分:
不难发现 \(j\) 与 \(i\) 的取值无关,因此转换为:
于是可以 \(\mathcal O\left(\max(n,m)\right)\) 计算,但是考虑到 \(1\leq n,m\leq10^9\),只需要使用数论分块优化即可。
参考代码
//#include<bits/stdc++.h>
#include<algorithm>
#include<iostream>
#include<cstring>
#include<iomanip>
#include<cstdio>
#include<string>
#include<vector>
#include<cmath>
#include<ctime>
#include<deque>
#include<queue>
#include<stack>
#include<list>
using namespace std;
constexpr const int P=19940417,inv2=9970209,inv6=3323403;
#define ll __int128
//typedef long long ll;
ll g(ll n,ll m){
ll ans=0;
for(ll l=1,r=n;l<=n;l=r+1,r=n){
ll t=m/l;
if(t){
r=min(n,m/t);
}
ans=(ans+t*(l+r)%P*(r-l+1)%P*inv2%P)%P;
}
return ans;
}
ll g2(ll n,ll m){
ll ans=0;
for(ll l=1,r=n;l<=n;l=r+1,r=n){
ll tn=n/l,tm=m/l;
r=min(n/tn,m/tm);
ans+=m*(l+r)%P*(r-l+1)*inv2*tn%P;
ans+=n*(l+r)%P*(r-l+1)*inv2*tm%P;
ans-=(r*(r+1)*(2*r+1)%P-(l-1)*l*(2*l-1)%P)*inv6*tn*tm%P;
ans%=P;
}
return ans;
}
int f(ll n,ll m){
ll ans=(n*n%P-g(n,n))*(m*m%P-g(m,m))%P;
ans=(ans-n*n%P*m%P)%P;
ans=(ans+g2(n,m))%P;
if(ans<0){
ans+=P;
}
return ans;
}
int main(){
/*freopen("test.in","r",stdin);
freopen("test.out","w",stdout);*/
ios::sync_with_stdio(false);
cin.tie(0);cout.tie(0);
int n,m;
cin>>n>>m;
if(n>m){
swap(n,m);
}
cout<<f(n,m)<<'\n';
cout.flush();
/*fclose(stdin);
fclose(stdout);*/
return 0;
}
积性函数
定义
若数论函数 \(f(n)\) 满足 \(f(1)=1\),且对于 \(\forall x,y\in\N^*\land x\perp y\),有 \(f(xy)=f(x)\cdot f(y)\),称 \(f(n)\) 为积性函数。
若数论函数 \(f(n)\) 满足 \(f(1)=1\),且对于 \(\forall x,y\in\N^*\),有 \(f(xy)=f(x)\cdot f(y)\),称 \(f(n)\) 为完全积性函数。
基础部分
若 \(f(n),g(n)\) 均为积性函数,则 \(f(x^p),f^p(x),f(x)g(x),(f*g)(x)\) 均为积性函数。其中 \(f*g\) 为迪利克雷卷积:
常见积性函数
-
单位函数:\(\varepsilon(n)=[n=1]\),完全积性函数。(
\varepsilon) -
恒等函数:\(\mathrm{id}_{k}(n)=n^k\),\(\mathrm{id}_1(n)\) 通常记为 \(\mathrm{id}(n)\),完全积性函数。
-
常数函数:\(1(n)=1\),完全积性函数。
-
除数函数:\(\sigma_k(n)=\sum\limits_{d\mid n}d^k\),\(\sigma_0(n)\) 通常记作 \(d(n)\) 或 \(\tau(n)\),\(\sigma_1(n)\) 通常记作 \(\sigma(n)\)。(
\sigma)$d(n)$ 的求法
令 $n=p_1^{c_1}p_2^{c_2}\cdots p_k^{c_k}$,其中 $p_1,p_2,\cdots,p_k$ 均为质数。
对于每一个指数 $c_i$ 的更改,都会产生新的因数。每一个 $i$,都有 $c_i+1$ 种选择。显然有:
$$ d(n)=\prod_{i=1}^k(c_i+1) $$时间复杂度 $\mathcal O\left(\sqrt n\right)$。
-
欧拉函数:\(\varphi(n)\)。(
\varphi) -
莫比乌斯函数:
\[\mu(n)= \begin{cases} 1&n=1\\ 0&\exist d>1,d^2\mid n\\ (-1)^k&k\text{ 为 }n\text{ 的本质不同质因子个数} \end{cases} \]\(\mu(n)=0\) 即 \(n\) 含有平方因子。(
\mu)
同时,除数函数 \(d(n)\) 即 \(n\) 的因数个数,且 \(d(ij)\) 满足:
证明
令 $i=p_1^{c_1}p_2^{c_2}\cdots p_k^{c_k},j={p'_1}^{c'_1}{p'_2}^{c'_2}\cdots {p'_{k'}}^{c'_{k'}}$。如果将指数对位相加,则可以得到 $ij$ 的质因数分解形式。
将每一个质因子 $p_l$ 枚举 $c_l+c'_l+1$ 次,便可以得到全部的因数。因此枚举 $x\mid i,y\mid j$,但是需要保证 $x,y$ 不包含重复质因子,即 $x\perp y$。
迪利克雷卷积
对于数论函数 \(f,g\),记 \(h=f*g\) 为 \(f,g\) 通过迪利克雷卷积相乘的到的结果,即:
常见迪利克雷卷积
- \(\varphi*1=\mathrm{id}\)。
- \(\mu*1=\varepsilon\)。
- \(\mu*1=\varphi\)。
莫比乌斯函数
莫比乌斯函数 \(\mu(n)\) 满足:
即 \(\mu*1=\varepsilon\)。
因此,\(\mu(n)\) 可用于处理互质相关信息:
莫比乌斯反演/莫比乌斯变换
若 \(g(n)=\sum\limits_{d\mid n}f(d)\),则有:
因为:
若 \(g(n)=\sum\limits_{n|d}f(d)\),则有:
详见此处。
莫比乌斯反演的本质其实就是高维差分。
线性筛
线性筛可以用于求解欧拉函数及莫比乌斯函数。实际上,线性筛几乎可以求解所有的积性函数,求解方法视具体函数而定。
欧拉函数
若 \(i\) 为质数,显然 \(\varphi(i)=i-1\)。
设质数 \(\textit{prime}_j\)。
若 \(i\bmod\textit{prime}_j\neq0\),则 \(i\perp\textit{prime}_j\),因为 \(\varphi\) 是一个积性函数,于是有:
否则,当 \(i\equiv0\pmod{\textit{prime}_j}\) 时,有:
因为与 \(i\cdot\textit{prime}_j\) 互质,即与 \(i\) 互质。\(1,2,3,\cdots,i\cdot\textit{prime}_j\) 在模 \(i\) 意义下以 \(i\) 为周期循环了 \(\textit{prime}_j\) 次。
莫比乌斯函数
若 \(i\) 为质数,显然 \(\mu(i)=-1\)。
设质数 \(\textit{prime}_j\)。
若 \(i\bmod\textit{prime}_j\neq0\) 时,则 \(i\perp\textit{prime}_j\),因为 \(\mu\) 是一个积性函数,于是有:
否则,当 \(i\equiv0\pmod{\textit{prime}_j}\) 时,有:
因为 \(i\) 含有 \(\textit{prime}_j\),则 \(i\cdot\textit{prime}_j\) 至少含有两个 \(\textit{prime}_j\),故 \(\mu\left(i\cdot\textit{prime}_j\right)=0\)。
参考代码
//ans1 为欧拉函数,ans2 为莫比乌斯函数
void pre(){
static int vis[N+1],prime[N+1],size;
ans1[1]=ans2[1]=1;
for(int i=2;i<=N;i++){
if(!vis[i]){
prime[++size]=i;
vis[i]=i;
ans1[i]=i-1;
ans2[i]=-1;
}
for(int j=1;j<=size&&i*prime[j]<=N;j++){
vis[i*prime[j]]=prime[j];
if(i%prime[j]==0){
//ans2 为 0,可以不写
ans1[i*prime[j]]=ans1[i]*prime[j];
break;
}
ans1[i*prime[j]]=ans1[i]*ans1[prime[j]];
ans2[i*prime[j]]=-ans2[i];
}
}
//前缀和,可不写
for(int i=1;i<=N;i++){
ans1[i]+=ans1[i-1];
ans2[i]+=ans2[i-1];
}
return;
}
杜教筛
杜教筛用于快速求解积性函数 \(f\) 的前缀和 \(S(n)=\sum\limits_{i=1}^nf(i)\)。(实际上,只要能够构造恰当的函数,均可以使用杜教筛求解。)
构造一组 \(h=f*g\),有:
可得:
因此,只要 \(h,g\) 适当,能够快速地求出来,就能够在较短时间内求出 \(S(n)\)。
数论分块
可以发现,\(\sum\limits_{i=2}^ng(i)S\left(\left\lfloor\dfrac ni\right\rfloor\right)\) 直接计算的复杂度较高,而 \(\left\lfloor\dfrac ni\right\rfloor\) 就很适合用于数论分块来加速。
记忆化
每当我们求出了 \(S(n)\) 之后,都应当将其记录下来,避免重复计算。
这样的杜教筛的时间复杂度时 \(\mathcal O\left(n^{\frac34}\right)\) 的。
线性筛预处理
可以预处理较小的一部分的 \(S(n)\),从而节约时间。
当预处理前 \(\mathcal O\left(n^\frac23\right)\) 的 \(S(n)\) 时,时间复杂度最小,为 \(\mathcal O\left(n^\frac23\right)\)。
实际常数影响
在 OI 代码中,递归求解 \(S(n)\) 会有常数的影响。
因此最好的方法应当是实际测试预处理部分大小 \(m\) 的值从而确定时间最小的 \(m\),这样即使复杂度劣一些,使用时间也小一些。
参考代码
//#include<bits/stdc++.h>
#include<algorithm>
#include<iostream>
#include<cstring>
#include<iomanip>
#include<cstdio>
#include<string>
#include<vector>
#include<cmath>
#include<ctime>
#include<deque>
#include<queue>
#include<stack>
#include<list>
#include<unordered_map>
#include <ext/pb_ds/assoc_container.hpp>
#include <ext/pb_ds/tree_policy.hpp>
using namespace std;
typedef long long ll;
__gnu_pbds::gp_hash_table<int,ll>mem1,mem2;
constexpr const int N=5e6/*1664510*/;
ll ans1[N+1],ans2[N+1];
ll S1(int n){
if(n<=N){
return ans1[n];
}
if(mem1[n]){
return mem1[n];
}
ll ans=n*(n+1ll)>>1;
for(ll l=2,r=n;l<=n;l=r+1,r=n){
ll t=n/l;
r=n/t;
ans-=(r-l+1ll)*S1(t);
}
return mem1[n]=ans;
}
ll S2(int n){
if(n<=N){
return ans2[n];
}
if(mem2[n]){
return mem2[n];
}
ll ans=1;
for(ll l=2,r=n;l<=n;l=r+1,r=n){
ll t=n/l;
r=n/t;
ans-=(r-l+1ll)*S2(t);
}
return mem2[n]=ans;
}
void pre(){
static int vis[N+1],prime[N+1],size;
ans1[1]=ans2[1]=1;
for(int i=2;i<=N;i++){
if(!vis[i]){
prime[++size]=i;
vis[i]=i;
ans1[i]=i-1;
ans2[i]=-1;
}
for(int j=1;j<=size&&i*prime[j]<=N;j++){
vis[i*prime[j]]=prime[j];
if(i%prime[j]==0){
ans1[i*prime[j]]=ans1[i]*prime[j];
break;
}
ans1[i*prime[j]]=ans1[i]*ans1[prime[j]];
ans2[i*prime[j]]=-ans2[i];
}
}
for(int i=1;i<=N;i++){
ans1[i]+=ans1[i-1];
ans2[i]+=ans2[i-1];
}
return;
}
int main(){
/*freopen("test.in","r",stdin);
freopen("test.out","w",stdout);*/
ios::sync_with_stdio(false);
cin.tie(0);cout.tie(0);
pre();
ll T;
cin>>T;
while(T--){
ll n;
cin>>n;
mem1.clear();
mem2.clear();
cout<<S1(n)<<' '<<S2(n)<<'\n';
}
cout.flush();
/*fclose(stdin);
fclose(stdout);*/
return 0;
}
群论
群论是近代数学的基础之一。
群即对称,对称即群。
“群”是一种代数结构,群 \((G,\circ)\) 由一个集合 \(G\) 和一个二元运算符 \(\circ\) 组成。
群 \((G,\circ)\) 需要满足以下性质:
-
封闭性。即对于 \(\forall a,b\in G\),有 \(a\circ b\in G\)。
-
结合律。即对于二元运算 \(\circ\) 满足 \(a\circ(b\circ c)=(a\circ b)\circ c\)。
-
存在单位元。即存在 \(e\in G\) 满足对于 \(\forall a\in G\) 有 \(e\circ a=a\circ e=a\)。
例如在 \(\circ\) 为乘法的情况下,单位元 \(e\) 为 \(1\)。
-
存在逆元。即对于 \(\forall a\in G\),存在 \(x\in G\) 满足 \(a\circ x=x\circ a=e\)。
记 \(x\) 为 \(a^{-1}\)。
同时,记:
有时为了强调单位元,也会将群 \((G,\circ)\) 记作群 \((G,\circ,e)\)。
例如这两个群:
-
群 \((\Z,+,0)\)。
显然满足上文所述性质。
无限群与整数集加法群同构。
-
群 \((\lbrace0,1,2,\cdots,p-1\rbrace,+_{\bmod p})\),其中 \(+_{\bmod p}\) 表示模 \(p\) 意义下的加法。
显然也满足上文所述性质。
有限群与模意义下的整数集加法群同构。
需要注意的是,群不一定需要满足交换律。
群的性质
已知 \((G,\circ,e)\) 是一个群。
单位元唯一
假设存在两个不同的单位元 \(e_1,e_2\in G\)。
那么有 \(e_1=e_1\circ e_2=e_2\),则 \(e_1=e_2\)。
故,单位元唯一。
消去律
- 左消去律:\(a\circ b=a\circ c\Rightarrow b=c\)。
- 右消去律:\(a\circ b=c\circ b\Rightarrow a=c\)。
逆元唯一
假设 \(a\) 存在两个不同的逆元 \(a_1,a_2\)。
则有 \(a\circ a_1=a\circ a_2=e\)。
由消去律有 \(a_1=a_2\)。
故,逆元唯一。
逆元的拆解
对于 \(\forall a,b\in G\),有 \((a\circ b)^{-1}=b^{-1}\circ a^{-1}\)。
由群的性质,可以得到:
因此,有 \((a\circ b)^{-1}=b^{-1}\circ a^{-1}\)。
群的有关概念
已知 \((G,\circ,e)\) 是一个群。
有限群
当 \(G\) 为有限集时。
有限群与模意义下的整数集加法群同构。
群的阶
群 \((G,\circ,e)\) 的阶即 \(\vert G\vert\)。
无限群
当 \(G\) 为无限集时。
无限群与整数集加法群同构。
阿贝尔群
若二元运算 \(\circ\) 满足交换律,即满足 \(a\circ b=b\circ a\),则称群 \((G,\circ,e)\) 为阿贝尔群、可交换群或加法群。
子群
设集合 \(H\in G\) 且 \(\vert H\vert>0\)。
如果满足 \((H,\circ,e)\) 为一个群,则称群 \((H,\circ,e)\) 是群 \((G,\circ,e)\) 的子群,记作 \(H\leq G\)。
平凡子群
对于任意群 \((G,\circ,e)\),总有平凡子群 \((G,\circ,e)\) 和 \((\lbrace e\rbrace,\circ,e)\)。

浙公网安备 33010602011771号