数论在ACM中的应用

转自：http://blog.csdn.net/danliwoo/article/details/48827813#t22

整除与剩余

整除在自然数范围内的引入是十分自然的，即要把一个整体平均分为若干份。若要求分得的结果必须是整数，则可能存在多余的情况，即余数。

若a能整除b，可记作：a∣b，否则a∤b

整数是对自然数的扩充，引入了负数的概念。负的余数是一个数学概念，在实际情形中说剩余的东西是负数，往往是可笑的。因为那意味着并没有剩余。利用同余，可以给出这类情况的意义。

整除判断多利用代数恒等变形，可尝试证明如下两题： 
1. 设a>1,m,n>0, 证明：(am−1,an−1)=a(m,n)−1. 
2. 设a>b,gcd(a,b)=1, 证明：(am−bm,an−bn)=a(m,n)−b(m,n).

同余，即除同一个数得到相同的余数。例如8=1∗5+3，13=2∗5+3，可知8，13关于5同余，一般可以记作：8≡13(mod5)

若a和b模d同余，则下列命题等价： 

a和b模d同余⇔存在整数n，使得a=b+nd⇔d∣(a−b)

当a是一个负数时，则将n置为负数。

几个看两眼的定理：

• 若a≡b(modm)且d∣m,则a≡b(modd)；
• 若a≡b(modm),则(a,m)=(b,m)(辗转相除法)；
• ∀1≤i≤n,a≡b(modmi)⇔a≡b(mod[m1,m2,...,mn])
• 若ac≡bc(modm),且(c,m)=d,则a≡b(modmd)

对于相同的模，同余式可以加、减、乘，可以通过上述等价的命题快速证明。在计算整数幂的同余式时可以根据乘法性质对小指数的先进行模运算，再代回原式，可以减小一定的计算量。例如计算22015模13，可以先得到 

24≡16≡3(mod13)

28≡32≡9≡−4(mod13)

212≡−12≡1(mod13)

22015≡211≡−4∗8≡−32≡7(mod13)

有些题目答案过大时，往往要求输出mod一个大数N的答案。可以在中间过程中先mod N，避免超int。

某些数被模除时具有特别的性质，下面看一道自己出的题： 
题目大意：从数字n到数字m(1≤n≤m≤108)之间有多少个数满足奇数位上的数字之和减去偶数位上的数字之和为2。 
由于数据范围不大，可以暴力试试看。机智的会用数位DP。而出题人其实是从数论的角度出发的。已知： 

102k≡100k≡(9∗11+1)k≡1(mod11)

102k+1≡100k∗10≡10≡−1(mod11)

