萌新学数论啦（二）（未完待续）

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.csdn.net/doyouseeman/article/details/52033204