Algorithm Notes 3.2 数论

BSGS

  • a x ≡ b (   m o d   p ) a^x \equiv b (\bmod p) axb(modp) 的最小自然数 x x x

  • a ⊥ p a \perp p ap

    • 由抽屉原理,只需找到在 [ 0 , p ) [0,p) [0,p) 内的解,设 t = ⌈ p ⌉ t = \lceil \sqrt p \rceil t=p ,则 [ 0 , p ) [0,p) [0,p) 内的整数均可以用 { x = i t − j , i ∈ [ 1 , t ] , j ∈ ( 0 , t ] } \{x = it - j,i\in [1,t], j \in(0,t]\} {x=itj,i[1,t],j(0,t]}
    • 因此可将原方程化为 a i t ≡ b a j (   m o d   p ) a^{it} \equiv ba^j(\bmod p) aitbaj(modp),将右侧插入哈希表,左侧暴力枚举即可。
  • a ⊥̸ p a \not \perp p ap,不断取 d n = ( a , p ∏ j = 1 n − 1 d j ) d_n = (a,\dfrac{p}{\prod\limits_{j = 1}^{n - 1}d_j}) dn=(a,j=1n1djp),直至 ( a , p ∏ j = 1 n d j ) = 1 (a,\dfrac{p}{\prod\limits_{j = 1}^{n}d_j}) = 1 (a,j=1ndjp)=1,设 D = ∏ j = 1 n d j D = \prod\limits_{j = 1}^{n}d_j D=j=1ndj,取 a ′ = a n D , b ′ = b D , p ′ = p D a' = \dfrac{a^{n}}{D},b'=\dfrac{b}{D},p'=\dfrac{p}{D} a=Dan,b=Db,p=Dp,原方程可化为:
    a x − n ≡ b ′ ( a ′ ) − 1 (   m o d   p ′ ) a^{x - n}\equiv b'(a')^{-1}(\bmod p') axnb(a)1(modp)

    • 暴力枚举 x < n x < n x<n 的情况, x ≥ n x \ge n xn 的情况用互质的做法解决。
  • a ⊥ p a\perp p ap,只调用函数 solve_BSGS 即可,否则沿用整个程序。

#include <bits/stdc++.h>

const int Maxn = 1e9;
int a, p, b; 

struct hashtable
{
    #define MOD 4999963
    #define H 5000005
    
    int adj[H], cst[H], to[H], nxt[H], stk[H];
    int top, T;
    bool vis[H];
    
    inline void Clear()
    {
        while (top)
        {
			int x = stk[top--];
			adj[x] = 0;
            vis[x] = false;
        }
        T = 0;
    }
    
    inline int Find(int y)
    {
        int x = y % MOD;
        for (int e = adj[x]; e; e = nxt[e])
            if (to[e] == y)
                return cst[e];
        return Maxn;        
    }
    
    inline void Insert(int y, int z)
    {
        int x = y % MOD;
        if (!vis[x])
            vis[stk[++top] = x] = true;
        nxt[++T] = adj[x]; adj[x] = T; to[T] = y; cst[T] = z;
    }
}ht;

inline int quick_pow(int x, int k, int mod)
{
	int res = 1;
	while (k)
	{
		if (k & 1)
			res = 1ll * res * x % mod;
		x = 1ll * x * x % mod;
		k >>= 1;
	}
	return res;
}

inline int ex_gcd(int a, int b, int &x, int &y)
{
	if (!b)
	{
		x = 1; 
		y = 0; 
		return a;
	}
	int res = ex_gcd(b, a % b, y, x); 
	y -= a / b * x;
	return res; 
}

inline int solve_equ(int a, int b, int c)
{
	int x, y; 
	int d = ex_gcd(a, b, x, y); 
	if (c % d != 0) return -1; 
	
	int mod = b / d; 
	return (1LL * c / d * x % mod + mod) % mod; 
}

inline int solve_BSGS(int a, int b, int p)
{
	int t = ceil(sqrt(p));
	ht.Clear(); 
	
	int tmp = b; 
	for (int i = 1; i <= t; ++i)
	{
		tmp = 1LL * tmp * a % p;
		ht.Insert(tmp, i); 
	}
	
	int pw = quick_pow(a, t, p); 
	tmp = pw; 
	for (int i = 1; i <= t; ++i)
	{
		int res = ht.Find(tmp);
		if (res != Maxn)
			return i * t - res;
		tmp = 1LL * tmp * pw % p; 
	}
	return -1; 
}

inline bool check()
{
	int k = 1 % p; 
	for (int i = 0; i <= 40; ++i)
	{
		if (k == b)
			return printf("%d\n", i), true; 
		k = 1LL * k * a % p; 
	}
	if (!a)
		return puts("No Solution"), true; 
	return false; 
}

int main()
{   // 求 a^x = b(mod p) 的最小自然数 x,无解输出 No Solution
	while (scanf("%d%d%d", &a, &p, &b), a || p || b)
	{
		a %= p, b %= p; 
		if (check())
			continue; 
		
		int d; 
		int ap = 1, n = 0; 
		bool flg = false; 
		
		while ((d = std::__gcd(a, p)) != 1)
		{
	  		++n; 
			ap = 1LL * ap * (a / d) % p; 
			p /= d; 
			
			if (b % d)
			{
				flg = true; 
				break; 
			}
			b /= d; 
		}
		
		if (flg)
			puts("No Solution"); 
		else
		{
			int res = solve_BSGS(a, 1LL * b * solve_equ(ap, p, 1) % p, p); 
			if (res == -1)
				puts("No Solution"); 
			else
				printf("%d\n", res + n); 
		}
	}
	return 0; 
}

质因数分解

Miller Rabin 素性测试

  • 根据费马小定理,我们能够得出一种检验素数的思路(即费马素性测试):
    • 验证 ∀ 1 < a < n , a n − 1 ≡ 1 ( m o d    n ) \forall 1 < a < n, a^{n - 1} \equiv 1(\mod n) ∀1<a<n,an11(modn)
  • 很遗憾,上述做法并不能准确地判断素数,对于满足上述条件的合数,我们称之为卡迈克尔数
  • Miller Rabin \text{Miller Rabin} Miller Rabin 素性测试是基于费马素性测试的优化。
  • 二次探测定理 P P P 为奇素数,则 x 2 ≡ 1 ⇔ x ≡ 1 ( m o d    P ) 或 x ≡ P − 1 ( m o d    P ) x^2 \equiv 1 \Leftrightarrow x\equiv 1(\mod P) 或 x \equiv P - 1(\mod P) x21x1(modP)xP1(modP)
  • n − 1 = u × 2 t n - 1 = u \times 2^t n1=u×2t,随机取一个 a a a,先求出 a u ( m o d    n ) a^u(\mod n) au(modn),对该数平方 t t t 次,若在这一过程中出现在模 n n n 意义下为 1 或为 n − 1 n - 1 n1,可认为通过了此轮测试。
  • 经过理论证明,单轮错误率大约 1 4 \frac{1}{4} 41 T T T 轮测试的错误率为 4 − T 4^{-T} 4T,一般取 T = 8 ∼ 12 T = 8 \sim 12 T=812 即可。
  • 若取 a = 2 , 3 , 5 , 7 , 11 , 13 , 17 , 19 , 23 , 29 , 31 , 37 a = 2,3,5,7,11,13,17,19,23,29,31,37 a=2,3,5,7,11,13,17,19,23,29,31,37(即前 12 个质数)可保证 n < 2 64 n < 2^{64} n<264 确定性判素。

生日悖论

  • k k k [ 1 , n ] [1, n] [1,n] 的随机整数两两不相同的概率,设该事件为 A A A,则
    P ( A ) = ∏ i = 1 k n − i + 1 n = ∏ i = 1 k ( 1 − i − 1 n ) P(A) = \prod \limits_{i = 1}^{k} \frac{n - i + 1}{n} = \prod\limits_{i = 1}^{k}(1 - \frac{i - 1}{n}) P(A)=i=1knni+1=i=1k(1ni1)
  • n n n 较大时,根据近似 1 − i − 1 n ≈ ( 1 − 1 n ) i − 1 1 - \frac{i - 1}{n} \approx (1 - \frac{1}{n})^{i - 1} 1ni1(1n1)i1,且 e ≈ ( 1 + 1 n ) n e \approx (1 + \frac{1}{n})^n e(1+n1)n,原式可化为
    P ( A ) ≈ ∏ i = 1 k ( 1 − 1 n ) i − 1 = ( 1 − 1 n ) k ( k − 1 ) 2 ≈ e − k ( k − 1 ) 2 n P(A) \approx \prod \limits_{i = 1}^{k}(1 - \frac{1}{n})^{i - 1} = (1 - \frac{1}{n})^{\frac{k(k - 1)}{2}} \approx e^{-\frac{k(k - 1)}{2n}} P(A)i=1k(1n1)i1=(1n1)2k(k1)e2nk(k1)
  • 观察该式可知,当 k k k 取到 O ( n ) \mathcal O(\sqrt n) O(n ) 级别时, P ( A ) P(A) P(A) 会发生骤降。
  • 因此 k k k [ 1 , n ] [1,n] [1,n] 的随机整数首次出现相同数字时 k k k 的期望为 O ( n ) \mathcal O(\sqrt n) O(n )

Pollard Rho 算法

  • Pollard Rho \text{Pollard Rho} Pollard Rho 算法是一种基于随机的质因数分解算法,可以以期望时间复杂度 O ( n 1 4 ) \mathcal O(n^{\frac{1}{4}}) O(n41) 找到合数 n n n 的某个非平凡因子。
  • 考虑选取适当的 k k k 使得 1 < d = gcd ⁡ ( k , n ) < n 1 < d = \gcd(k,n) < n 1<d=gcd(k,n)<n,则显然 d d d n n n 的一个约数,这样的 k k k 相对较多。
  • n n n 本身是质数,可直接用 Miller Rabin \text{Miller Rabin} Miller Rabin 素性测试判定。否则构造伪随机数列 x i + 1 = ( x i 2 + c ) m o d    n x_{i + 1} = (x_i^2 + c)\mod n xi+1=(xi2+c)modn,其中 c c c 为初始时指定的某一随机数,设 m m m n n n 的最小非平凡因子,令 y i = x i m o d    m y_i = x_i \mod m yi=ximodm,则由生日悖论,满足 ∃ i < j , y i = y j \exist i<j, y_i = y_j i<j,yi=yj 所需的序列 x x x 的期望长度为 O ( m ) ≤ O ( n 1 4 ) \mathcal O(\sqrt m) \le \mathcal O(n^{\frac{1}{4}}) O(m )O(n41),期望枚举 O ( m ) \mathcal O(\sqrt m) O(m ) i i i 我们便能得到 n n n 的一个约数 gcd ⁡ ( ∣ x i − x j ∣ , n ) \gcd(|x_i - x_j|, n) gcd(xixj,n)
  • 可实际情况并不允许我们去暴力枚举所有 i , j i,j i,j,考虑到 x i x_i xi 的周期性,我们采用一种基于倍增的实现方式,即正序枚举 t ≥ 0 t\ge 0 t0,每次检查 gcd ⁡ ( ∣ x i − x i + k ∣ , n ) ( 1 ≤ k ≤ 2 t ) \gcd(|x_i - x_{i+k}|,n)(1 \le k \le 2^t) gcd(xixi+k,n)(1k2t) 是否满足条件。
  • 注意到 gcd ⁡ ( a b m o d    n , n ) = gcd ⁡ ( a b , n ) ≥ gcd ⁡ ( a , n ) , 1 ≤ a , b < n \gcd(ab \mod n, n) =\gcd(ab, n) \ge \gcd(a,n), 1 \le a,b < n gcd(abmodn,n)=gcd(ab,n)gcd(a,n),1a,b<n,实际实现时我们并不需要每次都暴力求 gcd ⁡ \gcd gcd,而可以将若干次检查合为一次,以减小时间常数。
  • 上述步骤许多部分基于估计,实际上并不严谨,但 Pollard Rho \text{Pollard Rho} Pollard Rho 算法在实际环境中运行得相当不错。
using std::vector;
typedef long long ll;
typedef __int128 s128;
 
template <class T>
inline T Abs(T x) {return x < 0 ? -x : x;}
template <class T>
inline void CkMax(T &x, T y) {x < y ? x = y : 0;}
 
inline ll quick_pow(ll x, ll k, ll mod)
{
	ll res = 1;
	while (k)
	{
		if (k & 1)
			res = (s128)res * x % mod;
		x = (s128)x * x % mod;
		k >>= 1;
	}
	return res;
}

const int pri[] = {2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37};

inline bool millerRabin(ll n)
{
	if (n < 3 || !(n & 1))
		return n == 2;
	ll u = n - 1;
	int t = 0;
	while (!(u & 1))
		u >>= 1, ++t;
	for (int i = 0; i < 12; ++i)
	{
	    if (pri[i] >= n)
	    	break ;
		ll x = pri[i];
		ll res = quick_pow(pri[i], u, n);
		if (res == 1)
			continue ;
		int j = 0;
		for (j = 0; j < t; ++j)
		{
			if (res == n - 1)
				break ;
			res = (s128)res * res % n;		
		}
		if (j >= t)
			return false;
	}
	return true;
}

inline ll f(ll x, ll c, ll n)
{
	return ((s128)x * x + c) % n; 
}

inline ll pollardRho(ll n)
{
	ll rx = 0, x = 0, d, val = 1;
	ll c = rand() % (n - 1) + 1;
	for (int t = 1; ; t <<= 1, rx = x, val = 1)
	{
		for (int k = 1; k <= t; ++k)
		{
			x = f(x, c, n);
			val = (s128)val * Abs(x - rx) % n;
			if ((k % 127) == 0)
			{
				d = std::__gcd(val, n);
				if (d > 1)
					return d;
			}
		}
		ll d = std::__gcd(val, n);
		if (d > 1)
			return d;
	}
}

inline void Factor(ll n, vector<ll> &fac, int cnt)
{   // fac 存储所有质因子,顺序混乱 
    if (n == 1)
       return ;
    if (millerRabin(n))
    {
        while (cnt--)
            fac.push_back(n);
        return ;
    }
    ll p = pollardRho(n); // 求出 n 的某个非平凡正因子 
    int _cnt = 0;
    while (n % p == 0)
        n /= p, ++_cnt;
    Factor(n, fac, cnt);
    Factor(p, fac, cnt * _cnt);
}

类欧几里得

  • a , b , c , n a,b,c,n a,b,c,n 同阶,设

f ( a , b , c , n ) = ∑ i = 0 n ⌊ a i + b c ⌋ g ( a , b , c , n ) = ∑ i = 0 n i ⌊ a i + b c ⌋ h ( a , b , c , n ) = ∑ i = 9 n ⌊ a i + b c ⌋ 2 \begin{aligned} f(a, b, c, n) &= \sum \limits_{i = 0}^{n}\lfloor \frac{ai + b}{c} \rfloor \\ g(a, b, c, n) &= \sum \limits_{i = 0}^{n}i\lfloor \frac{ai + b}{c} \rfloor \\ h(a, b, c, n) &= \sum \limits_{i = 9}^{n}\lfloor \frac{ai + b}{c}\rfloor^2 \end{aligned} f(a,b,c,n)g(a,b,c,n)h(a,b,c,n)=i=0ncai+b=i=0nicai+b=i=9ncai+b2

  • 其中 f f f 可以独立计算,而 g , h g,h g,h 的计算则会导致 f , g , h f,g,h f,g,h 的交错递归,故将三项合并计算,时间复杂度 O ( log ⁡ n ) \mathcal O(\log n) O(logn)
inline void add(ll &x, ll y)
{
	x += y;
	x >= P ? x -= P : 0;
}

struct data 
{
	data() {f = g = h = 0;}
	ll f, g, h;
};

data calc(ll n, ll a, ll b, ll c) 
{
	ll ac = a / c, bc = b / c, m = (a * n + b) / c, n1 = n + 1, n21 = n << 1 | 1;
	data d;
	if (a == 0) 
	{  	
		d.f = bc * n1 % P;
		d.g = bc * n % P * n1 % P * i2 % P;
		d.h = bc * bc % P * n1 % P;
    	return d;
	}
	if (a >= c || b >= c) 
	{
	    d.f = (n * n1 % P * i2 % P * ac + bc * n1) % P;
    	d.g = (ac * n % P * n1 % P * n21 % P * i6 + bc * n % P * n1 % P * i2) % P;
    	d.h = (ac * ac % P * n % P * n1 % P * n21 % P * i6 +
          	   bc * bc % P * n1 + ac * bc % P * n % P * n1) % P;
		data e = calc(n, a % c, b % c, c); 
        add(d.f, e.f), add(d.g, e.g);
		d.h = (d.h + e.h + 2ll * bc % P * e.f + 2ll * ac % P * e.g) % P;	
    	return d;
	}
	data e = calc(m - 1, c, c - b - 1, a);
	d.f = (n * m + P - e.f) % P;
	d.g = (m * n % P * n1 + P - e.h + P - e.f) % P * i2 % P;
	d.h = (n * m % P * (m + 1) + (P - e.g << 1) + (P - e.f << 1) + P - d.f) % P;
	return d;
}

Lucas 定理

  • 用于快速计算 ( n m )   m o d   p \binom{n}{m} \bmod p (mn)modp
  • p p p 为质数,预处理阶乘 O ( p ) \mathcal O(p) O(p),单次询问 O ( log ⁡ n ) \mathcal O(\log n) O(logn)
inline int C(int n, int m)
{ 	//fac 和 ifac 预处理到 [0, mod) 
	if (n < 0 || m < 0 || n < m)
		return 0;
	return 1ll * fac[n] * ifac[n - m] % mod * ifac[m] % mod;
}

inline int Lucas(int x, int y)
{   //调用该函数返回结果 
	if (!y) return 1;
	return 1ll * C(x % mod, y % mod) * Lucas(x / mod, y / mod) % mod; 
}
  • p p p 不为质数,将其质因数分解,求出 ( n m )   m o d   p i c i \dbinom{n}{m} \bmod {p_i}^{c_i} (mn)modpici 后用 CRT 合并,时间复杂度 O ( ∑ ( p i c i + log ⁡ n + log ⁡ p ) ) \mathcal O(\sum (p_i^{c_i} + \log n + \log p)) O((pici+logn+logp))
  • 对质数 P P P,考虑到 ( n m )   m o d   P k = n ! P x m ! P y ( n − m ) ! P z P x − y − z   m o d   P k \dbinom{n}{m} \bmod P^k = \dfrac{\dfrac{n!}{P^x}}{\dfrac{m!}{P^y}\dfrac{(n - m)!}{P^z}}P^{x - y- z} \bmod P^k (mn)modPk=Pym!Pz(nm)!Pxn!PxyzmodPk,其中 x , y , z x,y,z x,y,z 为对应阶乘内 P P P 的次数,因而只需求 f ( n ) = n ! P x   m o d   P k f(n) = \dfrac{n!}{P^x} \bmod P^k f(n)=Pxn!modPk 以及对应的 x x x
  • 可将 n ! n! n! 拆解为以下形式:
    n ! = P ⌊ n P ⌋ ⌊ n P ⌋ ! ( ∏ i = 1 , i ≢ 0 (   m o d   P ) P k i ) ⌊ n P k ⌋ ( ∏ i = P k ⌊ n P k ⌋ , i ≢ 0 (   m o d   P ) n i ) n! = P^{\left\lfloor\frac{n}{P}\right\rfloor} \left\lfloor\dfrac{n}{P}\right\rfloor!\left(\prod\limits_{i = 1,i\not \equiv 0(\bmod P)}^{P^k}i\right)^{\left\lfloor\frac{n}{P^k}\right\rfloor}\left(\prod\limits_{i = P^k\left\lfloor\frac{n}{P^k}\right\rfloor,i\not \equiv 0(\bmod P)}^{n}i\right) n!=PPnPn! i=1,i0(modP)Pki Pkn i=PkPkn,i0(modP)ni
  • 两个乘积项显然可以预处理关于 P k P^k Pk 的循环节,循环节中前 i i i 项的积(模 P P P 为 0 视为乘 1)记为 s u m ( i ) sum(i) sum(i),则 f ( n ) f(n) f(n) 有以下递推式,且 x x x 可在递推中顺便计算:

f ( n ) = f ( ⌊ n P ⌋ ) s u m ⌊ n P k ⌋ ( P k − 1 ) s u m ( n   m o d   P k ) x = x + ⌊ n P ⌋ \begin{aligned} f(n) &= f\left(\left\lfloor\dfrac{n}{P}\right\rfloor\right)sum^{\left\lfloor\frac{n}{P^k}\right\rfloor}\left({P^k - 1}\right)sum\left(n \bmod P^k\right) \\ x &= x + \left\lfloor\dfrac{n}{P}\right\rfloor \end{aligned} f(n)x=f(Pn)sumPkn(Pk1)sum(nmodPk)=x+Pn

inline int quick_pow(int x, ll k, int mod)
{
	int res = 1;
	while (k)
	{
		if (k & 1) res = 1ll * res * x % mod;
		x = 1ll * x * x % mod; k >>= 1;
	}
	return res;
}

inline ll ex_gcd(ll a, ll b, ll &x, ll &y)
{
	if (!b)
	{
		x = 1; y = 0;
		return a;
	}
	ll e = ex_gcd(b, a % b, y, x);
	y -= a / b * x;
	return e;
}

inline int query_inv(int a, int b)
{
	ll x, y;
	ex_gcd(a, b, x, y);
	return (x % b + b) % b;
}

inline int solve_fac(ll n, int p, int mul_p, int opt)
{
	if (n < p)
		return fac[n];
	tot += n / p * opt;
	return 1ll * quick_pow(sum[mul_p - 1], n / mul_p, mul_p) * sum[n % mul_p] % mul_p
		       * solve_fac(n / p, p, mul_p, opt) % mul_p;
}

inline int C(ll n, ll m, int p, int mul_p)
{
	tot = 0;
	sum[0] = 1;
	for (int i = 1; i < mul_p; ++i)
	{
		if (i % p != 0) 
			sum[i] = 1ll * sum[i - 1] * i % mul_p;
		else 
			sum[i] = sum[i - 1];
	}
	fac[0] = 1;
	for (int i = 1; i < p; ++i)
		fac[i] = 1ll * fac[i - 1] * i % mul_p;

	return 1ll * solve_fac(n, p, mul_p, 1)
	           * query_inv(solve_fac(m, p, mul_p, -1), mul_p) % mul_p
	           * query_inv(solve_fac(n - m, p, mul_p, -1), mul_p) % mul_p
			   * quick_pow(p, tot, mul_p) % mul_p;
}

inline int exLucas(ll n, ll m, int p)
{
	int x = p, ans = 0;
	for (int i = 2, im = sqrt(x); i <= im && i <= x; ++i)
		if (x % i == 0)
		{
			int tmp = x;
			x /= i;
			while (x % i == 0)
				x /= i;
			tmp /= x;
			ans = (ans + 1ll * C(n, m, i, tmp) * (p / tmp) % p * query_inv(p / tmp, tmp)) % p;
		}
	if (x > 1)
		ans = (ans + 1ll * C(n, m, x, x) * (p / x) % p * query_inv(p / x, x)) % p;
	return ans;
}

本原勾股数组

  • 定义 三元组 ( a , b , c ) (a,b,c) (a,b,c) 满足 a 2 + b 2 = c 2 a^2 + b^2 = c^2 a2+b2=c2 ( a , b , c ) = 1 (a,b,c) = 1 (a,b,c)=1
  • 性质1 a , b a,b a,b 奇偶不同,且 c c c 为奇数

    证明 分情况讨论:

    • a , b a,b a,b 均为偶数, c c c 也为偶数, ( a , b , c ) = 2 (a,b,c) = 2 (a,b,c)=2,与定义矛盾。

    • a , b a,b a,b 均为奇数, c c c 为偶数,设 a = 2 x + 1 , b = 2 y + 1 , c = 2 z a = 2x + 1, b = 2y + 1, c = 2z a=2x+1,b=2y+1,c=2z,则
      a 2 + b 2 = c 2 ⇔ ( 2 x + 1 ) 2 + ( 2 y + 1 ) 2 = ( 2 z ) 2 ⇔ 2 x 2 + 2 y 2 + 2 x + 2 y + 1 = 2 z 2 a^2 + b^2 = c^2 \Leftrightarrow (2x + 1)^2 + (2y + 1)^2 = (2z)^2 \Leftrightarrow 2x^2+2y^2 + 2x + 2y + 1 = 2z^2 a2+b2=c2(2x+1)2+(2y+1)2=(2z)22x2+2y2+2x+2y+1=2z2
      等式左边为奇数,右边为偶数,矛盾。

    • 故只能满足 性质1

  • 性质2 a , b , c a,b,c a,b,c 两两互素。

    证明 反证法,若 ( a , b ) = d ( d > 1 ) (a,b) = d(d > 1) (a,b)=d(d>1),设 a = d x , b = d y a = dx, b = dy a=dx,b=dy,则
    c = ( d x ) 2 + ( d y ) 2 = d x 2 + y 2 c = \sqrt{(dx)^2 + (dy)^2} = d \sqrt{x^2 + y^2} c=(dx)2+(dy)2 =dx2+y2
    c c c 也含有约数 d d d ,与定义矛盾,其余情况同理。

  • 性质3 假定 a a a 为奇数, b b b 为偶数, c − b c - b cb c + b c + b c+b 均为平方数。

    证明 d ∣ ( c − b , c + b ) d | (c - b, c + b) d(cb,c+b),有

    • d ∣ ( c − b + c + b ) ⇔ d ∣ 2 c d | (c - b + c + b) \Leftrightarrow d | 2c d(cb+c+b)d∣2c
    • d ∣ [ c + b − ( c − b ) ] ⇔ d ∣ 2 b d | [c + b - (c - b)] \Leftrightarrow d | 2b d[c+b(cb)]d∣2b

    性质 2 b , c b,c b,c 互素 且 c − b c - b cb 为奇数,一定有 d = 1 d = 1 d=1,则 ( c − b , c + b ) = 1 (c - b, c + b) = 1 (cb,c+b)=1

    又因为 ( c − b ) ( c + b ) = a 2 (c - b)(c + b) = a^2 (cb)(c+b)=a2 c − b c - b cb c + b c + b c+b 均为平方数。

  • 性质3 c − b = t 2 c - b = t^2 cb=t2 c + b = s 2 c + b = s^2 c+b=s2 s > t ≥ 1 s > t \ge 1 s>t1,解得:
    a = s t , b = s 2 − t 2 2 , c = s 2 + t 2 2 a = st, b = \frac{s^2 - t^2}{2}, c = \frac{s^2 + t^2}{2} a=st,b=2s2t2,c=2s2+t2
  • 对于单位圆 x 2 + y 2 = 1 , x = b c , y = a c x^2 + y^2 = 1, x = \frac{b}{c}, y = \frac{a}{c} x2+y2=1,x=cb,y=ca,代入上式并取 m = t s = tan ⁡ θ 2 , m ∈ R m = \frac{t}{s} = \tan\frac{\theta}{2}, m \in \mathbb R m=st=tan2θ,mR,可得:
    ( x , y ) = ( 1 − m 2 1 + m 2 , 2 m 1 + m 2 ) = ( 1 − tan ⁡ 2 θ 2 1 + tan ⁡ 2 θ 2 , 2 tan ⁡ θ 2 1 + tan ⁡ 2 θ 2 ) = ( cos ⁡ θ , sin ⁡ θ ) (x,y) = (\frac{1 - m^2}{1 + m^2}, \frac{2m}{1 + m^2}) = (\frac{1 - \tan^2\frac{\theta}{2}}{1 + \tan^2\frac{\theta}{2}}, \frac{2\tan\frac{\theta}{2}}{1+\tan^2\frac{\theta}{2}}) = (\cos\theta, \sin\theta) (x,y)=(1+m21m2,1+m22m)=(1+tan22θ1tan22θ,1+tan22θ2tan2θ)=(cosθ,sinθ)

威尔逊定理

  • p 为素数 ⇔ ( p − 1 ) ! + 1 ≡ 0 ( m o d    p ) p为素数 \Leftrightarrow (p - 1)! + 1\equiv 0 (\mod p) p为素数(p1)!+10(modp)

证明

  • 先证 ( p − 1 ) ! + 1 ≡ 0 ( m o d    p ) ⇒ p 为素数 (p - 1)! + 1\equiv 0 (\mod p) \Rightarrow p为素数 (p1)!+10(modp)p为素数
  • 反证法,假设 p = a b , a , b ∈ ( 1 , p ) p = ab, a,b\in(1,p) p=ab,a,b(1,p)
    • a ≠ b a \not = b a=b,显然有 ( a b − 1 ) ! ≡ 0 ( m o d    a b ) (ab - 1)! \equiv 0(\mod ab) (ab1)!0(modab),矛盾。
    • a = b = p a = b = \sqrt p a=b=p ,当 p ≤ 4 p \le 4 p4 时可直接验证, p > 4 p > 4 p>4 必有 p , 2 p < p \sqrt p, 2\sqrt p < p p ,2p <p,因此 ( p − 1 ) ! ≡ 0 ( m o d    p ) (p - 1)! \equiv 0 (\mod p) (p1)!0(modp),同样矛盾。
  • 再证 p 为素数 ⇒ ( p − 1 ) ! + 1 ≡ 0 ( m o d    p ) p 为素数 \Rightarrow (p - 1)!+1 \equiv 0(\mod p) p为素数(p1)!+10(modp)
    • 对于集合 S = { x ∣ 1 ≤ x < p , x ∈ N } S = \{x|1 \le x < p,x\in \mathbb N\} S={x∣1x<p,xN},不断取出 x , y ∈ S , x ≠ y , x y ≡ 1 ( m o d    p ) x,y\in S,x\not = y, xy \equiv 1(\mod p) x,yS,x=y,xy1(modp) 配成一对,由逆元性质可知,最后必定剩下 1 , p − 1 1, p - 1 1,p1,因此有 ( p − 1 ) ! ≡ − 1 ( m o d    p ) (p - 1)! \equiv -1(\mod p) (p1)!1(modp)
好的,我知道你的问题了。首先,让我们来了解一下Numerical optimization书中Algorithm 3.2Algorithm 3.3的matlab例程的基本情况。 这两个算法都是针对最小化二元函数的,Algorithm 3.2是基于梯度下降法的,Algorithm 3.3是基于共轭梯度法的。它们的区别在于梯度下降法每次迭代都需要计算函数的梯度,而共轭梯度法则不需要。 下面是这两个算法的MATLAB代码示例: Algorithm 3.2(梯度下降法)的MATLAB代码: ```matlab function [x, fval, k] = gradient_descent(f, x0, eta, epsilon, max_iter) % f: 目标函数 % x0: 初始点 % eta: 步长 % epsilon: 精度 % max_iter: 最大迭代次数 % x: 最优解 % fval: 最优解对应的函数值 % k: 实际迭代次数 x = x0; k = 0; while k < max_iter grad = [diff(f, 'x1'); diff(f, 'x2')]; % 计算梯度 x_new = x - eta * grad; % 梯度下降 if norm(x_new - x) < epsilon % 判断是否满足精度要求 break end x = x_new; k = k + 1; end fval = double(subs(f, {'x1', 'x2'}, x)); end ``` Algorithm 3.3(共轭梯度法)的MATLAB代码: ```matlab function [x, fval, k] = conjugate_gradient(f, x0, epsilon, max_iter) % f: 目标函数 % x0: 初始点 % epsilon: 精度 % max_iter: 最大迭代次数 % x: 最优解 % fval: 最优解对应的函数值 % k: 实际迭代次数 x = x0; grad = [diff(f, 'x1'); diff(f, 'x2')]; % 计算梯度 d = -grad; k = 0; while k < max_iter alpha = subs(f, {'x1', 'x2'}, x + d)' * d / (d' * subs(diff(grad, 'x1'), {'x1', 'x2'}, x + d) + d' * subs(diff(grad, 'x2'), {'x1', 'x2'}, x + d)); % 计算步长 x_new = x + alpha * d; % 更新x if norm(x_new - x) < epsilon % 判断是否满足精度要求 break end grad_new = [diff(f, 'x1'); diff(f, 'x2')]; % 计算新的梯度 beta = (grad_new' * grad_new) / (grad' * grad); % 计算beta d = -grad_new + beta * d; % 更新d x = x_new; grad = grad_new; k = k + 1; end fval = double(subs(f, {'x1', 'x2'}, x)); end ``` 在这两个例程中,f代表的是目标函数,x0代表的是初始点,eta代表的是步长(Algorithm 3.2),epsilon代表的是精度,max_iter代表的是最大迭代次数。最后返回的是最优解x、最优解对应的函数值fval和实际迭代次数k。 以上是关于Numerical optimization书中Algorithm 3.2Algorithm 3.3的matlab例程的基本介绍和示例代码。
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值