数学知识
一、素数（质数）
假如给你以一个正整数N，你能否找到1~n之间所有的素数？素数的脾气，可不好惹。
1.双层循环暴力求值
这个求法没什么好讲的，就是两个for循环嵌套一下遍历所有的数，代码如下。
void Prime(int n)
{
vis[0] = vis[1] = 1;
for(int i = 2;i <= n;i++)
{
for(int j = 2;j < i;i++)
{
if(i % j == 0)vis[i] == 1;
}
}
}
如此一来，N以内的所有素数都被标记成了true，如果想要存在数组里，再写个for循环就可以做到。
2.素数筛
这个是最基础的素数筛，将N以内的所有素数都存到了数组prime里。这个筛法的原理是将所有的数的倍数都标记为ture，那么没有被标记的数i就说明从2~i-1之间没有它的约数，就说明它是质数。
void Prime(int n)
{
vis[0] = vis[1] = 1;
for(int i = 2;i < n;i++)
{
if(!vis[i]) prime[idx++] = i;
for(int j = i + i;j <= n;j += i)
vis[j] = 1;
}
}
3.埃氏筛
埃氏筛的原理是与上边的类似，但是埃氏筛只标记素数的倍数，时间复杂度大大降低。
void Prime(int n)
{
vis[0] = vis[1] = 1;
for(int i = 2;i <= n;i++)
{
if(!vis[i])
{
for(int j = i + i;j <= n;j += i) vis[j] = 1;
}
}
}
可见所有的素数都被标记了false，如果想要存在数组中，只要将一个数组放到if函数里存i即可。
4.欧拉筛（线性筛）
欧拉筛与埃氏筛的区别就是埃氏筛会对一个非素数重复标记，而欧拉筛会且仅访问一次每个数。在时间复杂度上进一步降低。
void Prime(int n)
{
vis[0] = vis[1] = 1;
for(int i = 2;i <= n;i++)
{
if(!vis[i]) prime[idx++] = i;
for(int j = 0;prime[j] <= n / i;j++)
{
vis[prime[j] * i] = 1;
if(i % prime[j] == 0) break;//当前记录的素数最大时i，能够被i整除的只有i本身。
}
}
}
二、约数
1.试除法求约数
当正整数num存在一个约数i的时候，num / i也是num的约数
int idx = 0;
cin >> num;
for(int i = 1;i * i <= num;i++)
{
if(num % i == 0)
{
a[idx++] = i;
if(num / i != i)a[idx++] = num / i;
}
}
如此一来，num的所有约数就储存到了数组a中。
2.约数的个数
输入几个数，求其乘积的约数的个数。我们假设乘积为N = a1 * a2 * a3，易知N的约数数量为a1,a2,a3的约数数量的积，再将三个数进行质因数分解即可得到答案。例如：a1 = p1b1 * p2b2 * p3b3 * ………… * pnbn，此时p1b1、p2b2、p3b3………… pnbn都是a1的约数，它们本身的约数数量为指数加一，将所有的相乘，我们就可以得到a1的约数的数量了。
unordered_map<int,int>primes;
int n;
cin >> c;
while(n--)
{
cin >> num;
for(int i = 2;i * i <= num;i++)
while(num % i == 0)
{
num /= i;
primes[i]++;
}
if(num > 1)primes[num]++;
}
for(auto p : primes)ans = ans * (p.second + 1) % mod;
3.约数之和
输入几个数，求其乘积的约数的和。方法与上边类似，首先将其质因数分解，得到a1 = p1b1 * p2b2 * p3b3 * ………… * pnbn，然后 我们可以知道，a1的约数之和为(p10 + p11 + …… + p1b1) * (p20 + p21 + …… + p2b2) * ……其他的部分于约数的个数计算方法类似
unordered_map<int,int>primes;
int n;
cin >> n;
while(n--)
{
cin >> num;
for(int i = 2;i * i <= num;i++)
while(num % i == 0)
{
num /= i;
primes[i]++;
}
if(num > 1)primes[num]++;
}
for(auto p : primes)
{
int t = 1;
int a = p.first,b = p.second;
while(b--)t = (t * a + 1) % mod;
ans = ans * t % mod;
}
4.最大公约数
利用辗转相除法，可以求得最大公约数。辗转相除法的实现基于gcd(a , b) = gcd(b , a mod b)，假设a与b的最大公约数为d，那么就有a % d == 0, b % d == 0，而a mod b = a - [a / b] * b = a - c * b，可见a mod b与可以被d整除。所以经过辗转相除即可得到一数为0的情况，此时另一个数就是最大公约数。
int gcd(int a,int b)
{
return b ? gcd(b, a % b) : a;
}
三、欧拉函数
1.欧拉函数
欧拉函数用于求自然数N内的所有于N互质的数的数量，即1 ~ N内于N互质的数的数量。
欧拉函数：φ(n) = n * (1 - 1 q 1 \frac{1}{q~1~} q 1 1) * (1 - 1 q 2 \frac{1}{q~2~} q 2 1) * …… * (1 - 1 q k \frac{1}{q~k~} q k 1)
利用欧拉函数就能够求出n内与其互质的自然数的数量
#include<iostream>
using namespace std;
int main()
{
int n;
cin >>n;
int res = n;
for(int i = 2;i * i <= n;i++)
if(n % i == 0)
{
res = res / i * (i - 1);
while(n % i == 0)n /= i;
}
if(n > 1)res = res / n * (n - 1);
cout << res << endl;
return 0;
}
2.用筛法求欧拉函数
我们可以利用线性筛对n之内的所有的欧拉函数进行求和处理。
首先看看线性筛的写法：
for(int i = 2;i <= n;i++)
{
if(!vis[i])primes[idx++] = i;
for(int j = 0;primes[j] <= n / i;j++)
{
vis[primes[j] * i] = 1;
if(i % primes[j] == 0)break;
}
}
下边，我们将对其改造，使其成为能够求得n范围内的欧拉函数的和的形式。
结合欧拉函数的定义，当i为质数时，它与其范围内的所有正整数互质，故其欧拉函数为i-1。
if(!vis[i])
{
primed[idx++] = i;
phi[i] = i - 1;
}
当i mod primes[j] == 0时，说明primes[j]为i的质因数，则phi[primes[j] * i] = phi[i] * primes[j];
当i mod primes[j] != 0时，说明primes[j]不是i的质因数，则phi[primse[j] * i ] = phi[i] * primes[j] * (1 - 1 p r i m e s [ j ] \frac{1}{primes[j]} primes[j]1) = phi[i] * (primes[j] - 1);
if(i % primes[j] == 0)
{
phi[primes[j] * i] = phi[i] * primes[j];
break;
}
phi[primes[j] * i] = phi[i] * (primes[j] - 1);
最后将其相加即可，看看最终的代码：
ll get_euler(int n)
{
phi[1] = 1;
for(int i = 2;i <= n;i++)
{
if(!vis[i])
{
primes[idx++] = i;
phi[i] = i - 1;
}
for(int j = 2;primes[j] * i <= n;j++)
{
int t = i * primes[j];
vis[t] = 1;
if(i % primes[j] == 0)
{
phi[t] = phi[i] * primes[j];
break;
}
primes[t] = phi[i] * (primes[j] - 1);
}
}
ll ans = 0;
for(int i = 1;i <= n;i++)
ans += phi[i];
return ans;
}
四、快速幂
1.低精度的快速幂
1.例如这题：875. 快速幂 - AcWing题库
由数学知识可知：a * b mod m = (a mod m) * (b mod m) mod m。那么就可以将ab mod p转换为若干ak(k = 2n) mod p的乘积。这一操作在二进制上进行，例如：35 mod 2，5的二进制为101，那么就是(3a mod 2) * (3b mod 2) mod 2【a = 20,b = 22】
#define ll long long
ll f(ll a,int b,int c)//a^b%c
{
ll ans = 1;
while(b)
{
if(b&1) ans = ans * a % c;
b >>= 1;
a = a * a % c;
}
return ans;
}
2.那么，如果要求的是最后n位数而不是取模呢？这里假设是三位数。
- 写法一：
ll f(ll a,int b)
{
ll ans = 1;
while(b)
{
if(b&1) ans = ans * a % 1000;
b >>= 1;
a = a * a % 1000;
}
return ans;
}
- 写法二：
ll f(ll a,int b)
{
ll ans = 1;
for(int i = 0;i < b;i++)
{
ans *= b;
ans %= 1000;
}
return ans;
}
2.高精度的快速幂
高精度快速幂常常用于求最后多位数字，例如这题：P1249 最大乘积 - 洛谷 | 计算机科学教育新生态 (luogu.com.cn)
这题涉及到了贪心算法，这里不加讨论，我们所需要知道的就是它的计算方法是从2开始往上逐步增加，最后减去相应的数即可。例如：15就是2+3+4+5+6-5。13就是2+3+4+5-2+1=2+3+4+6。
首先利用头文件cmath内置的log10()来计算数的位数。
log10(2)* p + 1;
它的快速幂分为三个模块，我们一一来看。
1.主体部分：
int ans[N],f[N],sav[N];
while(p)
{
if(p&1)f1();
p >>= 1;
f2();
}
2.函数f1()
void f1()
{
memset(sav, 0, sizeof sav);
for (int i = 1; i <= 500; i++)
for (int j = 1; j <= 500; j++)
sav[i + j - 1] += ans[i] * f[j];
for (int i = 1; i <= 500; i++)
{
sav[i + 1] += sav[i] / 10;
sav[i] %= 10;
}
memcpy(ans, sav, sizeof sav);
}
3.函数f2()
void f2()
{
memset(sav, 0, sizeof sav);
for (int i = 1; i <= 500; i++)
for (int j = 1; j <= 500; j++)
sav[i + j - 1] += f[i] * f[j];
for (int i = 1; i <= 500; i++)
{
sav[i + 1] += sav[i] / 10;
sav[i] %= 10;
}
memcpy(f, sav, sizeof f);
}
4.初始化与结果减1需注意：
//初始化
f[1] = 2;
ans[1] = 1;
…………
ans[1] -= 1;
下面看看完整代码：
#include<iostream>
#include<cmath>
#include<cstring>
using namespace std;
#define ll long long
const int N = 1005;
int p;
int ans[N], f[N], sav[N];
void f1()
{
memset(sav, 0, sizeof sav);
for (int i = 1; i <= 500; i++)
for (int j = 1; j <= 500; j++)
sav[i + j - 1] += ans[i] * f[j];
for (int i = 1; i <= 500; i++)
{
sav[i + 1] += sav[i] / 10;
sav[i] %= 10;
}
memcpy(ans, sav, sizeof sav);
}
void f2()
{
memset(sav, 0, sizeof sav);
for (int i = 1; i <= 500; i++)
for (int j = 1; j <= 500; j++)
sav[i + j - 1] += f[i] * f[j];
for (int i = 1; i <= 500; i++)
{
sav[i + 1] += sav[i] / 10;
sav[i] %= 10;
}
memcpy(f, sav, sizeof f);
}
int main()
{
cin >> p;
cout << int(log10(2) * p + 1) << endl;
ans[1] = 1;
f[1] = 2;
while (p)
{
if (p & 1)f1();
p >>= 1;
f2();
}
ans[1] -= 1;
for (int i = 500; i >= 1; i--)
{
if (i % 50 == 0 && i != 500)cout << endl;
cout << ans[i];
}
return 0;
}
3.指数较大的快速幂
先看题目要求（题目来源：J-大数乘法_哈尔滨理工大学第12届程序设计竞赛（同步赛） (nowcoder.com)）
我们可以看到，y的大小达到了惊人的10100000，所以常规的快速幂写法已经无法解决它了。所以我们可以将指数给拆分掉，分批次将其用运用到快速幂中。
来看看代码的实现：
#include<iostream>
using namespace std;
typedef long long ll;
ll x,p;
string y;
ll qpow(ll a,ll b,ll m)//快速幂模板，没什么好说的
{
ll res = 1;
while(b)
{
if(b & 1)res = res * a % p;
a = a * a % p;
b >>=1;
}
return res;
}
int main()
{
int t;
cin >> t;
while(t--)
{
cin >> x >> y >> p;
ll ans = 1;
for(int i = y.size() - 1;~i;--i)
{
ans = ans * qpow(x,y[i] - '0',p) % p;
x = qpow(x,10,p);//每一位计算完成之后，都将底数进行s
}
cout << ans << endl;
}
return 0;
}
4.快速幂求逆元
若整数 b，m互质，并且对于任意的整数 a，如果满足 b|a，则存在一个整数 x，使得 a b \frac{a}{b} ba ≡ a * x(mod m)，则称 x 为 b 的模 m 乘法逆元，记为 b−1(modm)。
b 存在乘法逆元的充要条件是 b与模数 m 互质。当模数 m 为质数时，bm−2 即为 b 的乘法逆元
然后就是直接套用快速幂的模板即可，返回即为逆元：
f(a, p - 2, p);
5.快速幂与矩阵【华东交通大学教学内容】
5.1矩阵快速幂的实现
在矩阵快速幂中，我们利用类来存储矩阵及其运算法则。
我们先来看看代码的具体实现，再来讲解其语法。
struct Matrix
{
ll a[N][N];
Matrix()//构造函数
{
memset(a,0,sizeof a);
}
void init()//单元矩阵
{
for(int i = 0;i <= n;i++)
a[i][i] = 1;
}
Matrix operator*(const Matrix& B)const//运算符*的重构
{
Matrix C;
for(int i = 0;i < n;i++)
for(int k = 0;k < n;k++)
for(int j = 0;j < n;j++)
C.a[i][j] = (C.a[i][j] + a[i][k] * B.a[k][j]) % MOD;
return C;
}
Matrix operator^(const int& t)const
{
Matrix A = (*this),res;
res.init();
int p = t;
while(p)
{
if(p & 1)res = res * A;
A = A * A;
p >>= 1;
}
return res;
}
};
-
该结构体使用到的陌生语法有：构造函数、运算符重载和this指针，下边我们来看看着两个语法：
-
构造函数：
构造函数的特点是在结构体或者类中，函数名与结构体名称相同且其声明前方没有返回值。它会在结构体创建的时候自动调用，完成对结构体的初始化。与之对应的是析构函数，作用是对结构体的清除，感兴趣的同学可以自行了解。
-
运算符重载：
- 运算发重载的含义：指的是将已有的运算符赋予新的功能，并通过识别使用是的数据类型和数目来判断采用哪种操作。c/c++中有许多运算符其实已经经过了重构，例如*，在将其作用于地址的时候，它的作用是解引用；将其用于两个数字的时候，它的作用是得到两数的乘积。运算符重载的存在，使得代码可以变得更加简洁易懂。
- 在这里，我们将 * 运算符重载出求两矩阵的积的作用，将 ^ 运算符重载出求矩阵的次方的作用。当然，在 ^ 中我们使用了重载之后的 * 运算符。
-
this指针
- this指针存在于成员函数中，是一个用来指向对象本身地址的指针。例如上面代码中，this指针将该指针所在的结构体本身赋值给了结构体A。
-
5.2例题讲解
5.2.1利用矩阵快速幂推导斐波那契数列
在这之前，我们推导斐波那契数列的方法是这样的：
f[1] = f[2] = 1;
for(int i = 3;i <= n;i++)
f[i] = f[i - 1] + f[i - 2];
但是，学了矩阵快速幂之后，我们可以利用另一种方法来推导它：
所以我们可以利用矩阵快速幂先得到状态转换矩阵的n次方，再将其于初始矩阵相乘
#include<iostream>
#include<cstring>
using namespace std;
typedef long long ll;
const int N = 2, MOD = 1e9 + 7;
struct Matrix
{
ll a[2][2];
Matrix()
{
memset(a, 0, sizeof a);
}
void init()
{
for (int i = 0; i < 2; i++)
a[i][i] = 1;
}
Matrix operator*(const Matrix& B)const
{
Matrix C;
for (int i = 0; i < 2; i++)
for (int k = 0; k < 2; k++)
for (int j = 0; j < 2; j++)
C.a[i][j] = (C.a[i][j] + a[i][k] * B.a[k][j]) % MOD;
return C;
}
Matrix operator^(const int& t)const
{
Matrix A = (*this), res;
res.init();
int p = t;
while (p)
{
if (p & 1)res = res * A;
A = A * A;
p >>= 1;
}
return res;
}
};
Matrix bace;
int main()
{
bace.a[0][0] = 0, bace.a[0][1] = 1;
bace.a[1][0] = 1, bace.a[1][1] = 1;
int n;
cin >> n;
if (n == 1 || n == 2)puts("1");
else
{
Matrix ans = bace ^ n;
cout << ans.a[0][1] << endl;//f[n] = ans.a[0][0] * f[0] + ans.a[0][1] * f[1],f[0] = 0,f[1] = 1
}
return 0;
}
5.2.2求Sn
我们先来看看式子的推导过程：
#include<iostream>
#include<cstring>
using namespace std;
typedef long long ll;
const int N = 2, MOD = 1e9 + 7;
struct Matrix
{
ll a[2][2];
Matrix()
{
memset(a, 0, sizeof a);
}
void init()
{
for (int i = 0; i < 2; i++)
a[i][i] = 1;
}
Matrix operator*(const Matrix& B)const
{
Matrix C;
for (int i = 0; i < 2; i++)
for (int k = 0; k < 2; k++)
for (int j = 0; j < 2; j++)
C.a[i][j] = (C.a[i][j] + a[i][k] * B.a[k][j]) % MOD;
return C;
}
Matrix operator^(const int& t)const
{
Matrix A = (*this), res;
res.init();
int p = t;
while (p)
{
if (p & 1)res = res * A;
A = A * A;
p >>= 1;
}
return res;
}
};
Matrix bace;
int main()
{
int a, b, n, m;
cin >> a >> b >> n >> m;
bace.a[0][0] = a, bace.a[0][1] = b;
bace.a[1][0] = 1, bace.a[1][1] = a;
Matrix t = bace ^ n;
int Xn = t.a[0][0];
cout << (Xn * 2) % MOD << endl;
return 0;
}
5.2.3梦想之弧
先看看推导过程：
#include<iostream>
#include<cstring>
using namespace std;
typedef long long ll;
const int MOD = 1e9 + 7;
int A0, AX, AY, B0, BX, BY;
struct Matrix
{
int a[6][6];
Matrix()
{
memset(a, 0, sizeof a);
}
void init()
{
for (int i = 0; i < 5; i++)
a[i][i] = 1;
}
Matrix operator*(const Matrix& B)const
{
Matrix C;
for (int i = 0; i < 5; i++)
for (int k = 0; k < 5; k++)
for (int j = 0; j < 5; j++)
C.a[i][j] = (C.a[i][j] + a[i][k] * B.a[k][j]) % MOD;
return C;
}
Matrix operator^(const int& t)const
{
Matrix A = (*this), res;
res.init();
int p = t;
while (p)
{
if (p & 1)res = res * A;
A = A * A;
p >>= 1;
}
return res;
}
};
Matrix bace, once;
int main()
{
int n;
cin >> n;
cin >> A0 >> AX >> AY;
cin >> B0 >> BX >> BY;
bace.a[0][0] = AX * BX, bace.a[0][4] = 1;
bace.a[1][0] = AX * BY, bace.a[1][1] = AX;
bace.a[2][0] = AY * BX, bace.a[2][2] = BX;
bace.a[3][0] = AY * BY, bace.a[3][1] = AY, bace.a[3][2] = BY, bace.a[3][3] = 1;
bace.a[4][4] = 1;
int a1 = A0 * AX + AY, b1 = B0 * BX + BY;
once.a[0][0] = a1 * b1, once.a[0][1] = a1, once.a[0][2] = b1, once.a[0][3] = 1, once.a[0][4] = A0 * B0;
Matrix ans = once * (bace ^ (n - 1));
cout << ans.a[0][4] << endl;
return 0;
}
5.3矩阵乘法的优化
在上边我们用到的矩阵乘法是三重循环的，时间复杂度为O(n3)，在进行较大的矩阵运算时便会显得吃力。下边我们介绍两种方法来对矩阵的乘法进行优化。
5.3.1省去为0的运算
Matrix operator*(const Matrix& B)const
{
Matrix C;
for(int i = 0;i < n;i++)
{
for(int k = 0;k < n;k++)
{
if(!a[i][k])continue;//为0时乘积为0，直接跳过
for(int j = 0;j < n;j++)
{
C.a[i][j] = (C.a[i][j] + a[i][k] * B.a[k][j]) % MOD;
}
}
}
}
5.3.2利用二进制优化
这个方法需要用到模板类bitset，我们先来学习一下待会要用到的语法。
- bitset是将数据用二进制储存下来并且进行编辑。在定义的时候，我们需要规定其长度，例如bitset<10>a，就是创建一个十位的二进制空间a。
bitset<10>a;//创建，a的初始状态为：0000000000
a.set(3);//将第3位变成1，此时a的状态为：0000000100
a.reset();//将所有位清零，此时a的状态为：0000000000
a.set(4),a.set(5),a.set(2);
a.count();//统计a中的1的数量，此时a的状态为：0000011010，a.count() = 3
在掌握以上语法之后，我们就可以来看看它的优化方法了：
首先我们创建两个数组：
bitset<1000>bt[805][3],ct[805][3];
然后，我们来看看它的运算过程：
如此如此便可得到结果，我们来看看最终的代码，这时候的时间复杂度为O(n2)：
#include<iostream>
#include<bitset>
using namespace std;
const int N = 805;
bitset<1000>bt[N][3], ct[N][3];
int main()
{
int n, c;
cin >> n;
for (int i = 1; i <= n; i++)
for (int j = 0; j < 3; j++)
bt[i][j].reset(), ct[i][j].reset();
for(int i = 1;i <= n;i++)
for (int j = 1; j <= n; j++)
{
cin >> c;
bt[i][c % 3].set(j);
}
for(int i = 1;i <= n;i++)
for (int j = 1; j <= n; i++)
{
cin >> c;
ct[j][c % 3].set(i);
}
for (int i = 1; i <= n; i++)
{
for (int j = 1; j <= n; j++)
{
int c = ((bt[i][1] & ct[j][1]).count()
+ 2 * (bt[i][1] & ct[j][2]).count()
+ 2 * (bt[i][2] & ct[j][1]).count()
+ 2 * (bt[i][2] & ct[j][2]).count()) % 3;
cout << c << ' ';
}
cout << endl;
}
return 0;
}
5.4矩阵快速幂加速DP
5.4.1Clarke and digits
这道题可是重量级。它运用的是矩阵快速幂和数位DP。
首先 ，我们要考虑的是如何表示状态，根据题干，我们可以提取三个信息：长度、模数、数位。
我们创建一个三位数组，每一维度分别表示一条信息。如==a[i] [j] [k]==表示的长度为i、模数为j且尾数为k的数的数量。
下面来看看它的表示原理和状态更新的方式：
但是矩阵应该由二维数组表示，那我们开看看，能不能通过某种方式来减去一个维度。答案当然是可以的，先来看看压缩之后的结果：
a[i] [j] [k] = a[i] [j * 10 + k]
原理何在？我们知道，j是对7取模之后的结果，所以它的范围为[0,7]；k是尾数，尾数即最后一位数，所以它的范围为[0,9]。所以数组a的两维可以表示为一个7 * 10的矩阵，而上图的压缩方法是将其变成了一个1 * 70的矩阵，0 ~ 9表示原来的第一行的状态，10 ~ 19表示原来的第二行的状态，以此类推。那么，我们就得到一个二维数组，即我们需要的矩阵了。
学到这里，结果就明了了，构建完状态矩阵和初始矩阵之后，其他的套模板即可：
状态矩阵：
void init(int k)
{
memset(bace.a,0,sizeof bace.a);
for(int i = 0;i < 70;i++)
{
int x = i / 10,y = i % 10;//余数和尾数
for(int j = 0;j < 70;j++)
{
int a = j / 10,b = j % 10;//上一状态的余数和尾数
if((a * 10 + y) == x && y + b != k)
bace.a[j][i] = 1;
}
}
}
初始矩阵：
for(int i = 1;i < 10;i++)once.a[0][(i % 7) * 10 + i] = 1;
完整代码如下：
#include<iostream>
#include<cstring>
using namespace std;
typedef long long ll;
const int N = 71, MOD = 1e9 + 7;
struct Matrix
{
int a[N][N];
Matrix()
{
memset(a, 0, sizeof a);
}
void init()
{
for (int i = 0; i < N; i++)
for (int j = 0; j < N; j++)
a[i][j] = (i == j);
}
Matrix operator*(const Matrix& B)const
{
Matrix C;
for (int i = 0; i < N; i++)
for (int k = 0; k < N; k++)
for (int j = 0; j < N; j++)
C.a[i][j] = (C.a[i][j] + a[i][k] * B.a[k][j]) % MOD;
return C;
}
Matrix operator^(const int& t)const
{
Matrix A = (*this), res;
res.init();
int p = t;
while (p)
{
if (p & 1)res = res * A;
A = A * A;
p >>= 1;
}
return res;
}
};
Matrix bace, once;
void init(int k)
{
memset(bace.a, 0, sizeof bace.a);
for (int i = 0; i < 70; i++)
{
int x = i / 10, y = i % 10;
for (int j = 0; j < 70; j++)
{
int a = j / 10, b = j % 10;
if ((a * 10 + y) % 7 == x && y + b != k)
bace.a[j][i] = 1;
}
}
for (int i = 0; i < 10; i++)bace.a[i][70] = 1;
bace.a[70][70] = 1;
}
int main()
{
memset(once.a, 0, sizeof once.a);
for (int i = 1; i < 10; i++)once.a[0][(i % 7) * 10 + i] = 1;
int t;
cin >> t;
while (t--)
{
int l, r, k;
cin >> l >> r >> k;
init(k);
Matrix ans1 = once * (bace ^ (l - 1));
Matrix ans2 = once * (bace ^ r);
ll ans = ans2.a[0][70] - ans1.a[0][70];
ans = (ans % MOD + MOD) % MOD;
cout << ans << endl;
}
return 0;
}
5.4.2我也不知道这题叫啥，因为找不到
这里没有什么编程上的语法知识，我们直接看看推导过程即可：
以下是费马小定理的衍生：
我们在求Kn时便可对p-1取模，防止数据过大
代码如下：
#include<iostream>
#include<cstring>
using namespace std;
typedef long long ll;
const int N = 3;
int n, a, b, c, p;
struct Matrix
{
int a[N][N];
Matrix()
{
memset(a, 0, sizeof a);
}
void init()
{
for (int i = 0; i < N; i++)
a[i][i] = 1;
}
Matrix operator*(const Matrix& B)const
{
Matrix C;
for (int i = 0; i < N; i++)
for (int k = 0; k < N; k++)
for (int j = 0; j < N; j++)
C.a[i][j] = (C.a[i][j] + a[i][k] * B.a[k][j]) % p;
return C;
}
Matrix operator^(const int& t)const
{
Matrix A = (*this), res;
int p = t;
res.init();
while (p)
{
if (p & 1)res = res * A;
A = A * A;
p >>= 1;
}
return res;
}
};
ll pow(int a, int b)
{
int ans = 1;
while (b)
{
if (b & 1)ans = (ans * a) % p;
a = (a * a) % p;
b >>= 1;
}
return 1;
}
Matrix bace, once;
void init()
{
memset(bace.a, 0, sizeof bace.a);
memset(once.a, 0, sizeof once.a);
bace.a[0][0] = c, bace.a[0][1] = 1, bace.a[0][2] = 1;
bace.a[1][0] = 1, bace.a[2][2] = 1;
once.a[0][0] = b, once.a[1][0] = 0, once.a[2][0] = b;
}
int main()
{
cin >> n >> a >> b >> c >> p;
if (n == 1)puts("1");
else if (n == 2)cout << pow(a, b) << endl;
else
{
p--;
Matrix ans = bace ^ (n - 2);
p++;
ans = once * ans;
cout << ans.a[0][0] << endl;
}
return 0;
}
五、拓展欧几里得算法
1.拓展欧几里得算法
裴蜀定理：对于任意整数a,b，它们的最大公约数d，关于未知整数x,y的线性不定方程成为裴蜀定理。即对于任意的整数a,b，都存在整数x,y使得ax + by = gcd(a,b)。
由快速幂可知，ax + by = by + (a - [ ⌊ a b ⌋ \lfloor\frac a b\rfloor ⌊ba⌋] * b)x
将式子提取之后，我么可以得到：ax + by = ax + b(y - [ ⌊ a b ⌋ \lfloor\frac a b\rfloor ⌊ba⌋]x)
由此可得：x = x,y = y - [ ⌊ a b ⌋ \lfloor\frac a b\rfloor ⌊ba⌋]x
int gcd(int a, int b, int& x, int& y)
{
if(!b)
{
x= 1,y = 0;
return a;
}
int d = gcd(b, a % b, y, x);
y -= a / b * x;
return d;
}
2.线性同余方程
给定a,b,m，求一个x使得$ a * x \equiv b $ (mod m)成立。
对于这个问题，我们可以将其先进行一定的变形： a ∗ x ≡ b + m ∗ y = a ∗ x − m ∗ ′ y ≡ b a * x \equiv b + m * y = a * x - m * 'y \equiv b a∗x≡b+m∗y=a∗x−m∗′y≡b
如此，该式子就被转化成了可以使用拓展欧几里得算法求结果的形式了。
#include<iostream>
using namespace std;
typedef long long ll;
int gcd(int a, int b, int& x, int& y)
{
if(!b)
{
x= 1,y = 0;
return a;
}
int d = gcd(b, a % b, y, x);
y -= a / b * x;
return d;
}
int main()
{
int a,b,m;
cin >> a >> b >> m;
int x,y;
int d = gcd(a,m,x,y);
if(b % d)puts("impossible");//如果b不是(a,m)的最大公约数的倍数，则说明无解，
else printf("%d\n",(ll)b / d * x % m);
return 0;
}
六、中国剩余定理
存在下边一元线性同余方程：
{ x m o d a 1 = m 1 x m o d a 2 = m 2 … … x m o d a k = m k \begin{cases}\\x mod ~ a1 = m1 \\ x ~ mod ~ a2 = m2 \\ …… \\ x ~ mod ~ ak = mk \end{cases} ⎩⎪⎪⎪⎨⎪⎪⎪⎧xmod a1=m1x mod a2=m2……x mod ak=mk
中国剩余定理给出了，已知a1-k和m1-k，求x的值求解方法。
首先我们来看看它的推导过程：
#include<iostream>
using namespace std;
typedef long long ll;
ll gcd(ll a,ll b,ll& x,ll& y)//求k1、k2和a1a2的最大公约数
{
if(!b)
{
x = 1,y = 0;
return a;
}
ll d = gcd(b,a % b,y,x);
y -= a / b * x;
return d;
}
int main()
{
int n;
cin >> n;
bool flag = 1;
ll a1,m1;
cin >> a1 >> m1;
for(int i = 1;i < n;i++)
{
ll a2,m2;
cin >> a2 >> m2;
ll k1,k2;
ll d = gcd(a1,a2,k1,k2);
if((m2 - m1) % d)
{
flag = 0;
break;
}
k1 *= (m2 - m1) / d;
ll t = a2 / d;
k1 = (k1 % t + t) % t;
m1 += k1 * a1;
a1 = abs(a1 / d * a2);
}
if(flag)cout << (m1 % a1 + a1) % a1 << endl;
else puts("-1");
return 0;
}
七、高斯消元
1.高斯消元解线性方程组
首先我们将等号左边的系数和右边的常数输入到二维数组中存起来，例如：
| 1 | 2 | -1 | -6 |
|---|---|---|---|
| 2 | 1 | -3 | -9 |
| -1 | -1 | 2 | 7 |
左边三列存的是等式左边各个a的值，最后一列储存的是各个b的值。
然后我们要做的就是通过矩阵的初等行列变换将其变成倒三角的形式：
- 枚举每一列，绝对值最大的那一列移到第一行
- 将第一行第一列最左边的数变为1，其他行的变为0，重复该过程，最终可以得到下列矩阵：
| 1 | 1/2 | -3/2 | -9/2 |
|---|---|---|---|
| 0 | 1 | 1/3 | -1 |
| 0 | 0 | 2/3 | 2 |
然后我们将其从下列往上处理，得到以下结果：
| 1 | 0 | 0 | 1 |
|---|---|---|---|
| 0 | 1 | 0 | -2 |
| 0 | 0 | 1 | 3 |
然后我们就可以得到结果，x1 = 1,x2 = -2,x3 = 3。
如果在出现了0=0的情况，那么说明具有无数解，如果出现了0=!0的情况，则说明无解。
const double eps = 1e-6;
int gauzz()
{
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]);
for(int i = n;i >= c;--i)a[r][i] /= a[r][c];
for(int i = r + 1;i < n;i++)
if(fabs(a[i][c]) > eps)
for(int j = c;j <= n;++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;--i)
for(int j = i + 1;j < n;j++)
a[i][n] -= a[i][j] * a[j][n];
return 1;//表示存在唯一解
}
2.高斯消元解异或线性方程
将操作换为异或之后，我们就可以得到高斯消元解异或方程组了：
#include<iostream>
using namespace std;
const int N = 110;
int n;
int a[N][N];
int gauss()
{
int r, c;
for (r = c = 0; c < n; c++)
{
int t = r;
for(int i = r;i < n;i++)
if (a[i][c])
{
t = i;
break;
}
if (!a[t][c])continue;
for (int i = 0; i < n; i++)swap(a[r][i], a[t][i]);
for (int i = r + 1; i < n; i++)
if (a[i][c])
for (int j = 0; j <= n; j++)
a[i][j] ^= a[r][j];
++r;
}
if (r < n)
{
for (int i = r; i < n; i++)
if (a[i][n])return 2;
return 1;
}
for (int i = 0; i < n; i++)
for (int j = i + 1; j < n; j++)
a[i][n] ^= a[i][j] & a[j][n];
return 0;
}
int main()
{
cin >> n;
for (int i = 0; i < n; i++)
for (int j = 0; j <= n; j++)
cin >> a[i][j];
int t = gauss();
if (t == 0)
{
for (int i = 0; i < n; i++)
cout << a[i][n] << endl;
}
else if (t == 1)puts("Multiple sets of solutions");
else puts("No solution");
return 0;
}
八、求组合数
1.int范围内数字较小时组合数的求法
组合数是指从M个物品中选出N个物品的方法个数，存在 C M N = C M − 1 N + C M − 1 N − 1 C_{M}^{N} = C_{M - 1}^{N} + C_{M - 1}^{N-1} CMN=CM−1N+CM−1N−1，如此我们可以利用dp思想得到求N内所有数的组合数：
void init()
{
for(int i = 0;i < n;i++)
for(int j = 0;j <= i;j++)
if(!f[i][j])f[i][j] = 1;
else f[i][j] = f[i - 1][j] + f[i - 1][j - 1];
}
2.int范围内数字较大时组合数的求法
由组合数的定义来看， C M N = M ! ( M − N ) ! ∗ N ! C_{M}^{N} = \frac{M!}{(M - N)! * N!} CMN=(M−N)!∗N!M!
我们可以利用逆元的知识，即 a b \frac{a}{b} ba ≡ a * x(mod m)，我们将除的运算换成乘逆元，然后遍历处理出范围内所有数的阶乘和阶乘的逆元即可。
int qpow(int a,int k,int p)
{
int ans = 1;
while(k)
{
if(k & 1)ans = (ll)ans * a % p;
a = (ll) a * a % p;
k >>= 1;
}
return ans;
}
void init()
{
fct[0] = fct[0] = 1;
for(int i = 1;i < N;i++)
{
fct[i] = (ll)fct[i - 1] * i % mod;
infct[i] = (ll)infct[i - 1] * qpow(i,mod - 2,mod) % mod;
}
}
然后在输出的时候：
printf("%d\n",(ll)fct[a] * infct[b] % mod * infct[a - b] % mod);
3.long long范围内组合数的求法
在数字范围需要用long long来表示的时候，如果使用遍历的方法来做的话，时间复杂度和市局范围会应为过大而变得麻烦。所以此时我们应该使用卢卡斯定理来求组合数。
卢卡斯定理： C M N C_{M}^{N} CMN = ∏ i = 1 k x i \prod_{i = 1}^k x_i ∏i=1kxi(mod p)
在程序中，我们可以使用下面的形式：
C M N C_{M}^{N} CMN = C M m o d p N m o d p C_{M~mod~p}^{N~mod~p} CM mod pN mod p * C M / p N / p C_{M / p}^{N / p} CM/pN/p(mod p)
int qpow(int a,int k)
{
int ans = 1;
while(k)
{
if(k & 1)ans = (ll)ans * a % mod;
a = (ll)a * a % mod;
k >>= 1;
}
return ans;
}
int C(int a,int b)
{
if(b > a)return 0;
int ans = 1;
for(int i = 1,j = a;i <= b;i++,j--)
{
ans = (ll)ans * j % mod;
ans = (ll)ans * qpow(i, mod - 2) % mod;
}
return ans;
}
int lucass(ll a,ll b)
{
if(a < mod && b < mod)return C(a,b);
else return (ll)C(a % mod,b % mod) * lucass(a / mod,b / mod) % mod;
}
4.高精度与组合数
从定义上来看，$\C_{M}^{N} $ = M ! ( M − N ) ! ∗ N ! \frac{M!}{(M - N)! * N!} (M−N)!∗N!M!，我们将三个阶乘数进行质因数分解，然后将质因子相消，最后将剩下的质因子进行高精度乘法即可。
#include<iostream>
#include<vector>
using namespace std;
const int N = 5010;
typedef vector<int> VI;
int primes[N],idx;
bool vis[N];
int num[N];
void get_prime(int n)
{
for(int i = 2;i <= n;i++)
{
if(!vis[i])primes[idx++] = i;
for(int j = 0;primes[j] * i <= n;j++)
{
vis[primes[j] * i] = 1;
if(i % primes[j] == 0)break;
}
}
}
int get(int num, int p)//质因子个数
{
int c = 0;
while(num)
{
c += num / p;
num /= p;
}
return c;
}
VI mul(VI a,int b)
{
VI 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;
}
return c;
}
int main()
{
int a,b;
cin >> a >> b;
get_prime(a);
for(int i = 0; i < idx; i++)
{
int p = primes[i];
num[i] = get(a, p) - get(b, p) - get(a - b, p);
}
VI ans;
ans.push_back(1);
for(int i = 0;i < idx;i++)
for(int j = 0;j < num[i];j++)
ans = mul(ans,primes[i]);
for(int i = ans.size() - 1;~i;--i)
cout << ans[i];
cout << endl;
return 0;
}
5.卡特兰数
卡特兰数的定义： C 2 n n − C 2 n n − 1 = C 2 n n n + 1 C_{2n}^{n} - C_{2n}^{n-1} = \frac{C_{2n}^{n}}{n+1} C2nn−C2nn−1=n+1C2nn
证明如下下：
代码实现如下：
#include<iostream>
using namespace std;
const int mod = 1e9 + 7;
typedef long long ll;
int qpow(int a, int k)
{
int res = 1;
while (k)
{
if (k & 1)res = (ll)res * a % mod;
a = (ll)a * a % mod;
k >>= 1;
}
return res;
}
int C(int a, int b)
{
int res = 1;
for (int i = 1, j = a; i <= b; i++, j--)
{
res = (ll)res * j % mod;
res = (ll)res * qpow(i, mod - 2) % mod;
}
return res;
}
int main()
{
int n;
cin >> n;
cout << (ll)C(2 * n, n) % mod * qpow(n + 1,mod - 2) % mod << endl;
return 0;
}
九、容斥原理
以上为容斥原理的原理图，将所有数相加，再减去所有的重复元素，让集合中的数字不重不漏。
看例题：https://www.acwing.com/problem/content/892/
我们将集合用不同的组合方式相乘，记下其中小于n的数，得到的数和它的倍数的数量就是这题的答案。
因为一共由m个p，m < 16，所以我们利用1 ~ 2m的二进制来表示所有的排列组合方式。
#include<iostream>
using namespace std;
typedef long long ll;
int n,m;
int p[20];
int main()
{
cin >> n >> m;
for(int i = 0;i < m;i++)cin >> p[i];
int ans = 0;
for(int i = 1;i < 1 << m;i++)
{
int t = 1,s = 0;
for(int j = 0;j < m;j++)
{
if(i >> j & 1)
{
if((ll)t * p[j] > n)
{
t = -1;
break;
}
t *= p[j];
s++;
}
}
if(t != -1)
{
if(s % 2)ans += n / t;
else ans -= n / t;
}
}
cout << ans << endl;
return 0;
}
十、博弈论
什么是博弈论？
博弈论是二人在平等的对局中各自利用对方的策略变换自己的对抗策略，达到取胜的目的。
博弈论不存在模板，类似于贪心算法，如何得到思路才是关键。下边我们看看一道例题来帮助理解博弈论。
例题：Nim游戏
给定 n 堆石子，两位玩家轮流操作，每次操作可以从任意一堆石子中拿走任意数量的石子（可以拿完，但不能不拿），最后无法进行操作的人视为失败。
问如果两人都采用最优策略，先手是否必胜。
Input
第一行包含整数 n。
第二行包含 n 个数字，其中第 i 个数字表示第 i 堆石子的数量。
Output
如果先手方必胜，则输出
Yes。否则，输出
No。
当轮到某位玩家时场上已经没有石子的时候，则该玩家失败，此时各堆石子的数量异或起来为0，即0 ^ 0 ^ 0 ^ 0 ^…… ^ 0 = 0。但是当所有堆的石子数异或起来不为0的时候，则其一定不为0。而且存在以下两个定理，当异或值为0时，从任意堆中拿走任意数，剩下的数异或值一定不为零；当异或值不为0时，从任意堆中拿走特定数（一定存在），剩下的数异或值一定为0。
证明1：
a1 ^ a2 ^ a3 ^ …… ^ ai ^ …… ^ an = 0，从ai中拿走任意数量的数使之变为ai’，假设其异或值为i，则a1 ^ a2 ^ a3 ^ …… ^ ai’ ^ …… ^ an = 0，将两式异或，由于除了i位其他位上均相同，所以就有ai ^ ai’ = 0，则ai = ai’，与假设矛盾，所以a1 ^ a2 ^ a3 ^ …… ^ ai’ ^ …… ^ an = !0。
证明2：
a1 ^ a2 ^ a3 ^ …… ^ ai ^ …… ^ an = !0，存在值x < ai，拿走ai - ai ^ x使第i位变为ai ^ x，则!0 ^ x = 0。
综上，当两人均不会出现失误的时候，初始数的异或值便已经决定了游戏的胜负，当异或值不为0的时候，先手胜，否则后手胜。
#include<iostream>
using namespace std;
int main()
{
int n;
cin >> n;
int res = 0;
while(n--)
{
int num;
cin >> num;
res ^= num;
}
if(res)puts("Yes");
else puts("No");
return 0;
}
十、博弈论
什么是博弈论？
博弈论是二人在平等的对局中各自利用对方的策略变换自己的对抗策略，达到取胜的目的。
博弈论不存在模板，类似于贪心算法，如何得到思路才是关键。下边我们看看一道例题来帮助理解博弈论。
例题：Nim游戏
给定 n 堆石子，两位玩家轮流操作，每次操作可以从任意一堆石子中拿走任意数量的石子（可以拿完，但不能不拿），最后无法进行操作的人视为失败。
问如果两人都采用最优策略，先手是否必胜。
Input
第一行包含整数 n。
第二行包含 n 个数字，其中第 i 个数字表示第 i 堆石子的数量。
Output
如果先手方必胜，则输出
Yes。否则，输出
No。
当轮到某位玩家时场上已经没有石子的时候，则该玩家失败，此时各堆石子的数量异或起来为0，即0 ^ 0 ^ 0 ^ 0 ^…… ^ 0 = 0。但是当所有堆的石子数异或起来不为0的时候，则其一定不为0。而且存在以下两个定理，当异或值为0时，从任意堆中拿走任意数，剩下的数异或值一定不为零；当异或值不为0时，从任意堆中拿走特定数（一定存在），剩下的数异或值一定为0。
证明1：
a1 ^ a2 ^ a3 ^ …… ^ ai ^ …… ^ an = 0，从ai中拿走任意数量的数使之变为ai’，假设其异或值为i，则a1 ^ a2 ^ a3 ^ …… ^ ai’ ^ …… ^ an = 0，将两式异或，由于除了i位其他位上均相同，所以就有ai ^ ai’ = 0，则ai = ai’，与假设矛盾，所以a1 ^ a2 ^ a3 ^ …… ^ ai’ ^ …… ^ an = !0。
证明2：
a1 ^ a2 ^ a3 ^ …… ^ ai ^ …… ^ an = !0，存在值x < ai，拿走ai - ai ^ x使第i位变为ai ^ x，则!0 ^ x = 0。
综上，当两人均不会出现失误的时候，初始数的异或值便已经决定了游戏的胜负，当异或值不为0的时候，先手胜，否则后手胜。
#include<iostream>
using namespace std;
int main()
{
int n;
cin >> n;
int res = 0;
while(n--)
{
int num;
cin >> num;
res ^= num;
}
if(res)puts("Yes");
else puts("No");
return 0;
}
8814
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
