整数分解

程序设计

发布日期: 2018-09-03

1.试除法

试除法是整数分解算法中最简单和最容易理解的算法。首次出现于意大利数学家斐波那契出版于1202年的著作。有合数n，n为待分解的正整数，从小于等于 $\sqrt{n}$ 的每个素数去试除待分解的整数。如果找到一个数能够整除除尽，这个数就是待分解整数的因子。如果找到一个数能够整除除尽，这个数就是待分解整数的因子。试除法一定能够找到n的因子，如果没有找到则证明n是素数，如果整数的末尾不是0或5那么可以不用计算5的倍数的素因子，如果n末尾为2，检查偶数因子即可。

某种意义上说，试除法是个效率非常低的算法，如果从2开始，一直算到 $\sqrt{n}$ 需要 $\pi ({\sqrt {n}})$ )次试除，这里pi(x)是小于x的素数的个数。这是不包括素性测试的。如果稍做变通——还是不包括素性测试——用小于 $\sqrt{n}$ 的奇数去简单的试除，则需要 ${ {\sqrt {n}} \over 2}$ )次。这意味着，如果n有大小接近的素因子（例如公钥密码学中用到的），试除法是不太可能实行的。但是，当n有至少一个小因子，试除法可以很快找到这个小因子。值得注意的是，对于随机的n，2是其因子的概率是50％，3是33％，等等，88％的正整数有小于100的因子，91％的有小于1000。 1.1卢卡斯-莱默检测法 梅森数是指形如 $2^{n}-1$ )的数，记为 $M_{n}$ )；如果一个梅森数是素数那么它称为梅森素数（英语：Mersenne prime）。卢卡斯－莱默检验法原理是这样：令梅森数 M__p = 2_p_− 1作为检验对象（预设p是素数，否则M__p_就是合数了）。定义序列{_s__i }：所有的i ≥ 0

$s_{i}={\begin{cases}\;\;\,4\\s_{ {i-1}}^{2}-2\end{cases}}$

，如果 $\displaystyle i=0$ ；

，如果 $\displaystyle i>0$ 。

这个序列的开始几项是4, 14, 194, 37634, … （OEIS中的数列A003010）那么_M__p_是素数当且仅当

$s_{ {p-2}}\equiv 0{\pmod {M_{p}}};$

否则, M__p_是合数。_s__p − 2模_M__p_的余数叫做p的卢卡斯－莱默余数。 1.2 埃拉托斯特尼筛法 埃拉托斯特尼筛法（希腊语：κόσκινον Ἐρατοσθένους，英语：sieve of Eratosthenes ），简称埃氏筛，也有人称素数筛。这是一种简单且历史悠久的筛法，用来找出一定范围内所有的素数。所使用的原理是从2开始，将每个素数的各个倍数，标记成合数。一个素数的各个倍数，是一个差为此素数本身的等差数列。此为这个筛法和试除法不同的关键之处，后者是以素数来测试每个待测数能否被整除。埃拉托斯特尼筛法是列出所有小素数最有效的方法之一，其名字来自于古希腊数学家埃拉托斯特尼，并且被描述在另一位古希腊数学家尼科马库斯所著的《算术入门》中。[1] 详细列出算法如下：

列出2以后的所有序列：
- 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25
标出序列中的第一个质数，也就是2，序列变成：
- 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25
将剩下序列中，划摽2的倍数（用红色标出），序列变成：
- 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25
如果现在这个序列中最大数小于最后一个标出的素数的平方，那么剩下的序列中所有的数都是质数，否则回到第二步。

本例中，因为25大于2的平方，我们返回第二步：
剩下的序列中第一个质数是3，将主序列中3的倍数划出（红色），主序列变成：
- 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25
我们得到的质数有：2，3
25仍然大于3的平方，所以我们还要返回第二步：
现在序列中第一个质数是5，同样将序列中5的倍数划出，主序列成了：
- 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25
我们得到的质数有：2, 3, 5 。
因为25等于5的平方，结束循环

