数学-剩余系
欧几里得算法与扩展欧几里得算法
gcd. 用于求最大公约数.
int gcd(int a,int b)
{ return !a ? b : gcd(b%a,a); }
扩欧. 用于求方程 $ax+by=gcd(a,b)$ 的解$(x,y)$ .
struct value{ int x,y; value(int x,int y):x(x),y(y){} }; value ExtEuclid(int a,int b) { if(a==0) return value(0,1); value r=ExtEuclid(b%a,a); return value(r.y-b/a*r.x,r.x); }
有解的充要条件:
$$ gcd(a,b)|c $$
它的解为
$$ (x,y)\frac{c}{gcd(a,b)} $$
$(x,y)$ 由 $ExEuclid(a,b)$ 的返回值给出.
非递归写法
gcd
int gcd(int a,int b)
{ while(a) b=b%a,swap(a,b); return b; }
扩欧的不会.....
求剩余系中某数的乘法逆元
考虑剩余系下的除法 $ a/b=ab^{-1} \space \space (mod \space MOD) $ ,有 $ a^{-1} a=1 \space \space (mod \space MOD) $
这个算法求出$a^{-1}$ ........
注意只有模数为质数时,整个剩余系的除零以外其它元素才都会有逆元.
int getrev(int p)
{ return p==0 ? -1 : (ExtEuclid(p,MOD).x%MOD+MOD)%MOD; }
具体的,要求 $ ax \equiv 1 \space \space (mod \space MOD) $ 中的x,有 $ ax+MODy=1 $
如果a与MOD的gcd不为1,则逆元不存在.
直接用扩欧来算(x,y)就行了......
线性同余方程组
考虑方程组
$$ \left\{\begin{matrix}
{x \equiv a_1 \; ( mod \; m_1 ) }\\
{x \equiv a_2 \; ( mod \; m_2 ) }\\
{ ...... }\\
{x \equiv a_n \; ( mod \; m_n ) }
\end{matrix}\right. $$
其中
$$ \forall_{i,j} \quad a_i \equiv a_j \; (mod \; gcd(m_i,m_j)) \quad (有解条件) $$
并且
$$ \forall_{i,j} \quad gcd(m_i,m_j)=1 \quad (中国剩余定理应用条件) $$
- 使用中国剩余定理构造解x的方法如下
设 $ M=\prod m_i $
设 $ M_i = M/m_i $
设 $ t_i = M_i ^{-1} (mod \; m_i) $ 即 $t_i$ 为 $M_i$ 在模 $m_i$ 剩余系下的逆元.
则 $ x=\sum t_i a_i M_i $
AC VIJOS 1164 裸的中国剩余定理.
1 int n;
2 ll a[20],m[20];
3 ll t;
4
5 struct value{ ll x,y; value(ll x,ll y):x(x),y(y){} };
6 value ExEuclid(ll a,ll b)
7 {
8 if(!a) return value(1,1);
9 value r=ExEuclid(b%a,a);
10 return value(r.y-b/a*r.x,r.x);
11 }
12 ll rev(ll k,ll m)
13 { return (ExEuclid(k%m,m).x%m+m)%m; }
14
15 int main()
16 {
17 n=getint();
18 for(int i=0;i<n;i++)
19 scanf("%I64d%I64d",m+i,a+i);
20
21 t=1;
22 for(int i=0;i<n;i++)
23 t*=m[i];
24
25 ll res=0;
26 for(int i=0;i<n;i++)
27 res+=t/m[i]*rev(t/m[i],m[i])*a[i];
28
29 printf("%I64d\n",res%t);
30 return 0;
31 }
注意 $M$ 的值会非常大,考虑是否需要使用高精度.
应用的时候不要混淆解的判定条件以及定理应用条件!
考虑方程组
$$ x \equiv 0 \; (mod \; 2) \\ x\equiv 2 \;(mod \; 4) $$
显然有解 $x=2,6,....$
而按照上述解法,有 $M_1=4$ ,因而 $M_1 \; mod \; m_1 = 0 $ ,这样 $M_1$ 就不存在逆元 $t_1$ .
- 如何得到一般方程组的解呢?
考虑把方程依次加入.
先看第一个方程 $$ x \; mod \; m_1 = a_1 $$
它的解是 $$ x=a_1+km_1 \; \; , \; k \in N \quad $$
现有解 $$ x_1 \leftarrow a_1+km_1 $$
当我们加入第$i$个方程时,
$$ x_i \leftarrow x_{i-1}+k \times lcm(m_1,m_2,...,m_{i-1}) $$
其中的 $k$ 满足 $$ x_{i-1} + k \times lcm(m_1,m_2,...,m_{i-1}) \; mod \; m_i = a_i $$
亦即 $$ x_{i-1} + k \times lcm(m_1,m_2,...,m_{i-1}) + tm_i = a_i $$
使用扩展欧几里得算出k即可.
为什么算法是正确的? lcm的意义是什么?
@iwtwiioi 指出lcm是必要的,这里没有给出lcm的必要性证明.
呃,这个东西丢到周末去做吧.........在我彻底理解那玩意之后......
update(3.29):于是还是没改......
这里有归纳证明: $$ 若 x_i 满足了前i个方程,则由上述算法得出的 x_{i+1} 满足了前i+1个方程. $$
我们有 $$ k \times lcm(m_1,m_2,...,m_i) \; mod \; m_p = 0 $$ 这是因为 $1 \leq p \leq i$ ,从而 $m_p \mid lcm(m_1,m_2,...,m_i)$
对于前$i$个方程中的某个方程 $p$ ,必定有 $$ x_i \; mod \; m_p = a_p $$ 这是因为 $x_i$ 满足前i个方程.
根据上边的式子以及递推过程我们得到
$$ x_{i+1} \; mod \; m_p \\ = ( x_i + k \times lcm(m_1,m_2,...,m_i)) \; mod \; m_p \\ = x_i \; mod \; m_p + k \times lcm(m_1,m_2,...,m_i) \; mod \; m_p \\ = a_p + 0 \\ = a_p $$
这样 $x_{i+1}$ 就满足前$i$个方程了.我们还剩第$i+1$个方程. 注意我们是按照第 $i+1$ 个方程的约束求出 $k$ 的,很明显 $k$ 的值会使 $x_{i+1}$ 满足方程 $i+1$.
这样就结束了$=\omega=$
一题. AC POJ 2891 解一般的线性模方程组.
1 inline ll modc(ll k,const ll MOD)
2 { return (k%MOD+MOD)%MOD; }
3
4 int n;
5 ll t;
6
7 struct value{ ll x,y; value(ll x,ll y):x(x),y(y){} };
8 value ExEuclid(ll a,ll b)
9 {
10 if(!a) return value(1,1);
11 value r=ExEuclid(b%a,a);
12 return value(r.y-b/a*r.x,r.x);
13 }
14
15 ll rev(ll k,ll m)
16 { return modc(ExEuclid(k%m,m).x,m); }
17
18 ll gcd(ll a,ll b)
19 { while(a) b=b%a,swap(a,b); return b; }
20
21 int main()
22 {
23 while(scanf("%d",&n)>0)
24 {
25 bool ok=1;
26 ll a,m,x;
27 m=getll();
28 a=getll();
29
30 x=a;
31 ll lcm=m;
32
33 for(int i=1;i<n;i++)
34 {
35 m=getll();
36 a=getll();
37
38 if(!ok) continue;
39
40 ll d=gcd(lcm,m);
41
42 if((a-x)%d!=0) { ok=false; continue; }
43
44 ll k=ExEuclid(lcm,m).x;
45 k=(a-x)/d*k%(m/d);
46
47 x+=k*lcm;
48 lcm=m/d*lcm;
49 x=x%lcm;
50 }
51
52 if(!ok) printf("-1\n");
53 else printf("%lld\n",modc(x,lcm));
54 }
55
56 return 0;
57 }
注意精度问题. 能取模就取模......
由于是同余恒等式,切记不要搞错符号......
ExEuclid的变量不是可以随便换的......
再来一题. AC HDU 1573 也是比较裸的同余方程,要求给定范围内解的个数.
没仔细看题导致WA了一个小时 TAT
X是正整数啊..... 我把0算进去了啊...... TAT
View Code
Baby Step Giant Step (BSGS)
考虑有一个完全剩余系:$$a^0,a^1,a^2 \; ... \;,a^{p-1} (mod \; p)$$
我们想求一个 $x$ 使得 $$a^x\; \equiv \; b \; \; (mod \; p) $$
那么,考虑把剩余系按原顺序分成 $m=\lceil{\sqrt{p}}\rceil$ 块,每一块包含元素
$$ a^{km},a^{km+1},...,a^{km+m-1} $$
,即
$$ a^{km},a^{km}a^1,...,a^{km}a^{m-1} $$
我们把 $$ ba^{-1},ba^{-2},...,ba^{-m+1} $$ 放入一个双射的表(hash,或者map).
然后枚举 $k$ 计算 $ a^{km} $ ,在表里面找到与之相等的元素 $ba^{-i}$ ,那么 $km+i$ 就是答案了.......
嗯,这么做只是因为有如下等式 $$ a^{km} \; \equiv \; ba^{-i} \;\;(mod\; p) $$这是显然的吧.....
AC BZOJ 2242
涵盖了常用的数值算法....
task1: 快速幂
task2: 扩欧
task3: 分块(Baby Step Giant Step)
1 #include <cstdio>
2 #include <fstream>
3 #include <iostream>
4
5 #include <cstdlib>
6 #include <cstring>
7 #include <algorithm>
8 #include <cmath>
9
10 #include <queue>
11 #include <vector>
12 #include <map>
13 #include <set>
14 #include <stack>
15 #include <list>
16
17 typedef unsigned int uint;
18 typedef long long int ll;
19 typedef unsigned long long int ull;
20 typedef double db;
21
22 using namespace std;
23
24 inline int getint()
25 {
26 int res=0;
27 char c=getchar();
28 bool mi=false;
29 while(c<'0' || c>'9') mi=(c=='-'),c=getchar();
30 while('0'<=c && c<='9') res=res*10+c-'0',c=getchar();
31 return mi ? -res : res;
32 }
33 inline ll getll()
34 {
35 ll res=0;
36 char c=getchar();
37 bool mi=false;
38 while(c<'0' || c>'9') mi=(c=='-'),c=getchar();
39 while('0'<=c && c<='9') res=res*10+c-'0',c=getchar();
40 return mi ? -res : res;
41 }
42
43 //==============================================================================
44 //==============================================================================
45 //==============================================================================
46 //==============================================================================
47
48 ll MOD;
49
50 struct pl{ ll x,y; pl(ll a,ll b):x(a),y(b){} };
51 pl ExEuclid(ll a,ll b)
52 {
53 if(!a) return pl(1,1);
54 pl r=ExEuclid(b%a,a);
55 return pl(r.y-b/a*r.x,r.x);
56 }
57
58 ll FastMulti(ll a,ll x)
59 {
60 ll res=1;
61 while(x)
62 {
63 if(x&1) (res*=a)%=MOD;
64 (a*=a)%=MOD;
65 x>>=1;
66 }
67 return res;
68 }
69
70 ll gcd(ll a,ll b)
71 { if(a>b)swap(a,b); while(a)b=b%a,swap(a,b); return b; }
72
73 ll rev(ll a)
74 { return (ExEuclid(a,MOD).x%MOD+MOD)%MOD; }
75
76 map<ll,int> g;
77
78 int main()
79 {
80 int T=getint();
81 int type=getint();
82 if(type==1) while(T--)
83 {
84 ll a=getint();
85 ll b=getint();
86 MOD=getint();
87 printf("%lld\n",FastMulti(a,b));
88 }
89 else if(type==2) while(T--)
90 {
91 ll a=getint();
92 ll c=getint();
93 MOD=getint();
94 if(c%gcd(a,MOD)!=0) printf("Orz, I cannot find x!\n");
95 else
96 printf("%lld\n",((c/gcd(a,MOD)*ExEuclid(a,MOD).x)%MOD+MOD)%MOD);
97 }
98 else if(type==3) while(T--)
99 {
100 g.clear();
101 ll a=getint();
102 ll b=getint();
103 MOD=getint();
104 ll lim=(ll)(ceil(sqrt((long db)MOD)));
105 ll c=1;
106 for(int i=0;i<lim;i++)
107 g[(rev(c)*b%MOD+MOD)%MOD]=i,(c*=a)%=MOD;
108 ll d=1;
109 bool found=false;
110 ll res=0;
111 for(int i=0;i<=lim;i++)
112 {
113 map<ll,int>::iterator v=g.find(d);
114 if(v!=g.end())
115 { res=(lim*i+v->second)%MOD; found=true; break; }
116 else (d*=c)%=MOD;
117 }
118 if(!found) printf("Orz, I cannot find x!\n");
119 else printf("%lld\n",(res%MOD+MOD)%MOD);
120 }
121
122 return 0;
123 }
原根
考虑一个模 $p$ 剩余系. 我们知道对于 $\leq a < p-1$ , $a^1,a^2,...,a^{\phi(p)}$ 是 $a$ 在剩余系下的一个循环节.但不一定是最短的.
我们把这个剩余系下循环节长度刚好是 $\phi(p)$ 的数 $r$ 称为这个数的原根.
模 $p$ 剩余系的充要条件: $p=2$ 或 $p=4$ 或 $p=s^k$ 或 $p=2s$ ,其中 $s$ 是奇素数, $k$ 是任意正整数.
如何判断 $x$ 是否是原根? 注意 $x^{\phi(p)}=p$ 是必然的,不论 $x$ 是什么值.
于是乎,直接枚举 $d|\phi(p)$ ,判断 $x^d=p$ 是否成立.如果有除了 $\phi(d)$ 以外的数成立,那么这个数就不是原根.否则就是.
如何求原根? 从 $2$ 开始逐一枚举然后判断....... 原根一般都很小.....
未完待续
模意义下的高斯消元
吐槽
写LaTeX要抓狂了啊啊啊啊啊 $QAQ$
谁告诉我公式左对齐怎么搞Wiki大法好.....
参考及拓展
