NOIp 数学基础
目录
- 素数判断
- gcd/exgcd
- 逆元
- crt（中国剩余定理）
- 线性代数
-快速幂- 矩阵相关
- 高斯消元
- 组合数学
- 卡特兰数
- 排列组合
- Lucas定理
- 二项式定理
- 康托展开
- 容斥原理
- 进制相关
- 高精度运算
一、素数判断
1、朴素的埃氏筛法
bool isp[MAXN];
int p[MAXN],tot;
void Prime(){
for(int i=0;i<=n;i++)
isp[i]=1;
isp[0]=isp[1]=0;
for(int i=2;i<=n;i++){
if(isp[i])
p[++tot]=i;
for(int j=i*2;j<=n;j+=i)
isp[j]=0;
}
}
可以发现一个合数被判断的次数非常多，所以有如下优化：
2、欧拉筛（线性筛素数）
bool notp[MAXN];
int p[MAXN],tot;
void Prime(){
memset(notp,0,sizeof notp);
notp[0]=notp[1]=1;
for(int i=2;i<=n;i++){
if(!notp[i])
p[++tot]=i;
for(int j=1;j<=tot;j++){
if(p[j]*i>n)break;
notp[p[j]*i]=1;
if(i%p[j]==0)break;
}
}
}
保证了每个合数只被它最小的质因数筛去，提高了时间效率。
二、gcd、exgcd
1、(gcd)
(gcd(a,b))表示(a)与(b)的最大公因数
(lcm(a,b))表示(a)与(b)的最小公倍数
性质
① (gcd(a,b)=gcd(a-b,b)=gcd(a\%b,b))
证明：设(a>b)
设 (gcd(a-b,b)=g)
(∴g|a-b,g|b)
(∵a=a+(a-b))
(∴g|a)
(∴g|gcd(a,b))
(∴gcd(a-b,b)|gcd(a,b))
同理可证，(gcd(a,b)|gcd(a-b,b))
(∴gcd(a,b)=gcd(a-b,b))
`倒推可知，(gcd(a,0)=gcd(a,a)=a)
②(gcd(a,b)×lcm(a,b)=a×b)
证明：设(g=gcd(a,b)，a=k1 * g , b=k2 * g)
(lcm(a,b)=lcm(k1×g,k2×g))
(=g * lcm(k1,k2)=g×k1×k2)
(gcd(a,b)×lcm(a,b)=g×g×k1×k2)
(=(g×k1)×(g×k2)=a×b)
性质②也可通过唯一分解定理证明。
int Gcd(int a, int b){
return b==0?a:gcd(b,a%b);
}
int Lcm(int a, int b){
return(int)((long long)a*b/gcd(a,b));
}
2.(exgcd)
引理：贝祖定理（裴蜀定理）
对于任意整数(a，b)和它们的最大公约数(d)，关于(x，y)的方程(ax+by=m)当且仅当(m)是(d)的倍数时有整数解
证明（来自这位dalao）：
当(b=0)时，(gcd（a，b）=a)，显然存在一对整数解(x=1，y=0)
若(b！=0)
设(ax_1+by_1=gcd(a,b)=gcd(b,a mod b )=bx_2+(a mod b)y_2)
又因 (a mod b=a-a/b*b)
则 (ax_1+by_1=bx_2+(a-a/b*b)y_2)
(ax_1+by_1=bx_2+ay_2-a/b*by_2)
(ax_1+by_1=ay_2+bx_2-b*a/b*y_2)
(ax_1+by_1=ay_2+b(x_2-a/b*y_2))
解得 (x_1=y_2 , y_1=x_2-a/b*y_2)
因为当(b=0) 时存在 (x , y) 为最后一组解
而每一组的解可根据后一组得到
所以第一组的解 (x , y) 必然存在
证毕
代码如下
void Exgcd(int a,int b,int &x,int &y){
if(!b){
x=1;
y=0;
return;
}
Exgcd(b,a%b,y,x);
y=y-a/b*x;
}
(Exgcd)可用于解不定方程。如NOIp2012的 同余方程
void exgcd(int a,int b,int &x,int &y){
if(b==0){
x=1;y=0;return;
}
exgcd(b,a%b,y,x);
y-=a/b*x;
}
int main(){
cin>>a>>b;
exgcd(a,b,x,y);
cout<<(x+b)%b;
return 0;
}
如果(gcd(a,b))不为1(假设为(c))，可以先求出(gcd)为1的情况，再让(x)乘上(c/gcd(a,b))
如何求最小正整数解呢？
还是一样，只是把((x+b))中的(b)换成(b/gcd(a,b))即可
三、逆元
前备知识
1、模
(a≡b(mod m)) 表示在模m意义下，a等于b。
在模m意义下，两个小于m的数a,b四则运算如下：
加法：((a+b) \% m)
减法：((a-b+ m ) \%m)
乘法：(a * b \% m)
除法—— 就是逆元啦
2、欧拉定理
(Φ(n))表示1-n中与n互质的数的个数
性质：
①对于质数p，Φ(p)=p-1
② $ Φ(n)=Σ^n_{i=1}[(i,n)==1]=n∏(1-frac{1}{p_i})$
(其中pi是n的质因子)
③ (Σ_{d|n}Φ(d)=n)
也就是说，一个数所有因数的Φ值的和等于这个数本身。
费马小定理：如果(a,n)=1，即a与n互质，那么(a^{Φ(n)}=1(mod n))
所以:(a^{Φ(n)-1}*a^1=1(mod n))
由性质①可知，当n是质数p时，对于任意不是p的倍数的a，都有:(frac{1}{a}=a^{p-2})
我们称1/a为a在模n意义下的逆元，经常记作(inv(a))。
顺便提一下，欧拉筛素数的时候可以把(Φ(n))也求出来
phi[1]=1;
memset(check,0,sizeof(check));
for(int i=2;i<=n;i++){
if(!check[i]){
p[++tot]=i;
phi[i]=i-1;
}
for(int j=1;j<=tot;j++){
if(i*p[j]>n)
break;
check[i*p[j]]=1;
if(i%p[j]==0){
phi[i*p[j]]=phi[i]*p[j];
break;
}
else
phi[i*p[j]]=phi[i]*(p[j]-1);//积性函数的性质
}
}
3、求逆元的n种方法
-
扩欧求逆元
见上面例题’同余方程‘
-
快速幂求逆元
(frac{1}{a}=a^{p-2}(mod n))前提是 (n)是质数
- 递推求逆元
void inverse(int n, int p) {
inv[1] = 1;
for (int i=2; i<=n; ++i) {
inv[i]=(ll)(p-p/i)*inv[p%i]%p;
}
}
- 阶乘逆元
已知(n!)的逆元，怎么求((n-1)!)的逆元呢？
很简单，因为(frac{1}{(n-1)!})与(frac{1}{n!})只差分母上的一个(n)，所以只要将(n)的逆元(*n)，即可得到(n-1)的逆元
一个小小的应用——排列组合
(C^m_n = frac{n!}{m!(n-m)!})
用逆元可以优化计算过程
代码如下：
#include<cstdio>
#include<iostream>
#define ll long long;
ll fac[MAXN],inv[MAXN];
int n,mod;
ll pow(ll x,ll n,ll m){
ll res=1;
while( n ){
if( n % 2== 1 ) {
res = res * x;
res = res % m;
}
x = x * x;
x = x % m;
n >>= 1;
}
return res;
}
void pree(){//预处理
fac[0] = inv[0] = 1;
for( int i = 1 ; i <= n ; i++ ){
fac[i] = fac[i - 1] * i % mod;
inv[n] = pow( fac[n] , mod-2 , mod );//快速幂
for( int i = n; i >= 2 ; i-- )
inv[i - 1] = inv[i] * i % mod;
}
}
int C(int n,int m){//求组合数
if( n<0 || m<0 || n-m <0 )
return 0;
return fac[n]*inv[m] % mod *inv[n-m] % mod;
}
四、中国剩余定理
五、快速幂
int pow(int a,int p){//(a^p)%mod
int ret=1;
while(p){
if(p&1)
ret=ret*a%mod;
a=a*a%mod;
p>>=1;
}
return ret;
}
六、矩阵相关
在这里推荐3B1B的线性代数视频，微积分系列的也很不错。
1.定义
(A_{nm})表示一个 (n) 行 (m) 列的矩阵。
用 $a_{ij} 表示矩阵 (A_{nm}) 第 (i) 行第 (j) 列的元素。
特别的，一个 1 行 (n) 列的矩阵，可以被称为行向量，列向量
同理。
一个 (n) 行 (n) 列的矩阵，可以被称为 (n) 阶方阵 (A_n)。
对角矩阵表示除了主对角线上的元素外其他元素全部为 0
的矩阵。
主对角线以下全部是 0 的方阵是上三角矩阵。
单位矩阵是主对角线上的元素全部为 1 的对角矩阵。
矩阵(A_{nm}) 和 (B_{mp}) 相乘可以得到矩阵 (C_{np}) ，且
(c_{ij} =∑^m_{k=1}a_{ik}b_{kj})
一般使用的矩阵乘法是方阵之间相乘，即 (A_{nn}) 乘 (B_{nn}) 。
实现：利用结构体存储矩阵，矩阵乘法使用 O((n^3)) 暴力即可。
2.性质
矩阵乘法有结合律，但没有交换律
3.根据结合律性质，可以仿照快速幂写出矩阵快速幂
#include<bits/stdc++.h>
using namespace std;
#define ll long long
const int MOD=1e9+7;
struct M{
ll a[101][101];
}qwq,e;
int n;ll k;
M cheng(M x,M y){
M c;
for(int i=1;i<=n;i++)
for(int j=1;j<=n;j++)
c.a[i][j]=0;
for(int i=1;i<=n;i++)
for(int j=1;j<=n;j++)
for(int k=1;k<=n;k++)
c.a[i][j]=c.a[i][j]%MOD+x.a[i][k]*y.a[k][j]%MOD;
return c;
}
M pow(M x,ll y){
for(int i=1;i<=n;i++)
e.a[i][i]=1;
M ans=e;
while(y){
if(y&1)ans=cheng(ans,x);
x=cheng(x,x);
y>>=1;
}
return ans;
}
int main(){
scanf("%d%lld",&n,&k);
for(int i=1;i<=n;i++)
for(int j=1;j<=n;j++)scanf("%lld",&qwq.a[i][j]);
M ans=pow(qwq,k);
for(int i=1;i<=n;i++){
for(int j=1;j<=n;j++)
cout<<ans.a[i][j]%MOD<<" ";
cout<<endl;
}
return 0;
}
矩阵快速幂的一个基础运用是加速斐波那契数列的计算。
矩阵快速幂也经常应用在DP优化中。
七、高斯消元
1、行列式
定义:
一个方阵 A 的行列式表示为 (|A|)
(|A| =∑_p(−1)^{σ(p)}
∏^n
_{i=1}
a_{i,p_i})
2、基础性质：
- 一个上三角矩阵的行列式值是所有对角线上元素的乘积。
- 交换矩阵的两行/两列，行列式值取相反数。
- 将矩阵的一行/一列乘上一个固定的常数 k，行列式值也乘
上 k。 - 将矩阵的一行加到另外一行上去，行列式值不变，列同理。
- 存在两行完全一样的矩阵，则行列式值为 0。
3、高斯消元
根据行列式的性质，可以通过高斯消元解多元一次方程组。
理解起来非常容易，就是模拟手算。
bool work(){
for(int i=1;i<=n;i++){
int r=i;
for(int j=i+1;j<=n;j++){
if(fabs(mapp[j][i])>fabs(mapp[r][i]))
r=j;
}
if(fabs(mapp[r][i])<eps)
return 0;
if(r!=i)
swap(mapp[i],mapp[r]);
double tmp=mapp[i][i];
for(int j=i;j<=n+1;j++){
mapp[i][j]/=tmp;
}
for(int j=i+1;j<=n;j++){
tmp=mapp[j][i];
for(int k=i;k<=n+1;k++)
mapp[j][k]-=tmp*mapp[i][k];
}
}
ans[n]=mapp[n][n+1];
for(int i=n-1;i>=1;i--){
ans[i]=mapp[i][n+1];
for(int j=i+1;j<=n;j++)
ans[i]-=(ans[j]*mapp[i][j]);
}
return 1;
}
4.矩阵求逆
由于将一个矩阵变成单位矩阵的过程与将单位矩阵变成该矩阵的逆矩阵的过程相同，我们可以把初始矩阵和单位矩阵放在一起，对初始矩阵通过高斯消元思想变成单位矩阵，变换结束后的’原单位矩阵‘就是初始矩阵的逆矩阵。
#include<iostream>
#include<cstdio>
#include<cstring>
#include<algorithm>
using namespace std;
#define ll long long
const ll p = 1e9 + 7;
int n;
ll f[405][810], tmp;
ll pow( ll x, ll mod ,ll ans ){
while( mod ){
if( mod % 2 ) ans = ans * x % p;
x = x * x % p, mod = mod >> 1;
}
return ans;
}
bool gsxy(){
for( int i = 1; i <= n; i++ ){
if( !f[i][i] ) return 0;
tmp = pow( f[i][i], p-2 ,1 );
for( int j = i + 1; j <= n * 2; j ++ ) f[i][j] = f[i][j] * tmp % p;
f[i][i] = 1;
for( int j = 1; j <=n ; j++ ){
if( j != i ){
for( int k = i + 1; k <= n * 2; k++ )
f[j][k] = ( f[j][k] - f[i][k] * f[j][i] % p + p ) % p;
f[j][i] = 0;
}
}
}
return 1;
}
void print(){
for( int i = 1; i <= n; i ++ ){
for( int j = 1 + n; j <= n * 2; j++ )
printf( "%lld ",f[i][j] );
cout<<endl;
}
}
int main(){
scanf( "%d",&n );
for( int i = 1; i <= n;i++ )
for( int j = 1; j <= n; j++ ){
scanf( "%lld",&f[i][j] );
f[i][i+n] = 1;
}
if( !gsxy() ) cout<<"No Solution"<<endl;
else print();
return 0;
}
八、卡特兰数
(P=a_1×a_2×a_3×……×a_n)，依据乘法结合律，不改变其顺序，只用括号表示成对的乘积，试问有几种括号化的方案？
一个栈(无穷大)的进栈序列为(1，2，3，…，n)，有多少个不同的出栈序列?
在一个凸多边形中，通过若干条互不相交的对角线，把这个多边形划分成若干个三角形，有多少种划分方式？
给定(N)个节点，能构成多少种不同的二叉搜索树？
给定(n)对括号，有多少种括号正确配对的字符串？
令(h(0)=1,h(1)=1)
公式一：
(h(n)= h(0)*h(n-1)+h(1)*h(n-2) + ... + h(n-1)h(0) (n>=2))
公式二：
(h(n)=h(n-1)*(4*n-2)/(n+1))
公式三：
(h(n)=C(2n,n)/(n+1) (n=0,1,2,...))
公式四：
(h(n)=c(2n,n)-c(2n,n-1)(n=0,1,2,...))
九、排列组合
组合数取模问题
计算 (C_n^m(mod P))
(Case I: n,m ≤ 10^3 ,P ≤ 10^9)
(Case II: n,m ≤ 10^5 ,P ≤ 10^9) 且为质数
(Case III: n,m ≤ 10^6 ,P ≤ 10^5)
(Case IV: n,m ≤ 10^9 ,P ≤ 10^5) 且为质数
(Case V: n,m ≤ 10^9 ,P ≤ 10^9) ，(P)的质因数分解
(P =p_1^{e_1}p_2^{e_2}p_3^{e_3}……p_k^{e_k})
满足(∀i = 1,2,...,k)(p_i^{e_i}≤ 10^5)
-
公式(C_n^m=C_{n-1}^{m}+C_{n-1}^{m-1})
-
因为(P)是质数，可以通过上面“逆元”部分的逆元求组合数求得
-
对每个乘数暴力分解质因数，乘除化加减即可
-
应用(Lucas)定理，见下章
5)扩展(Lucas),对(P)每个质因数使用(Lucas)，用中国剩余定理合并。因为我不会所以不讲了。
十、Lucas定理
如果 P 是一个质数，则(C_n^m(mod P)=C(n mod P,m mod p)*C(n/p,m/p) mod P)
ll lucas(ll x,ll y){
if(x<y)return 0;
return x<p?fz[x]*fm[y]*fm[x-y]%p:lucas(x%p,y%p)*lucas(x/p,y/p)%p;
}
void init(){
fz[0]=fm[0]=fz[1]=fm[1]=1;
for(ll i=2;i<=n+m;i++){
fz[i]=fz[i-1]*i%p;
fm[i]=(p-p/i)*fm[p%i]%p;
}
for(ll i=2;i<=n+m;i++)fm[i]=fm[i]*fm[i-1]%p;
}
十一、二项式定理
应该是常识吧...小学就应该学过杨辉三角。
#include<bits/stdc++.h>
using namespace std;
int a,b,n,m,k;
const int MOD=10007;
int pow(int a,int b){
int ans = 1;
while (b){
if (b & 1)ans=ans*a%MOD;
a = a*a%MOD;
b>>=1;
}
return ans;
}
int C[1500][1500];
void init(){
C[1][1]=1;
C[1][0]=1;
for(int i=2;i<=k;i++)C[i][i]=C[i][0]=1;
for(int i=2;i<=k;i++)
for(int j=1;j<=i;j++){
C[i][j]=(C[i-1][j-1]+C[i-1][j])%MOD;
}
}
int main(){
scanf("%d%d%d%d%d",&a,&b,&k,&n,&m);
a%=MOD;
b%=MOD;
init();
int ans=1;
ans=(ans*pow(a,n))%MOD;
ans=(ans*pow(b,m))%MOD;
ans=ans*C[k][m]%MOD;
printf("%d
",ans);
return 0;
}
十二、康托展开
由于这次初赛考了，有点小紧张orz
(X=a[n]*(n-1)!+a[n-1]*(n-2)!+...+a[i]*(i-1)!+...+a[1]*0!)
(a[i])表示(i)之后比(i)小的元素个数。
排列(P)的排名就是(P)的康托展开值(+1)
康托展开求逆(参考这位dalao)：
假设求4位数中第19个位置的数字。
① (19)减去(1 → 18)
② (18) 对 (3!) 作除法 (→) 得(3)余(0)
③ (0)对 (2!) 作除法 → 得(0)余(0)
④ (0)对 (1!) 作除法 → 得(0)余(0)
据上面的可知：
我们第一位数（最左面的数），比第一位数小的数有3个，显然 第一位数为→ 4
比第二位数小的数字有0个，所以 第二位数为→1
比第三位数小的数字有0个，因为1已经用过，所以第三位数为→2
第四位数剩下 3
该数字为 4123 (正解)
void reverse_contor(int x){
memset(vis,0,sizeof vis);
x--;
int j;
for(int i=1;i<=n;i++){
int t=x/fac[n-i];
for(j=1;j<=n;j++){
if(!vis[j]){
if(!t) break;
t--;
}
}
printf("%d ",j);
vis[j]=1;
x%=fac[n-i];
}
puts("");
}
十三、容斥原理
一个比较玄学的专题，我也不是很会。放一些基础理论吧。
考到就挂了
十四、进制相关
发现自己进制专题比较薄弱，NOIp2000的一道黄题也不会做。当初学二进制与十进制转化的时候没有很好地理解。现在补一下。
我们看到NOIp这道题出现了负的进制，然而c++里模负数和我们想要的结果并不一样。题解中dalao如是说道：
（商+1）* 除数+（余数-除数）=商 * 除数+除数+余数-除数=商*除数+余数=被除数
char jz[20]={'0','1','2','3','4','5','6','7','8','9','A','B','C','D','E','F','G','H','I','J'};
void zh(int x,int y){
if(x==0)return;
else{
if(x>0||x%y==0){
zh(x/y,y);
cout<<jz[x%y];
return;
}
else {
zh(x/y+1,y);
cout<<jz[-y+x%y];
return;}
}
}
int main(){
int n,r;
while(~scanf("%d",&n)){
cin>>r;
cout<<n<<"=";
zh(n,r);
cout<<"(base"<<r<<")"<<endl;
}
return 0;
}
十五、高精度运算
前备知识：重载运算符
C++ 允许在同一作用域中的某个函数和运算符指定多个定义，分别称为函数重载和运算符重载。
重载声明是指一个与之前已经在该作用域内声明过的函数或方法具有相同名称的声明，但是它们的参数列表和定义（实现）不相同。
当您调用一个重载函数或重载运算符时，编译器通过把您所使用的参数类型与定义中的参数类型进行比较，决定选用最合适的定义。选择最合适的重载函数或重载运算符的过程，称为重载决策。
也就是说，当我们给高精度数定义加减乘除时，由于操作的是结构体，所以对之前默认存在的数字加减乘除不会产生影响，二者并存。
拿自定义高精加举例：
data operator + (const data &a, const data &b)
data是结构体的类型，在重载的+之前需要加上operator(记住这个单词qwq)
struct data
{
int a[1000], len;
};//用结构体内数组储存高精度数字，len储存长度
data operator + (const data &a, const data &b)
{//重载运算符
data c;
c.len = max(a.len, b.len);//得数的最小位数
for (int i = 0; i <= c.len; i++)
c.a[i] = a.a[i] + b.a[i];//逐位相加。
//注意，此时可能出现一位上超过十的情况
for (int i = 0; i <= c.len; i++)
{//进位
c.a[i + 1] += c.a[i] / 10;
c.a[i] %= 10;
}
if (c.a[c.len + 1]) c.len++;//判断是否最高位有进位
return c;
}
data operator - (const data &a, const data &b)
{
data c;
c.len = a.len;//前提是a>b，得数最大的位数是被减数位数
for (int i = 0; i <= c.len; i++)
c.a[i] = a.a[i] - b.a[i];//逐位相减
for (int i = 0; i <= c.len; i++)
if (c.a[i] < 0)
{//退位
c.a[i + 1]--;
c.a[i] += 10;
}
while (c.a[c.len] == 0) c.len--;
return c;
}
data operator * (const data &a, const int &b)
{//高精乘int
data c;
c.len = a.len;//至少是高精数的位数
for (int i = 0; i <= c.len; i++)
c.a[i] = a.a[i] * b;
for (int i = 0; i <= c.len + 10; i++)
//注意：此时每一位上的数字可能大于20，甚至非常大
{
c.a[i + 1] += c.a[i] / 10;
//进位
c.a[i] %= 10;
}
for (int i = c.len + 10; i >= 0; i--)
if (c.a[i])
{
c.len = i;
break;
}//寻找最高位
return c;
}
data operator * (const data &a, const data &b)
{//高精乘高精
data c;
c.len = a.len + b.len;//划重点qwq
for (int i = 0; i <= a.len; i++)
for (int j = 0; j <= b.len; j++)
c.a[i + j] += a.a[i] * b.a[j];
for (int i = 0; i <= c.len + 10; i++)
{
c.a[i + 1] += c.a[i] / 10;
c.a[i] %= 10;
}
for (int i = c.len + 10; i >= 0; i--)
if (c.a[i])
{
c.len = i;
break;
}
return c;
}
data operator / (const data &a, const int &b)
{
data c;
int d = 0;
for (int i = a.len; i >= 0; i--)
{//从高往低除
d = d * 10 + a.a[i];
c.a[i] = d / b;
d %= b;
}
for (int i = a.len; i >= 0; i--)
if (c.a[i])
{
c.len = i;
break;//寻找最高位
}
return c;
}
//高精除高精...老师说oi不会考（可以尝试二分答案，用猜到的数去乘除数），勉强信了qwq
【例题】
给出正整数a，b，求a^b的最后500位。
a<=10^500 , b<=100
【分析】
分析a的数据范围可知此题需要用高精度。
根据常识，较大位上的数不会对较小位上的数产生影响，所以只需记录最后500位的数即可。（结合快速幂）