结论：去掉红色的数字，2到25之间的质数是：2, 3, 5, 7, 11, 13, 17, 19, 23。 1.2.1伪代码

Input: an integer n > 1

Let A be an array of Boolean values, indexed by integers 2 to n,
initially all set to true.

for i = 2, 3, 4, …, not exceeding √n:
if A[i] is true:
for j = i2, i2+i, i2+2i, i2+3i, …, not exceeding n :
A[j] := false

Output: all i such that A[i] is true.

1.3.1 C++

bool flag[MAXN] = {1};    //将标识初始化为true
void erat(int maxn)
{
    flag[0]=0;            //0不是素数
    flag[1]=0;            //1不是素数
    for(int i=2;i<=maxn;++i)
    {
        /*当i为素数时，i的所有倍数都不是素数*/
        if(flag[i])
        {
            for(int j=i*i;j<=maxn;j+=i)
            {
                flag[j]=0;
            }
        }
    }
 }


`#include <iostream>`

`#include <cstdio>`

`#include <cmath>`

`#include <cstring>`

`#include <vector>`

`#include <algorithm>`

`using` `namespace` `std;`

`const` `long` `long` `maxn = 10000007 + 10;`

`const` `long` `long` `maxp = 700000;`

`int` `vis[maxn];    ``// i是合数vis为1，i是素数，vis为0`

`long` `long` `prime[maxp];`

`void` `sieve(``long` `long` `n){`

`long` `long` `m = (``long` `long``)``sqrt``(n + 0.5);`

`memset``(vis, 0, ``sizeof``(vis));`

`vis[2] = 0;`

`for` `(``long` `long` `i = 3; i <= m; i += 2) {`

`if``(!vis[i])`

`for` `(``long` `long` `j = i * i; j <= n; j += i)`

`vis[j] = 1;`

`if``(i * i > n)`

`break``;`

`}`

`}`

`long` `long` `gen(``long` `long` `n){`

`sieve(n);`

`long` `long` `c = 1;`

`prime[0] = 2;`

`for` `(``long` `long` `i = 3; i <= n; i += 2)`

`if``(!vis[i])`

`prime[c++] = i;`

`return` `c;`

`}`

`int` `main()`

`{`

`long` `long` `n, c;`

`cout << ``"刷素数到n："``;`

`cin >> n;`

`c = gen(n);`

`for` `(``long` `long` `i = 0; i < c; i++)`

`printf``(``"%lld"``, prime[i]);`

`cout << endl << c;`

`return` `0;`

`}`

1.3 米勒-拉宾素性检验 米勒-拉宾素性检验是一种素数判定法则，利用随机化算法判断一个数是合数还是_可能是_素数。卡内基梅隆大学的计算机系教授Gary Lee Miller首先提出了基于广义黎曼猜想的确定性算法，由于广义黎曼猜想并没有被证明，其后由以色列耶路撒冷希伯来大学的Michael O. Rabin教授作出修改，提出了不依赖于该假设的随机化算法。 1.3.1 概念首先介绍一个相关的引理。 $1^{2}{\bmod {p}}$ 和 $(-1)^{2}{\bmod {p}}$ 总是得到 $1$ )，称这两个数为1的“平凡平方根”。当 $p$ 是素数且 $p>2$ 时，不存在 $1{\bmod {p}}$ 的“非平凡平方根”。为了证明该引理，我们假设 $x$ 是 $1{\bmod {p}}$ 的平方根，于是有

$x^{2}\equiv 1{\pmod {p}}$

$(x+1)(x-1)\equiv 0{\pmod {p}}$

也就是说，素数 $p$ 能够整除 $(x-1)(x+1)$ ，根据欧几里得引理， $x-1$ 或者 $x+1$ )能够被 $p$ 整除，即 $x\equiv 1{\pmod {p}}$ 或 $x\equiv -1{\pmod {p}}$ , 即 $x$ 是 $1{\bmod {p}}$ 的平凡平方根。现在假设 $n$ 是一个奇素数，且 $n>2$ )。于是 $n-1$ 是一个偶数，可以被表示为 $2^{s}*d$ 的形式， $s$ )和 $d$ )都是正整数且 $d$ 是奇数。对任意在 $(Z/nZ)^{*}$ 范围内的 $a$ 和 $0\leq r\leq s-1$ ，必满足以下两种形式的一种：

