Skip to main content

数学知识

质数

tip

什么是质数? 在大于1的整数中,如果只包含1和本身这两个约数,就被称为质数(素数)

质数的判定—试除法

从定义出发的写法 O(n)

bool is_prime(int n)
{
if(n<2) return false;
for(int i=2;i<n;i++)
if(n % i == 0) return false;
return true;
}
tip

质数的性质:

如果 d | n, 那么 n/d | n; 所以我们在枚举的时候可以只枚举d <= r/d 的数

优化后的代码 O(n^(1/2))

bool is_prime(int n)
{
if(n<2) return false;
for(int i=2;i<=n/i;i++) // 不会溢出,且速度最快的写法
if(n % i == 0) return false;
return true;
}

分解质因数—试除法

tip

暴力思路:从小到大枚举所有质因数,如果i等于0就输出;

优化思路:和优化判断质因数的方法一样

void divide(int n)
{
for(int i=2;i<=n/i;i++)
{
int s = 0;
while(n % i == 0)
{
n /= i;
s++;
}
if(s) printf("%d %d\n",i,s);
}
if(n > 1) printf("%d %d\n",n,1);
puts("");
}

筛法

tip

信安数论的内容:埃式筛法,如果你想要求1-n中所有的质数那么就从1-n中剔除掉1-n\sqrt{n}中所有质数的倍数,剩下的就都是质数啦;

线性筛法

每次都用i的最小质因子来进行筛除即可,如果不是最小质因子就break掉 ;

int primes[N],cnt;
bool st[N]; // 记录i是否被筛除掉
// 埃式筛法
void get_primes(int n)
{
for(int i=2;i<=n;i++)
{
// 在把这个数的倍数给筛除掉, 也就是说这个数是一个质数
if(!st[i])
{
primes[cnt++] = i; // 记录一下素数
// 将这个质数的所有质数删除掉
for(int j = 2*i;j <= n/i;j+=i) st[j] = true;
}
}
}

// 线性筛法
void get_primes2(int n)
{
for(int i=2;i<=n;i++)
{
if(!st[i]) primes[cnt++] = i;
// 枚举所有的质数
for(int j=0;primes[j]<=n/i;j++)
{
// 把当前质数和primes[j]的乘积给筛除掉
st[primes[j]*i] = true;
if(i % primes[j] == 0) break;
}
}
}

约数

试除法求一个数所有的约数

tip

和试除法求质数是同样的道理,只枚举较小的约数就可以了,较大的约数可以直接算 时间复杂度:O(nn)O(n\sqrt{n})

vector<int> get_divisors(int n)
{
vector<int> res;
for(int i=1;i<=n/i;i++)
{
if(n % i == 0)
{
res.push_back(i);
if(i != n / i) res.push_back(n/i); // 特判一下平方的情况
}
}
sort(res.begin(),res.end());
return res;
}

约数个数

tip

约数个数可以使用一个定理进行计算: 假设数x可以写作:p1x1p2x2...pkxkp_1^{x_1}p_2^{x_2}...p_k^{x_k}, 那么x约数的个数就是(x1+1)(x2+1)(xn+1)(x_1+1)(x_2+1)\cdot \cdot \cdot (x_n+1)

typedef long long LL;
const int mod = 1e9 + 7;
int main()
{
int n;
cin >> n;
unordered_map<int,int> primes;
while(n--)
{
int x;
cin >> x;
for(int i = 2; i<= x/i;i++)
{
while(x % i == 0)
{
x /= i;
primes[i] ++;
}
}
if(x > 1) primes[x] ++;
}
LL res = 1;
for(auto prime:primes) res= res * (prime.second + 1) % mod;
}

约数之和

tip

约数之和也可以使用一个定理进行计算:

假设数x可以写作:p1x1p2x2...pkxkp_1^{x_1}p_2^{x_2}...p_k^{x_k}, 那么x约数之和就是:

(p10+p21+...+p1x1)(pk0+pk1+...+pkxk)(p_1^0 + p_2^1 + ... + p_1^{x_1}) \cdot \cdot \cdot (p_k^0+p_k^1+...+p_k^{x_k})

// 前面的部分都是和求约数个数一样的,只有在最后枚举的时候不一样
for(auto prime:primes)
{
int p = prime.first, a = prime.second;
LL t = 1;
while(a--) t= (t*p+1) % mod; // 秦九昭公式
res = res * t % mod;
}

欧几里得算法(辗转相除法)

tip