对于任意的正整数 N=a2na2n−1…a2a1¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯(a2n与a2n−1至少有一个不为0），并且设x = 奇数位之和-偶数位之和，则原正整数都可拆解为：

a2na2n−1...a2a1¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯=∑i=12nai⋅10i=∑i=1na2i−1⋅102i−1+∑i=1na2i⋅102i

≡∑i=1na2i−1⋅1+∑i=1na2i⋅(−1)≡x≡d(mod11)

现在要求x=2,则d=2，则N≡d≡2(mod11),可知满足题意的数都有这个性质，再进行判断。但要注意的是d=2时，x不一定等于2，比如13，-9等都能得到d=2. 因此对于mod 11为2的数中仍要进行筛查。

---

素数问题

基本描述

素数仅能被1和自身整除（除1外），可以看作是组成整数的基本单位。 
要知道一个整数的阶乘能被素数p整除多少次，即p的幂为多少次，可利用：

[np]+[np2]+[np3]+...+[npt]

其中t=logpn.

素数定理：定义π(x)为1~x内素数的个数（即欧拉函数），则：

limx→∞π(x)xln(x)=1⇔x→∞时，π(x)∼xln(x)

由此可以估算出在1~x范围内素数的个数。

素数的形成没有成型的规律，但有许多相关的猜想。

• 伯特兰猜想：∀n>1,∃素数p,s.t.n<p<2n
• 孪生素数：存在无穷多素数对p，p+2.
• 哥德巴赫猜想：任何偶数都可以拆分成两个素数的和。
• 构造猜想：存在无穷多个素数满足p=n2+1

素数筛法

• 埃拉托斯尼斯筛法 
  先将2~N内所有的数标记为素数，从最小的素数开始筛去其倍数，再找到下一个素数，依次筛去非素数的数，剩余的即为素数。 
  需要筛的素数范围只需在2~N−−√.

int num[N] = {0,0,1}, prime[N] = {2}, pr = 1;
void set()
{
    for(int i = 2;i < N;i++)
        num[i] = 1;
    for(int i = 2;i <= sqrt((double)N);i++)
        if(num[i]) for(int j = i*i;j < N;j += i)
            num[j] = 0;
    for(int i = 3;i < N;i += 2)
        if(num[i])
            prime[pr++] = i;
};
1
2
3
4
5
6
7
8
9
10
11
12
1
2
3
4
5
6
7
8
9
10
11
12

在POJ 2689，可以节省空间大小。

• 6N±1筛法 
  由于6N，6N+2，6N+3，6N+4（N>0）均不是素数，故在筛素数时可以直接从6N±1中判断和筛选素数。

素数判定

• 朴素的暴力判断算法O(n−−√) 
  从0~n−−√依次尝试能否整除
• 素数筛法后直接判定O(1)
• Lucas-Lehmer判定法（只对梅森数1作用） 
  设p是素数，第p个梅森数为Mp=2p−1，r1=4.对于k≥2,rk≡rk−12−2(modMp),(0≤rk<Mp),得到序列{rk} 
  
  Mp是素数 ⇔ rp−1≡0(modMp)
• Miller-Robin测试 
  判断n是否为素数 
  
  • 由费马小定理，取与n互质的a作为底数，验证an−1≡1(modn)是否成立。
  • 当上式成立时，n是基于a的伪素数。
  • 多取几个小于n的不同的a，反复验证，若都成立，则认为n是素数。取5次即可保证比较高的正确率。
  • 不适用范围：卡麦克尔数2。
• 二次探测 
  对卡麦尔数（伪素数）进行测试

二次探测定理：对素数p,满足x2≡1(modp)的小于p的正整数解x只有1或p−1.

对n-1除2，检验an−12≡1 or n−1(modn)是否成立,一旦不成立则可认定n不是素数。反复除2直到除不尽为止。

//从n-1的最大奇因子为指数开始判定
bool recheck(LL a, LL n, LL m, LL j)
{
    LL d = po(a, m, n);//对大数运算时里面的乘法要另外算mult_mod
    if(d == 1 || d == n-1)
        return true;
    for(int i = 0;i < j;i++)
    {
        d = d*d % n;
        if(d == 1 || d == n-1)
            return true;
    }
    return false;
}
bool Miller_Rabin(LL n)
{
    if( n < 2)return false;
    if( n == 2)return true;
    if( (n&1LL) == 0)return false;//偶数
    srand((unsigned) time (NULL));
    LL m = n-1, j = 0;
    while( !(m & 1LL) )
    {
        m >>= 1;
        j++;
    }
    for(int i = 0;i < 5;i++)
    {
        LL a  = rand()%(n-2) + 2;
        if( !recheck(a,n,m,j) )
            return false;
    }
    return true;
}
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34

剩余系配对

素数p的完全剩余系为0~p-1，其中12≡(p−1)2≡1(modp)，其他的数2,3,...,p-2可以两两配对(i,j)，满足ij≡1(modp).可以先从中任取一个数a,方程ax≡1(modp)必有唯一解。

---

整数分解

将整数化为质因子的幂的乘积的唯一形式。常用试除法、筛选法，比较简单不举例了。对于非常大的数而言，两者用处不大。

Pollard ρ 整数分解法

原理还不太懂： 
1. 生成两个整数a、b，计算p=(a-b,n)，直到p不为1或a,b出现循环为止。 
2. 若p=n，则n为质数；否则为n的一个约数。

分解n的具体步骤如下： 
1. 选取一个小的随机数x1，迭代生成xi=x2i−1+k，取k=1，若序列出现循环则退出。 
2. 计算p=gcd(xi−1−xi,n)直到p>1,否则返回上一步。 
3. 若p=n,则n为素数。否则p为n的约数，继续分解p和n/p

//找出一个因子
LL pollard_rho(LL x,LL c)
{
    LL i = 1, k = 2;
    srand(time(NULL));
    LL x0 = rand()%(x-1) + 1;
    LL y = x0;
    while(1)
    {
        i++;
        x0 = (mult_mod(x0,x0,x) + c)%x;//迭代公式为(x0*x0+c)%x
        LL d = gcd(y - x0,x);//gcd返回一个绝对值
        if( d != 1 && d != x)
            return d;
        if(y == x0)
            return x;
        if(i == k)
        {
            y = x0;
            k += k;
        }
    }
}

//对n进行素因子分解，存入factor. k设置为107左右即可
void findfac(LL n,int k)
{
    if(n == 1)
        return;
    if(Miller_Rabin(n))
    {
        factor[tol++] = n;
        return;
    }
    LL p = n;
    int c = k;
    while( p >= n)
        p = pollard_rho(p,c--);//c即k表示最大搜索次数？
    //值变化，防止死循环k
    findfac(p,k);
    findfac(n/p,k);
}
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42

如POJ 1811 1061

---

扩展欧几里得

辗转相除法可以来求a，b两个数的最大公约数。而辗转相除的过程中，附带也可以得到满足ax+by=c(△)的解系(X,Y)

先来解决ax+by=gcd(a,b)(∗)这一特殊情况。

1. 最简情况：b=0时，gcd(a,b)=a，方程化为ax=a,则一组特殊解为(1,0)
2. 当b≠0时，由辗转相除法
   gcd(a,b)=gcd(b,a%b)
   可以改写原方程得到 
   
   bx′+(a%b)y′=gcd(b,a%b)
   又a%b=a−⌊ab⌋×b，代入前一方程，得到 
   
   ay′+b(x′−⌊ab⌋y′)=gcd(a,b)
   与方程(*)对比ax+by=gcd(a,b)，可知：x=y′，y=(x′−⌊ab⌋y′)得到一组递推关系。
3. 方程依次向下递归可以得到最简情况，得到解（1，0），再由递推关系回代到上一组解，最终可以得到方程（*）的解。 注意到，这是一组特解，并不唯一。

//解方程ax+by=gcd(a,b) 返回gcd(a,b) 得到一组特解(x,y)
int extend_Euclid(int a, int b, int &x, int &y){
    if(b==0){
    x = 1; y = 0;
    return a;
    }
    int r = extend_Euclid(b, a%b, y, x);
    y -= a/b*x;
    return r;
}
1
2
3
4
5
6
7
8
9
10
1
2
3
4
5
6
7
8
9
10

• 回到方程(△)首先，判断gcd(a,b)是否是c的约数。若非，则方程无解。

• 下面考虑 gcd(a,b)∣c，则方程有解。 
  由扩展欧几里得算法算出方程ax+by=gcd(a,b)的特解(x0，y0) 
  原方程(△)的特解(X0，Y0)=(x0,y0)⋅cgcd(a,b) 
  为保证整数解则ΔY=aΔXb∈Z,则ΔX=bgcd(a,b) 
  原方程(△)的解系为{X=X0+bgcd(a,b)⋅K[2ex]Y=Y0−agcd(a,b)⋅KK为任意整数。

• 令K从0取到[gcd(a,b)−1],求出方程所有解中的代表元素共有gcd(a,b)个。是X模b意义下的所有解。其中最小非负整数解x1=(x0⋅cgcd(a,b)%ran+ran)%ran,其中ran=bgcd(a,b)（在模ran范围内有唯一的解）x1是模ran范围内最小的正整数解，又ran∣b则必是模b范围内最小的正整数解。

• 注意：解系中解的间隔与c无关。也可以从坐标系上的直线看出来。直线与网格的交点为一组整数解，c只决定了直线的平移位置，a，b决定了斜率。

拉梅定理：用欧几里得算法计算两个正整数的最大公因子时，所需的除法次数不会超过两个整数中较小的那个十进制数的位数的5倍。

推论：求两个正整数a,b,a>b的最大公因子需要O(log2a)3次的位运算。

扩展欧几里得的应用 
主要有：求解不定方程、模的逆元、同余方程。如POJ 1061，只要化为ax+by=c的形式即可求解。

---

解二元一次线性方程

注意：为保证gcd(a,b)为正，要让a，b均为正，避免错误，此时解应该反号。

//求解ax+by=c 返回0为无解,否则返回gcd(a,b)个解,用X[],Y[]存;
int X[N], Y[N];
int equation(int a, int b, int c)
{
    int x, y;
    int g = extend_Euclid(a, b, x, y);
    if(c % g)
        return 0;    //表示无解
    x *= c/g, y *= c/g;
    for(int k = 0;k < g;k++)
    {
        X[k] = (x+b/g*k)%b;
        Y[k] = (c-a*X[k])/b;
    }
    return g;
}
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16

---

解一元线性同余方程

Thm 对于方程ax≡b(modm),其中a,b,m∈Z and m>0,(a,m)=d.如果d∣b,则方程恰有d个模m不同余的解。否则方程无解。

容易发现，方程ax≡b(modm)可以写为方程ax+my=b的形式，则利用前面得到的结论就很容易理解了。

//求解ax=b(mod m) 返回0为无解，否则返回gcd(a,m)个mod m意义下的解,用X[]存
int mod(int a, int b, int m)
{
    return equation(a, m, b);
}
1
2
3
4
5
1
2
3
4
5

---

求逆元

特别地，当有(a,m)=1时，ax≡1(modm)得到的x即是a模m的逆。

---

解一元线性同余方程组

当有方程组 

⎧⎩⎨⎪⎪⎪⎪⎪⎪x≡b1(modm1)x≡b2(modm2)...x≡bn(modmn)

可知每个方程必有满足的一个特解，并存在模mi的解系 

{x=b1+m1y1x=b2+m2y2

相减后得到： 

m1y1−m2y2=b2−b1

化为二元一次线性方程，得到特解y01,再化为模lcm(m1,m2)意义下的特解(保证了在原来两个方程的解系中)，代入到第一个方程，则两个方程等价为

x≡b1+m1y10(modlcm(m1,m2))

再继续进行两两运算，直至最后得到满足的答案。

//在模m范围内只有可能有一个解或无解 要求最小整数解时equation()内只需保存一个解代回
LL X[N], Y[N];
int equation(int a, int b, int c)
{
    int x, y;
    int g = extend_Euclid(a, b, x, y);
    if(c % g)
        return 0;    //表示无解
    int ran = b/g;
    X[0] = (c/g*x%ran+ran)%ran;    //只需存最小整数解
    return g;
}
int mo[N], bo[N];
int eq_set()
{
    int n, b1, m1, b2, m2, m, r = 1;
    scanf("%lld", &n);
    for(int i = 0;i < n;i++)
        scanf("%lld", &mo[i]);
    for(int i = 0;i < n;i++)
        scanf("%lld", &bo[i]);
    m1 = mo[0], b1 = bo[0];
    for(int i = 1;i < n;i++)
    {
        m2 = mo[i], b2 = bo[i];
        r = equation(m1, m2, b2-b1);
        if(!r)
            return -1;  //返回无解
        b1 += m1*X[0];
        m1 *= m2/r;
        b1 %= m1;
    }
    return (b1%m1+m1-1)%m1+1; //返回正整数解，0被化为最小公倍数
}
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34

如POJ 1006，由题意可列出满足的方程组，再求出相应的答案。当结果为非正数时利用同余性质化为正数。

---

中国剩余定理

Thm 若m1,m2,...,mr是两两互素的正整数，则同余方程组（同上面的解一元线性同余方程组），x≡ai(modmi)有模M=m1m2...mn的唯一解。

令Ni=Mmi,则(Ni,mi)=1,故存在xi,满足Nixi+miyi=1,即Nixi≡1(modmi)，易得

∑i=1raiNixi≡ai(modmi)

即左式为方程组的解。

//解方程组x=ai(mod mi) mi之间两两互质
int China(int r)
{
    int M = 1, ans = 0;
    for (int i = 0; i < r; ++i)
        M *= m[i];
    for(int i = 0;i < r;i++)
    {
        int N = M/m[i];
        int x, y;
        extend_Euclid(N, m[i], x, y);
        ans = (ans+a[i]*N*x)%M;
    }
    ans = (ans - a[r])%M;
    return (ans+M)%M;
}
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16

---

解多元线性同余方程

Thm 多元线性同余方程a1x1+a2x2+...+anxn+b≡0(modm)有解⇔(a1,a2,...,an,m)∣b 
若有解，则解（模m的剩余类）的个数为mn−1(a1,a2,...,an,m)

令d1=(a1,a2,...,an−1,m),d=(a1,a2,...,an−1,an,m),则d=(d1,an)且d1∣m,由同余性质，可得

a1x1+a2x2+...+an−1xn−1+anxn+b≡0(modd1)

又d1∣a1x1+a2x2+...+an−1xn−1,所以

anxn+b≡0(modd1)

设anxn+b=b1d1，代入原方程，得

a1x1+a2x2+...+an−1xn−1+b1d1≡0(modm)

原方程可依次消元并解出所有解。

---

解高次同余方程

仅讨论Ax≡B(modC)的情况。

当A与C互质时，由费马小定理，AC−1≡1(modC)且A0≡1(modC)则可得解的循环节小于C。因此，若选择暴力，则直接令x从0~C-1变化，检验其是否是方程的解。当C比较大时，则采用BSGS算法。

• BSGS(Baby Step Giant Step)算法

  • 先把x=i×m+j，其中m=ceil(C−−√)这样原式就变为Aim+j≡B(modC)再变为Aj≡B⋅A−mi(modC)
  • 先循环j=0→(m−1)，把(Aj%C,j)加入hash表中，这个就是Baby Steps
  • 然后我们再枚举等号右边B⋅A−mi(modC)，（这就是Giant Step）,从hash表中找看看有没有，有的话就得到了一组(i,j)，x=i*m+j，得到的这个就是正确解

//A^X=B(mod C)求X C为整数
#define MOD 76543
int hs[MOD],head[MOD],next[MOD],id[MOD],top;
void insert(int x, int y)
{
    int k = x%MOD;
    hs[top] = x, id[top] = y, next[top] = head[k], head[k] = top++;
}
int find(int x)
{
    int k = x%MOD;
    for(int i = head[k]; i != -1; i = next[i])
        if(hs[i] == x)
            return id[i];
    return -1;
}
int BSGS(int a,int b,int c)
{
    memset(head, -1, sizeof(head));
    top = 1;
    if(b == 1)
        return 0;
    int m = sqrt(c*1.0), j;
    long long x = 1, p = 1;
    for(int i = 0; i < m; ++i, p = p*a%c)
        insert(p*b%c, i);//存的是(a^j*b, j)
    for(long long i = m; ;i += m)
    {
        if( (j = find(x = x*p%c)) != -1 )
            return i-j;  //a^(ms-j)=b(mod c)
        if(i > c)
            break;
    }
    return -1;
}
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35

当A与C不互质时，先进行消因子处理再进行BSGS算法求解。

• 扩展BSGS算法

  • 原式等价为Ax+Cy=B，消因子处理，使得方程化为aAz+C′y=B′满足C′与A不再有可以约的因子,cnt=x−z。每次消因子过程中，方程右边也同时消去一样的因子。若不能，则说明方程不可解。将Az作为一个整体，利用扩展欧几里得算法得到其解系。假设出来的一组原始特解为x0，则Az=x0B′+KC′,K为任意整数。

  • 化为高次同余方程Az≡x0B′(modC′)利用BSGS求解得到z0，加回cnt即为原方程的解。

  • 注意：这样得到的只有不大于cnt的解。可能会漏掉小于cnt的解。因此先从1~log(MaxN)(一般设定为50即可)暴力查找一遍。

LL xiao(LL &a, LL &b, LL &c)
{
    LL r = 0, d, x, y, a1 = 1;
    while((d = extend_Euclid(a, c, x, y)) != 1)
    {
        if(b % d)
            return -1;
        c /= d;
        b /= d;
        a1 = a1*(a/d)%c;
        r++;
    }
    if(r == 0)
        return 0;
    extend_Euclid(a1, c, x, y);
    b = (b*x%c+c)%c;
    return r;
}
LL extend_BSGS(LL a,LL b,LL c)
{
    a %= c;   //简化运算
    LL r = 1;
    for(int i = 0;i < 50;i++)
    {
        if((r-b) % c == 0)
            return i;
        r *= a; r %= c;
    }
    memset(head, -1, sizeof(head));
    LL cnt = xiao(a, b, c);
    if(cnt == -1)
        return -1;
    top = 1;
    if(b == 1)
        return cnt;
    LL m = ceil(sqrt(c*1.0)), j;
    LL x = 1, p = 1;
    for(LL i = 0; i < m; ++i, p = p*a%c)
        insert(p*b%c, i);
    for(LL i = m; ;i += m)
    {
        if( (j = find(x = x*p%c)) != -1 )
            return i-j+cnt;
        if(i > c)
            break;
    }
    return -1;
}
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48

---

几类不定方程

毕达哥拉斯三元组

满足方程x2+y2=z2的三元组(x,y,z)为毕达哥拉斯三元组。当gcd(x,y,z)=1时，称其为本原的。

Thm x,y,z构成本原毕达哥拉斯三元组(y为偶数)⇔∃互素的正整数m,n(m>n分别为一奇一偶)且满足

⎧⎩⎨x=m2−n2y=2mnz=m2+n2

很容易证明后者的充分性。在《什么是数学》3书中证明了其必要性。从后者可以看出，本原毕达哥拉斯三元组中的最大数一定是奇数。 
特别地，有fnfn+3, 2fn+1fn+2, f2n+1+f2n+2构成毕达哥拉斯三元组。将fn=fn+2−fn+1, fn+3=fn+2+fn+1即得。 
如POJ1305

//求解n以内本原的毕达哥拉斯三元组有多少个
int Bida(int n)
{
    int m = floor(sqrt(n)+1e-6), r = 0;
    for(int i = 2;i < m;i++)
    {
        for(int j = 1;j < i;j++)
        {
            if((i-j) % 2==0 || gcd(i,j) != 1)//i,j要互质
                continue;
            int x = i*i-j*j, y = 2*i*j, z = i*i+j*j;
            if(z > n)
                continue;
            r++;
        }
    }
    return r;
}
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18

---

佩尔方程

形如整数方程

x2−dy2=1

其中d大于1且不为完全平方数。（当d为完全平方数时显然无解）

迭代公式 
Def 设r，s为整数，并且满足r2−Ns2=T，则称a=r−sN−−√给出该方程的解。

Thm 设a=x1−y1N−−√,b=x2−y2N−−√给出方程x2−Ny2=T的解，则ab=A−BN−−√给出方程x2−Ny2=T2的解。其中

{A=x1x2+Ny1y2B=x1y2+x2y1

取T=1，若最小特解为(x1,y1),则有迭代公式 

{xn=xn−1x1+dyn−1y1yn=xn−1y1+yn−1x1

可写为矩阵

(xnyn)=(x1y1dy1x1)n−1(x1y1)

利用快速幂也可快速求解。

求最小特解

• 暴力枚举
• 连分数法（待补坑）

\正整数方程x^2-dy^2=1,已知最小的特解x1,y1,求之后的每组解
#define N 100000
int X[N]={x1}, Y[N]={y1}, cnt;
int peir(int x1, int y1, int d, int n)
{
    cnt = 1;
    int x2 = x1, y2 = y1, x3, y3;
    for(int i = 1;i <= 10;i++)
    {
        x3 = x2*x1+d*y2*y1;
        y3 = x2*y1+y2*x1;
        x2 = x3, y2 = y3;
        X[cnt] = x3; y[cnt++] = y3;
    }
    return 0;
}
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16

\正整数方程x^2-dy^2=1,已知最小的特解x1,y1,求第n个解
void peir(int x1, int y1, int d, int n)
{
    int a[2][2] = {x1, d*y1, y1, x1};
    int b[2][2] = {x1, 0, y1, 0};
    A.set(a, 2, 2);
    B.set(b, 2, 1);
    C = po(A, n-1)*B;
    printf("%d %d
", C.a[0][0], C.a[1][0]);
}
1
2
3
4
5
6
7
8
9
10
1
2
3
4
5
6
7
8
9
10

---

同余式定理

费马小定理

Thm：若素数p满足p∤a，则ap−1≡1(modp)

利用构造法可以证明定理成立。

比较常用的做法是把除法转换成乘法

ap−2≡1a(modp)⇒ba≡b⋅ap−2(modp)

可由快速幂得到结果。为了调用方便，可以预处理数组，之后直接调用。 
可以试试做HDU 4869。相似的还有特别的斐波那契数列 HDU 4549 
当斐波那契数列的递推公式由加号改为乘号时，设

F(0)=a,F(1)=b,F(n+2)=F(n+1)×F(n)

则a,b的指数分别是前两项为{1,0} {0,1}的斐波那契数列。要求第n项模一个大数（给定的是一个质数）的值。当指数很大时，可以用费马小定理降低指数。

---

威尔逊定理

Thm 若p是素数，则(p−1)!≡−1(modp)

主要是证明在2,3,...,p−2中间存在两两的逆元对(x1,y1)(x1≠y1)，使得x1y1≡1(mod1),最后剩余的是1,p−1，则显然阶乘模p得到p−1

欧拉定理

Thm 设m是一个正整数，a是一个整数，(a,m)=1，则aϕ(x)≡1(modm)

特别的

若(a,b)=1,a,b为正整数,则aϕ(b)+bϕ(a)≡1(modab)

---

欧拉函数

ϕ(x)表示从1到x中与x互质的数的个数。 

n=pa11pa22...pann

ϕ(n)=n(1−1p1)(1−1p2)...(1−1pn)

void genPhi()
{
    for(int i = 1;i < max;i++)
        minDiv[i] = i;
    for(int i = 2;i*i < max;i++)
    {
        if(minDiv[i] == i)
        {
            for(int j = i*i;j < max;j += i)
                minDiv[j] = i;
        }
    }
    phi[1] = 1;
    for(int i = 2;i < max;i++)
    {
        phi[i] = phi[i/minDiv[i]];
        if((i/minDiv[i]) % minDiv[i] == 0)
            phi[i] *= minDiv[i];
        else
            phi[i] *= minDiv[i] - 1;
    }
}
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22

---

1. Mp=2p−1为第p个梅森数。若Mp为素数，则p必为素数，反之不成立。 ↩
2. 对于一个非素数n，任意一个小于n并与它互质的数a，an−1≡n(mod1)都成立 ↩
3. 【美】R.科朗 著 第一章补充第三节 ↩