$a^{d}\equiv 1{\pmod {n}}$

$a^{2^{r}d}\equiv -1{\pmod {n}}$

因为由于费马小定理，对于一个素数 $n$ ，有

$a^{n-1}\equiv 1{\pmod {n}}$

又由于上面的引理，如果我们不断对 $a^{n-1}$ 取平方根，我们总会得到 1 或 -1。如果我们得到了 -1，意味着②式成立， $n$ 是一个素数。如果我们从未得到 -1，那么通过这个过程我们已经取遍了所有2的幂次，即①式成立。米勒-拉宾素性测试就是基于上述原理的逆否，也就是说，如果如果我们能找到这样一个 $a$ ，使得对任意 $0\leq r\leq s-1$ 以下两个式子均满足：

$a^{d}\not \equiv 1{\pmod {n}}$

$a^{2^{r}d}\not \equiv -1{\pmod {n}}$

那么 $n$ 就不是一个素数。这样的 $a$ 称为 $n$ 是合数的一个凭证（witness）。否则 $a$ 可能是是一个证明 $n$ 是素数的“强伪证”（strong liar），即当 $n$ )确实是一个合数，但是对当前选取的 $a$ 来说上述两个式子均不满足，这时我们认为 $n$ )是基于 $a$ 的大概率素数。每个奇合数 $n$ 都有很多个对应的凭证 $a$ ，但是要生成这些 $a$ 并不容易。当前解决的办法是使用概率性的测试。我们从集合 $(Z/nZ)^{*}$ )中随机选择非零数 $a$ ，然后检测当前的 $a$ 是否是 $n$ 为合数的一个凭证。如果 $n$ 本身确实是一个合数，那么大部分被选择的 $a$ 都应该是 $n$ 的凭证，也即通过这个测试能检测出 $n$ 是合数的可能性很大。然而，仍有极小概率的情况下我们找到的 $a$ 是一个“强伪证”（反而表明了 $n$ 可能是一个素数）。通过多次独立测试不同的 $a$ ，我们能减少这种出错的概率。对于测试一个大数是否是素数，常见的做法随机选取基数 $a$ ，毕竟我们并不知道凭证和伪证在 $[1,n-1]$ 这个区间如何分布。典型的例子是 Arnault 曾经给出了一个397位的合数 $n$ )，但是所有小于307的基数 $a$ )都是 $n$ 是素数的“强伪证”。不出所料，这个大合数通过了 Maple 程序的isprime() 函数（被判定为素数）。这个函数通过检测 $a=2,3,5,7,11$ 这几种情况来进行素性检验。 1.3.2 伪代码


Input #1: n > 3, an odd integer to be tested for primality;
Input #2: k, a parameter that determines the accuracy of the test
Output: composite if n is composite, otherwise probably prime

write n − 1 as 2r·d with d odd by factoring powers of 2 from n − 1
WitnessLoop: repeat k times:
    pick a random integer a in the range \[2, n − 2\]
    x ← ad mod n
    if x = 1 or x = n − 1 then
        continue WitnessLoop
    repeat r − 1 times:
    x ← x2 mod n
    if x = 1 then
        return composite
    if x = n − 1 then
        continue WitnessLoop
    return composite
return probably prime

1.3.3 c++