y总的代码太简洁了,我佩服的五体投地

int gcd(int a,int b)
{
return b ? gcd(b, a%b) : a;
}

欧拉函数

tip

什么是欧拉函数?

φ(n)\varphi(n)表示1-n中与n中的互质的数的个数

tip

如何求欧拉函数?

想要求φ(n)\varphi(n), 我们首先需要将数进行质因数分解,然后用如下的公式进行计算:

φ(n)=N(11p1)(11p2)(11pn))\varphi(n) = N(1-\frac{1}{p_1})(1-\frac{1}{p_2}) \cdot \cdot \cdot (1-\frac{1}{p_n}))
note

时间复杂度 O(n)O(\sqrt{n})

int a;
cin >> a;
int res = a;
for(int i=2;i<=a/i;i++)
{
if(a%i == 0)
{
res = res / i * (i - 1);
while(a % i == 0) a /= i;
}
if(a>1) res = res / a * (a - 1);
cout << res << endl;
}

欧拉筛法

tip

在计算线性筛法的时候顺便把欧拉函数都求出来(doge)

typedef long long LL;
int primes[N],cnt;
int phi[N];
bool st[N];

LL get_eulers(int n)
{
phi[1] = 1;
for(int i=2;i<=n;i++)
{
if(!st[i])
{
primes[cnt++] = i;
phi[i] = i-1; // 情况1
}
for(int j=0;j<=n/i;j++)
{
st[primes[j]*i] = true;
if(i % primes[j] == 0)
{
phi[primes[j]*i] = primes[j] * phi[i]; // 情况2
break;
}else{
phi[primes[j]*i] = phi[i] * (primes[j] - 1); // 情况3
}
}
}
LL res = 0;
for(int i=1;i<=n;i++) res += phi[i];
return res;
}

快速幂

快速幂代码模板

tip

快速的求出 ak mod pa^k \ mod \ p 的结果

算法思路见代码,还是比较好理解的

int qmi(int a,int k,int p)
{
int res = 1;
while(k)
{
if(k & 1) res = (LL)res * a % p;
k >>= 1;
a = (LL)a * a % p;
}k
return res;
}

快速幂求逆元

tip

比如你要计算aa1a a^{-1} mod p中a的逆元那么就等价于你计算: ap2a^{p-2} mod p

int res = qmi(a,p-2,p);
if(a % p != 0) cout << res;
else cout << "impossible";

扩展欧几里得算法

tip

裴蜀定理:

对于任意一对正整数a, b, 那么存在非零整数x, y, 使得ax + by = (a,b)

欧几里得定理就是为了求出左边等式的系数

// 扩展欧几里得算法 ax + by = (a,b)
int exgcd(int a,int b,int &x,int &y)
{
if(!b)
{
x = 1,y = 0;
return a;
}
int d = exgcd(b,a%b,y,x);
y -= a / b * x;
return d;
}

求解线性同余方程

tip

求解形如ax = b(mod m)的方程

使用扩展欧几里得算法进行求解:

ax = b(mod m) => ax = b + my => ax + m(-y) = b, 此时b如果是gcd(a,m)的倍数,则有解,否则无解。

tip

如何使用扩展欧几里得定理求解?

由欧几里得定理我们可以得到:

ax0+by0=gcd(a,b)=tax_0 + by_0 = gcd(a,b) = t

我们需要求解的目标是:

ax+my=bax + my = b

上下两个式子相除可以得到(y前面的系数有m在取余的时候就没有了):

x=x0bt(mod m)x = x_0 \cdot \frac{b}{t} (mod \ m)

int a,b,m;
int x,y;
int d = exgcd(a,m,x,y);
if(b % d) puts("impossible");
else printf("%d\n",(LL)x * (b / d) % m);

求解逆元

tip

b = 1 的时候就是求解逆元了

int inverseElement(int a,int mod)
{
int x,y;
int t = exgcd(a,mod,x,y);
t = (t%mod + mod)%mod;
}

中国剩余定理

tip

给定两两互质的数,可以解决一个线性同余方程组:

{x=a1(modm1)x=a2(modmn).x=ak(modmk)\left\{\begin{aligned} x &=a_{1}\left(\bmod m_{1}\right) \\ x &=a_{2}\left(\bmod m_{n}\right) \\ & . \\ x &=a_{k}\left(\bmod m_{k}\right) \end{aligned}\right.

对于这样的同余方程组有一个通用的公式:

设: M = m1m2...mk , Mi = M / mi (不是除法,是逆元运算)

那么

x=a1M1M11+a2M2M21+akMkMk1\begin{aligned} x=& a_{1} \cdot M_{1} \cdot M_{1}^{-1}+a_{2} \cdot M_{2} \cdot M_{2}^{-1} \\ & \cdots+a_{k} \cdot M_{k} \cdot M_{k}^{-1} \end{aligned}

求解逆元可以用扩展欧几里得算法来求逆: 求解 ax = 1 mod m

tip

但是我们在写代码的时候还是采用推导然后进行合并的一个方式

#include <iostream>
#include <algorithm>

using namespace std;

typedef long long LL;

LL exgcd(LL a,LL b,LL &x,LL &y)
{
if(!b)
{
x = 1;
y = 0;
return a;
}
LL d = exgcd(b,a%b,y,x);
y -= a/b*x;
return d;

}

int main()
{
int n;
cin >> n;
LL a1,m1;

bool has_answer = true;
cin >> a1 >> m1; // 取出第一个等式
for(int i=1;i<n;i++)
{
LL a2,m2;
cin >> a2 >> m2;
LL k1,k2;
LL d = exgcd(a1,a2,k1,k2);
if((m2-m1)%d) // 判断是否有解
{
has_answer = false;
break;
}
// 对合并后的等式中的a和m进行更新
k1 *= (m2-m1) / d;
LL t = a2/d;
k1 = (k1 % t + t) % t;

m1 = a1 * k1 + m1;
a1 = abs(a1 * a2 / d);

}
if(has_answer)
{
cout << (m1%a1 + a1) % a1 << endl;
}
else puts("-1");
return 0;
}

高斯消元

tip

高斯消元可用于求解多元方程组

算法步骤:

  1. 枚举每一列c
    1. 找到绝对值最大的一行
    2. 把这一行换到最上面去
    3. 将该行第一个数变为1
    4. 将下面所有行的第c列消成0
  2. 逆过来将上半个三角也消掉
tip

初等行列变换

  1. 把某一行成一个非零的数
  2. 交换某两行
  3. 把某行的若干倍加到另一行
#include <iostream>
#include <algorithm>
#include <cmath>
using namespace std;

const int N = 110;
const double eps = 1e-8;


int n;
double a[N][N];

int gauss()
{
int c,r;
for(c=0,r=0;c<n;c++)
{
int t=r;
// 找到绝对值最大的一行
for(int i = r;i<n;i++)
if(fabs(a[i][c]) > fabs(a[t][c]))
t = i;
if(fabs(a[t][c]) < eps) continue;
// 将绝对值最大的行换到第一行
for(int i=c;i<=n;i++) swap(a[t][i],a[r][i]);
// 将这行的第一个数变为1
for(int i=n;i>=c;i--) a[r][i] /= a[r][c];
// 将下面所有行的第c列消成0
for(int i=r+1;i<n;i++)
if(fabs(a[i][c]) > eps)
for(int j=n;j>=c;j--)
a[i][j] -= a[r][j] * a[i][c];
r++;
}
if(r < n)
{
for(int i=r;i<n;i++)
if(fabs(a[i][n]) > eps)
return 2; // 无解
return 1; // 有无穷多组解
}
// 倒着将方程消一下
for(int i=n-1;i>=0;i--)
for(int j=i+1;j<n;j++)
a[i][n] -= a[j][n] * a[i][j];
return 0; // 有唯一解
}

int main()
{
cin >> n;
for(int i=0;i<n;i++)
for(int j=0;j<n+1;j++)
cin >> a[i][j];
int t = gauss();
if(t == 0)
{
for(int i = 0;i<n;i++)
{
if(fabs(a[i][n]) < eps) cout << "0.00" <<endl;
else printf("%.2lf\n", a[i][n]);
}
}
else if(t == 1) puts("Infinite group solutions");
else puts("No solution");
return 0;
}

组合数

组合数1

tip

根据公式:Cab=Ca1b+Ca1b1C_a^b = C_{a-1}^b + C_{a-1}^{b-1} 对组合数进行计算,使用递推的方法得到;

tip

使用情况,当查询次数非常多,且a和b大小不大的情况

#include <iostream>

using namespace std;
typedef long long LL;
const int N = 2010, mod = 1e9 + 7;
LL c[N][N];

void init()
{
for(int i=0;i<N;i++)
for(int j=0;j<=i;j++)
{
// 特判一下
if(!j) c[i][j] = 1;
// 使用公式进行计算
else c[i][j] = (c[i-1][j] + c[i-1][j-1]) % mod;
}
}

int main()
{
int n;
cin >> n;
// 把0-N的方阵全部都初始化好
init();
while(n--) // 执行n轮查询
{
int a,b;
cin >> a >> b;
cout << c[a][b] << endl;
}
return 0;
}

组合数2

tip

根据如下公式进行计算:

Cab=a!b!(ab)!C_a^b = \frac{a!}{b!(a-b)!}

ps:这里的除法用逆元的形式进行求解,使用快速幂去求逆元;

tip

适用情况:a和b比较大的时候, 但是询问的次数并不是很多的时候

#include <iostream>
#include <algorithm>
using namespace std;

typedef long long LL;
const int N = 1e5 + 10, mod = 1e9 + 7;


LL fact[N],infact[N];

LL qmi(int a,int k,int p)
{
LL res = 1;
while(k)
{
if(k & 1) res = (LL)res * a % p;
a = (LL)a * a % p;
k >>= 1;
}
return res;
}

int main()
{
fact[0] = infact[0] = 1;
for(int i=1;i<N;i++)
{
fact[i] = fact[i-1] * i % mod;
infact[i] = (LL)infact[i-1] * qmi(i,mod-2,mod) % mod;
}
int n;
cin >> n;
while(n--)
{
int a,b;
cin >> a >> b;
cout << (LL)fact[a] * infact[a-b] % mod * infact[b] % mod << endl;
}
}

组合数3

tip

使用Lucas定理进行求解

Cab=Ca mod pb mod pCa/pb/p(mod p)C_a^b = C_{a \ mod \ p}^{b \ mod \ p} \cdot C_{a/p}^{b/p} (mod \ p)
tip

适用情况: 询问的次数很少,但是询问的数据范围非常非常大

#include <iostream>
#include <algorithm>

using namespace std;
typedef long long LL;
int p;

int qmi(int a,int k)
{
int res = 1;
while(k)
{
if(k & 1) res = (LL)res * a % p;
a = (LL)a * a % p;
k >>= 1;
}
return res;
}

int C(int a,int b)
{
if(b > a) return 0; // 边界条件
int res = 1;
// 根据定义 a! / b!(a-b)! 进行计算
for(int i=1,j=a;i<=b;i++,j--)
{
res = (LL)res * j % p;
res = (LL)res * qmi(i,p-2) % p;
}
return res;
}

int lucas(LL a, LL b)
{
if(a<p and b <p) return C(a,b);
// 递归计算卢卡斯定理
return (LL)C(a%p,b%p) * lucas(a/p,b/p) % p;
}

int main()
{
int n;
cin >> n;
while(n--)
{
LL a,b;
cin >> a >> b >> p;
cout << lucas(a,b) << endl;
}
return 0;
}

组合数4

tip

高精度+

#include<iostream>
#include<algorithm>
#include<vector>

using namespace std;

const int N=5010;

int primes[N],cnt;
int sum[N];
bool st[N];

void get_primes(int n)
{
for(int i=2;i<=n;i++)
{
if(!st[i])primes[cnt++]=i;
for(int j=0;primes[j]*i<=n;j++)
{
st[primes[j]*i]=true;
if(i%primes[j]==0)break;//==0每次漏
}
}
}

// 对p的各个<=a的次数算整除下取整倍数
int get(int n,int p)
{
int res =0;
while(n)
{
res+=n/p;
n/=p;
}
return res;
}

//高精度乘
vector<int> mul(vector<int> a, int b)
{
vector<int> c;
int t = 0;
for (int i = 0; i < a.size(); i ++ )
{
t += a[i] * b;
c.push_back(t % 10);
t /= 10;
}
while (t)
{
c.push_back(t % 10);
t /= 10;
}
// while(C.size()>1 && C.back()==0) C.pop_back();
return c;
}

int main()
{
int a,b;
cin >> a >> b;
get_primes(a);

for(int i=0;i<cnt;i++)
{
int p = primes[i];
sum[i] = get(a,p)-get(a-b,p)-get(b,p);//是a-b不是b-a
}

vector<int> res;
res.push_back(1);

for (int i = 0; i < cnt; i ++ )
for (int j = 0; j < sum[i]; j ++ )//primes[i]的次数
res = mul(res, primes[i]);

for (int i = res.size() - 1; i >= 0; i -- ) printf("%d", res[i]);
puts("");

return 0;
}