本文深入探讨数论中的斯特林数、原根与离散对数、高斯消元法及二次剩余等核心概念，提供算法实现与实例解析。
md数论，我快吐了，花Q~花Q.....
（一）斯特林数：
第一类斯特林数：把n个对象分成m个非空循环排列的方案数：
S(n,m) = (n-1)*S(n-1, m) + S(n-1, m-1).
第二类斯特林数：把n个对象划分到k个不可区分的非空盒子的划分数：
S(n, m) = m*S(n-1, m) + S(n-1, m-1).
把n个对象划分到k个不可区分的可以空的盒子的划分数：
S(n, 1) + S(n, 2) + .... + S(n, k).
（二）原根与离散对数：
对于两个互质的正整数a, m由欧拉定理可知，存在正整数d <= m-1, 使得a^d ≡ 1(mod m), 定义Ordm(a)为最小的d
由定义可知一定有Ordm(a) <= φ(m), 如果有Ordm(a) = φ(m), 则称a为模m的一个原根, 具有以下性质：
具有原根的数字只有以下几种形式:2, 4, 2p^n, p^n, 其中p为奇素数, n为正整数
m的最小原根的大小是O(m^1/4)
若g为m的原根，则g^d是m的原根的充要条件是(d, φ(m)) = 1, 且m的原根的个数等于φ(φ(m)).
给定n ,a ,b 求解 a^x ≡ b(mod n), 称为离散对数问题，这里只讨论n为素数的情况，下面给出一种O(√nlogn)算法：
首先，算出 x = 0, 1, 2, ..., m时a^x(mod n)的值。由对于任何一个大于m的数字k，存在唯一的一对p和q, q < m,
使得 k = pm+q, 假设a^k ≡ b(mod n), 可以改写为 a^q ≡ b*(a^m)^p (mod n). 因此，只需枚举p的大小，依次
计算出b*(a^-m)^p (mod n)的值, 然后在x = 0,1,2, ... ,m时a^m(mod n)的值中进行查找即可, 在m = √n时，算法较好。
经过多方对比，从网上找到的一个既能当做普通BSGS使用，又能作为扩展BSGS使用的板子：
/*计算a^x ≡ b(mod c), 无解返回-1*/
typedef long long ll;
const int N = 65536;
struct Node{
ll idx, val;
}baby[N];
bool cmp(Node x, Node y){
if(x.val != y.val) return x.val < y.val;
else return x.idx < y.idx;
}
ll gcd(ll x, ll y){
if(y == 0) return x;
else return gcd(y, x%y);
}
ll exgcd(ll a, ll b, ll &x, ll &y){
if(b == 0){
x = 1;
y = 0;
return a;
}
ll d = exgcd(b, a%b, x, y);
ll t = x;
x = y;
y = t - a/b * y;
return d;
}
ll inval(ll a, ll b, ll n)
{
ll e, x, y;
exgcd(a, n, x, y);
e = x*b % n;
return e < 0 ? e+n : e;
}
ll powmod(ll a, ll b, ll mod)
{
ll ret = 1%mod, t = a%mod;
while(b){
if(b & 1) ret = ret*t % mod;
t = t*t % mod;
b >>= 1;
}
return ret;
}
ll BinSearch(ll num, ll m)
{
ll low = 0, high = m-1, mid;
while(low <= high){
mid = (low+high) >> 1;
if(baby[mid].val == num) return baby[mid].idx;
if(baby[mid].val < num) low = mid+1;
else high = mid-1;
}
return -1;
}
ll BSGS(ll A, ll B, ll C)
{
ll tmp, D = 1%C, temp;
for(ll i = 0, tmp = 1%C; i < 100; i++, tmp = tmp*A%C){
if(tmp == B) return i;
}
ll d = 0;
while((temp = gcd(A, C)) != 1){
if(B%temp) return -1;
d++;
C /= temp;
B /= temp;
D = (A / temp) * D % C;
}
ll m = (ll)ceil(sqrt( (double)C ) );
for(ll i = 0, tmp = 1%C; i <= m; i++, tmp = tmp*A%C){
baby[i].idx = i;
baby[i].val = tmp;
}
sort(baby, baby+m+1, cmp);
ll cnt = 1;
for(ll i = 1; i <= m; i++){
if(baby[i].val != baby[cnt-1].val) baby[cnt++] = baby[i];
}
ll am = powmod(A, m, C);
for(ll i = 0; i <= m; i++, D = D*am%C){
ll tmp = inval(D, B, C);
if(tmp >= 0){
ll pos = BinSearch(tmp, cnt);
if(pos != -1) return i*m+pos+d;
}
}
return -1;
}
（三）高斯消元法：
这个我还是看这位大佬的博客就好了：http://www.cppblog.com/menjitianya/archive/2014/06/08/207226.html
(从网上找的模板怎么感觉参差不齐啊，代码大多比较乱，这个的注释还不错)
在贴一个模板：
//高斯消元法
#include<cstdio>
#include<iostream>
#include<cstring>
#include<cmath>
#include<algorithm>
using namespace std;
const int maxn = 50;
int a[maxn][maxn]; //增广矩阵
int x[maxn]; //解集
bool free_x[maxn]; //用标记是否是不确定的变元
inline int gcd(int a, int b){
if(b == 0) return a;
else return gcd(b, a%b);
}
inline int lcm(int a, int b){
return a/gcd(a, b)*b;
}
/*
-2表示有浮点数解，无整数解
-1表示无解，0表示有唯一解
大于0表示无穷解，并返回自由变元的个数
有equ个方程，var个变元，增广矩阵行数为equ，从0->equ-1, 列数为var+1, 从0->var;
*/
int Gauss(int equ, int var)
{
int i, j, k;
int max_r; //当前这列的绝对值最大的行
int col; //当前处理这列
int ta, tb, LCM, temp;
int free_x_num, free_index;
for(int i = 0; i <= var; i++){
x[i] = 0;
free_x[i] = true;
}
//转化为阶矩阵
col = 0; //当前处理的列
for(k = 0; k < equ && col < var; k++, col++){
//找到该col列绝对值最大的那行并与第k行交换(为了在除法时减小误差)
max_r = k;
for(i = k+1; i < equ; i++){
if(abs(a[i][col]) > abs(a[max_r][col]))
max_r = i;
}
if(max_r != k){ //与第k行交换
for(j = k; j < var+1; j++)
swap(a[k][j], a[max_r][j]);
}
if(a[k][col] == 0){
//说明该col列第k行以下全为0了,则处理当前行的下一列
k--;
continue;
}
for(i = k+1; i < equ; i++){
//枚举要删去的行
if(a[i][col] != 0){
LCM = lcm(abs(a[i][col]), abs(a[k][col]));
ta = LCM / abs(a[i][col]);
tb = LCM / abs(a[k][col]);
if(a[i][col] * a[k][col] < 0) tb = -tb;
for(j = col; j < var+1; j++){
a[i][j] = a[i][j]*ta - a[k][j]*tb;
}
}
}
}
// 1.无解的情况: 化简的增广阵中存在(0, 0, ..., a)这样的行(a != 0).
for(i = k; i < equ; i++){
if(a[i][col] != 0) return -1;
}
/* 2.无穷解的情况: 在var * (var + 1)的增广阵中出现(0, 0, ..., 0)这样的行
且出现的行数即为自由变元的个数.*/
if (k < var){
//自由变元有var-k个
for(i = k-1; i >= 0; i--){
free_x_num = 0;
for(j = 0; j < var; j++){
if(a[i][j] != 0 && free_x[j]){
free_x_num++;
free_index = j;
}
}
if(free_x_num > 1) continue; //无法求解出确定的变元
temp = a[i][var];
for(j = 0; j < var; j++){
if(a[i][j] != 0 && j != free_index) temp -= a[i][j]*x[j];
}
x[free_index] = temp / a[i][free_index]; //求出该变元
free_x[free_index] = 0; //该变元是确定的
}
return var-k;
}
/* 3.唯一解的情况: 在var * (var + 1)的增广阵中形成严格的上三角阵.
计算出Xn-1, Xn-2 ... X0.*/
for(i = var-1; i >= 0; i--){
temp = a[i][var];
for(j = i+1; j < var; j++){
if(a[i][j] != 0) temp -= a[i][j]*x[j];
}
if(temp % a[i][i] != 0) return -2;
x[i] = temp / a[i][i];
}
return 0;
}
int main()
{
int i, j;
int equ, var;
while(~scanf("%d %d", &equ, &var)){
//清存，输入
memset(a, 0, sizeof(0));
for(i = 0; i < equ; i++){
for(j = 0; j <= equ; j++){
scanf("%d", &a[i][j]);
}
}
int free_num = Gauss(equ, var);
if(free_num == -1) cout << "无解" << endl;
else if(free_num == -2) cout << "有浮点数解，无整数解" << endl;
else if(free_num > 0){
cout << "自由变元个数为: " << free_num << endl;
for(i = 0; i < var; i++){
if(free_x[i]) printf("x%d不确定\n", i+1);
else printf("x%d = %d\n", i+1, x[i]);
}
}
else{
for(i = 0; i < var; i++){
printf("x%d = %d\n", i+1, x[i]);
}
}
cout << endl;
}
return 0;
}求自由变元那一部分貌似很多其他模板是直接返回var-k处理的.....
（四）二次剩余：
有一位真*大佬，搞出了Cipolla算法，来学习膜拜一下。
博客推荐：https://blog.youkuaiyun.com/doyouseeman/article/details/52033204
扫一扫
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