#include <iostream>
#include <cstdlib>
using namespace std;
typedef long long llong;
//求取(x * y) % n
llong mod(llong x, llong y,llong n)
{
    llong res = 0;
    llong temp = x % n;
    while(y)
   {
        if(y & 0x1)
            if((res += temp) > n)
                res -= n;
        if((temp <<= 1) \>  n)
            temp -= n;
        y >>= 1;
    }
    return res;
}
//求取(x ^ y) % n
llong get_mod(llong x, llong y, llong n)
{
    llong res = 1;
    llong temp = x;
    while(y)
    {
        if(y & 0x1)
            res = mod(res, temp, n);
        temp = mod(temp, temp, n);
        y >>= 1;
    }
    return res;
}
//编写bool函数，判定是否为素数
bool is_prime(llong n, int t)
{
    if(n < 2)
        return false;
    if(n == 2)
        return true;
    if(!(n & 0x1))
        return false;
    llong k = 0, m, a, i;
    for(m = n -1; !(m & 0x1); m >>= 1, ++k);
    while(t--)
    {
        a = get_mod(rand() % (n - 2) \+ 2, m, n);
        if(a != 1)
        {
            for(i = 0; i < k && a != n-1; ++i)
            {
                cout << a << endl;
                a = mod(a, a, n);
            }
            //根据二次探测定理，只要不满足(a == 1) || (a == n - 1)，就会一直遍历下去，直到最后返回false
            if(i >= k
                return false;
        }
    }
    return true;
}
//主函数
int main()
{
    int times;
    llong num;
    cin >\> times;
    while(times--)
    {
        cin >\> num;
        if(is_prime(num, 1))
            cout << "Yes" << endl;
        else
            cout << "No" << endl;
    }
    return 0;
}

2.轮式因子分解法


### 3.Pollard's rho algorithm

`/* C++ program to find a prime factor of composite using`

`Pollard's Rho algorithm */`

`#include<bits/stdc++.h>`

`using` `namespace` `std;`

`/* Function to calculate (base^exponent)%modulus */`

`long` `long` `int` `modular_pow(``long` `long` `int` `base,` `int` `exponent,`

`long` `long` `int` `modulus)`

`{`

`/* initialize result */`

`long` `long` `int` `result = 1;`

`while` `(exponent > 0)`

`{`

`/* if y is odd, multiply base with result */`

`if` `(exponent & 1)`

`result = (result * base) % modulus;`

`/* exponent = exponent/2 */`

`exponent = exponent >> 1;`

`/* base = base * base */`

`base = (base * base) % modulus;`

`}`

`return` `result;`

`}`

`/* method to return prime divisor for n */`

`long` `long` `int` `PollardRho(``long` `long` `int` `n)`

`{`

`/* initialize random seed */`

`srand` `(``time``(NULL));`

`/* no prime divisor for 1 */`

`if` `(n==1)` `return` `n;`

`/* even number means one of the divisors is 2 */`

`if` `(n % 2 == 0)` `return` `2;`

`/* we will pick from the range [2, N) */`

`long` `long` `int` `x = (``rand``()%(n-2))+2;`

`long` `long` `int` `y = x;`

`/* the constant in f(x).`

`* Algorithm can be re-run with a different c`

`* if it throws failure for a composite. */`

`long` `long` `int` `c = (``rand``()%(n-1))+1;`

`/* Initialize candidate divisor (or result) */`

`long` `long` `int` `d = 1; `

`/* until the prime factor isn't obtained.`

`If n is prime, return n */`

`while` `(d==1)`

`{`

`/* Tortoise Move: x(i+1) = f(x(i)) */`

`x = (modular_pow(x, 2, n) + c + n)%n;`

`/* Hare Move: y(i+1) = f(f(y(i))) */`

`y = (modular_pow(y, 2, n) + c + n)%n;`

`y = (modular_pow(y, 2, n) + c + n)%n;`

`/* check gcd of |x-y| and n */`

`d = __gcd(``abs``(x-y), n);`

`/* retry if the algorithm fails to find prime factor`

`* with chosen x and c */`

`if` `(d==n)` `return` `PollardRho(n);`

`}`

`return` `d;`

`}`

`/* driver function */`

`int` `main()`

`{`

`long` `long` `int` `n = 10967535067;`

`printf``(``"One of the divisors for %lld is %lld."``,`

`n, PollardRho(n));`

`return` `0;`

`}`

Jinzhengxu

http://jinzhengxu.github.io/2018/09/03/zheng-shu-fen-jie/