使用GMP高效地分解大数

9
我需要获取所有大于1k位的数字的质因数。这些数字几乎是随机的,所以应该不难。我用C++和GMP库。如何高效地完成这个任务?
编辑:我想你们都误解了我的意思。我所说的质因数是指得到数字的所有质因数。很抱歉我的英语有问题,在我的语言中,prime和factor是一样的。
澄清(来自OP的其他帖子):
我需要一种使用C++和GMP(Gnu Multiple Precession lib)或者其他方式高效地分解(找出数字的质因数)大数字(可能达到2048位)的方法。这些数字几乎是随机的,所以即使数字很难分解,我也可以重新生成数字(但不能选择)。

4
在这个语境中,"prime"是什么意思?你是想生成大素数吗?还是事先准备数字以便特定用途? - DGH
2
你按下数字旁边的小软按钮,将其提高几个档位,使其变为质数,引擎开始平稳运转。 - Potatoswatter
我同意,这篇文章的英语糟糕得让我能够猜到原帖作者的意思,但肯定不够清晰。 - Omnifarious
1
在我看来,具有两个质因数的合成数的困难情况的概率是不可忽略的(我的快速估计表明,在具有1024位的随机数中,大约有1/100000的概率)。 - starblue
@starblue:有趣!你能详细解释一下你的数学吗? - Jason S
显示剩余3条评论
5个回答

11
一个好的开始是使用小质数进行预过滤,比如所有小于100000的质数。尝试将其除以每一个质数(创建一个表格,然后在运行时加载它或将其作为静态数据放入您的代码中)。这可能看起来很慢和愚蠢,但如果数字是完全随机的,这将快速地给出一些因子,并且有着极高的概率。然后查看剩余的数字并决定下一步该怎么做。如果它相当小(“小”是由你自己定义的),您可以尝试进行素性测试(我认为GMP中有这样的东西),如果它是质数,那么在大多数情况下,您可以信任它。否则,您需要进一步分解它。
如果您的数字真的很大,并且您关心性能,那么您绝对需要实现比愚蠢的除法更复杂的算法。看看二次筛法(尝试维基百科)。它非常简单但非常强大。如果您准备好挑战,请尝试MPQS,它是二次筛选算法的变种。这个论坛是一个很好的信息来源。甚至已经存在您所需工具的现有实现 - 例如这个
请注意,1000位的数字在所有方面都是巨大的。分解这样一个数字(即使使用MPQS或其他方法)可能需要数年时间,如果您有运气的话,否则将永远无法完成。我认为,如果它们由两个几乎相等的质数组成(这当然是最困难的情况),那么MPQS对大约100-400位的数字表现良好。

对于超过约100位数字的数字,二次筛法开始被广义数域筛法(GNFS)所击败。 - Chris Card
1
@Chris Card - 一百位数字是在哪个进制下?请具体说明。 - Omnifarious
当然是以十进制为基础。对于这么大的数字,Chris是正确的,QS(或更明确地说,MPQS,因为QS在这样的数字上要慢得多)开始变得比GNFS慢。 - PeterK
1
当前记录的因数分解是RSA768(http://www.crypto-world.com/FactorRecords.html),这需要大量的工作。使用标准现代PC上的GNFS可以对大约150位数字进行因数分解,但是1k位目前甚至连顶尖研究人员也无法达到。 - Chris Card
是的,对于分解巨大数字的任何尝试,大多数情况下都是在分布式计算机系统上完成的,其中每个节点都为结果做出了一小部分的贡献。即使如此,完成任务也需要数月甚至数年的时间。 - PeterK
你不需要逐个除以小质数来筛选它们。只需将所有小质数的乘积(这部分是预处理)取出,然后计算最大公约数即可。 - RghtHndSd

6
以下是Java中的一个示例算法(不是使用GMP的C++,但转换应该非常简单),它实现了以下功能:
  1. 生成一个长度为Nbits的随机数x
  2. 尝试分解所有小于100的质因数,并保留能够整除x的质因数列表。
  3. 使用Java的isProbablePrime方法测试剩余的因子是否为质数
  4. 如果剩下的因子是质数且概率足够高,我们已经成功地将x分解。(停止
  5. 否则,剩下的因子一定是合数(请参阅isProbablePrime文档)。
  6. 在还有时间的情况下,运行Pollard rho算法直到找到一个除数d。
  7. 如果我们用完了时间,我们就失败了。(停止
  8. 我们找到了一个除数d。因此,我们分解出d,将d的质因数添加到x的质因数列表中,并转到步骤4。
此算法的所有参数都在程序清单的开头附近。我寻找1024位的随机数,超时时间为250毫秒,并且一直运行程序,直到我获得至少有4个质因数的数字x(有时程序会先找到1、2或3个质因数的数字)。使用这个参数设置,我的2.66Ghz iMac通常需要大约15-20秒的时间。
Pollard's rho算法并不是非常高效,但与二次筛选(QS)或广义数域筛选(GNFS)相比,它很简单 - 我只想看看这个简单算法如何工作。

为什么这个方法可行:(尽管你们中的许多人认为这是一个难题)

事实是,质数并不是那么稀有。对于1024位数,素数定理表明大约每1024 ln 2(约710)个数字中就有一个是素数。

因此,如果我生成一个随机数x,它是素数,并且我接受概率素数检测,那么我就成功地分解了x。

如果它不是素数,但我快速地分解出一些小因子,并且剩余因子是(概率上)素数,那么我就成功地分解了x。

否则,我会放弃并生成一个新的随机数。(这是OP认为可接受的)

大多数成功分解的数字将具有1个大素数因子和几个小素数因子。

难以分解的数字是那些没有小素数因子且至少有2个大素数因子的数字(其中包括两个大数的积作为加密密钥;OP没有提到加密),当我时间不够时,我可以跳过它们。

package com.example;

import java.math.BigInteger;
import java.util.ArrayList;
import java.util.List;
import java.util.Random;

public class FindLargeRandomComposite {
    final static private int[] smallPrimes = {2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 
        31, 37, 41, 43, 47, 53, 59, 61, 67, 71, 
        73, 79, 83, 89, 97};

    final static private int maxTime = 250;
    final static private int Nbits = 1024;
    final static private int minFactors = 4;
    final static private int NCERTAINTY = 4096;

    private interface Predicate { public boolean isTrue(); }

    static public void main(String[] args)
    {
        Random r = new Random();
        boolean found = false;
        BigInteger x=null;
        List<BigInteger> factors=null;
        long startTime = System.currentTimeMillis();
        while (!found)
        {
            x = new BigInteger(Nbits, r);
            factors = new ArrayList<BigInteger>();
            Predicate keepRunning = new Predicate() {
                final private long stopTime = System.currentTimeMillis() + maxTime;
                public boolean isTrue() {
                    return System.currentTimeMillis() < stopTime;
                }
            };
            found = factor(x, factors, keepRunning);
            System.out.println((found?(factors.size()+" factors "):"not factored ")+x+"= product: "+factors);
            if (factors.size() < minFactors)
                found = false;
        }
        long stopTime = System.currentTimeMillis();
        System.out.println("Product verification: "+(x.equals(product(factors))?"passed":"failed"));
        System.out.println("elapsed time: "+(stopTime-startTime)+" msec");
    }

    private static BigInteger product(List<BigInteger> factors) {
        BigInteger result = BigInteger.ONE;
        for (BigInteger f : factors)
            result = result.multiply(f);
        return result;
    }

    private static BigInteger findFactor(BigInteger x, List<BigInteger> factors,
            BigInteger divisor)
    {
        BigInteger[] qr = x.divideAndRemainder(divisor);
        if (qr[1].equals(BigInteger.ZERO))
        {
            factors.add(divisor);
            return qr[0];
        }
        else
            return x;
    }

    private static BigInteger findRepeatedFactor(BigInteger x,
            List<BigInteger> factors, BigInteger p) {
        BigInteger xprev = null;
        while (xprev != x)
        {
            xprev = x;
            x = findFactor(x, factors, p);
        }
        return x;
    }

    private static BigInteger f(BigInteger x, BigInteger n)
    {
        return x.multiply(x).add(BigInteger.ONE).mod(n);
    }

    private static BigInteger gcd(BigInteger a, BigInteger b) {
        while (!b.equals(BigInteger.ZERO))
        {
            BigInteger nextb = a.mod(b);
            a = b;
            b = nextb;
        }
        return a;
    }
    private static BigInteger tryPollardRho(BigInteger n,
            List<BigInteger> factors, Predicate keepRunning) {
        BigInteger x = new BigInteger("2");
        BigInteger y = x;
        BigInteger d = BigInteger.ONE;
        while (d.equals(BigInteger.ONE) && keepRunning.isTrue())
        {
            x = f(x,n);
            y = f(f(y,n),n);
            d = gcd(x.subtract(y).abs(), n);
        }
        if (d.equals(n))
            return x;
        BigInteger[] qr = n.divideAndRemainder(d);
        if (!qr[1].equals(BigInteger.ZERO))
            throw new IllegalStateException("Huh?");
        // d is a factor of x. But it may not be prime, so run it through the factoring algorithm.
        factor(d, factors, keepRunning);
        return qr[0];
    }

    private static boolean factor(BigInteger x0, List<BigInteger> factors,
            Predicate keepRunning) {

        BigInteger x = x0;
        for (int p0 : smallPrimes)
        {
            BigInteger p = new BigInteger(Integer.toString(p0));
            x = findRepeatedFactor(x, factors, p);          
        }
        boolean done = false;
        while (!done && keepRunning.isTrue())
        {
            done = x.equals(BigInteger.ONE) || x.isProbablePrime(NCERTAINTY);
            if (!done)
            {
                x = tryPollardRho(x, factors, keepRunning);
            }
        }
        if (!x.equals(BigInteger.ONE))
            factors.add(x);
        return done;
    }
}

4

如果你要分解的数字具有小质因子,可以使用Pollard p-1分解算法。它已经将数字2 ^ 740 + 1的30位质因子分解出来。ECM是一种类似但次指数的算法,但实现更加困难。算法所需时间取决于设定的边界b。它将分解任何具有因子p且p - 1为b平滑的数字。

//Pollard p - 1 factorization algorithm

void factor(mpz_t g, mpz_t n, long b)
{
    //sieve for primes
    std::vector<bool> r;

    for(int i = 0; i < b; i++)
        r.push_back(true);


    for(int i = 2; i < ceil(sqrt(b - 1)); i++)
        if(r.at(i) == true)
            for(int j = i * i; j < b; j += i)
                r.at(j) = false;

    std::vector<long> p;
    std::vector<long> a;
    for(int i = 2; i < b; i++)
        if(r[i] == true)
        {
            p.push_back(i);//Append the prime on to the vector
            int temp = floor(log(b) / log(i)); //temp = logb(i)

            // put primes in to sieve
            // a = the maximum power for p ^ a < bound b
            if(temp == 0)
                a.push_back(1);
            else
                a.push_back(temp);                
        }

    int m = p.size();//m = number of primes under bound b

    mpz_t c;// c is the number Which will be exponated
    mpz_init(c);
    long two = 2;
    mpz_set_ui(c, two);// set c to 2

    int z = 0;
    long x = 2;

    // loop c until a factor is found
    for(;;)
    {
    mpz_set_si( c, x);

    //powering ladder
    for(long i = 0; i < m; i++)
        for(long j = 0; j < a[i]; j++)
            mpz_powm_ui(c , c, (p[i]), n);

    //check if a factor has been found;
    mpz_sub_ui(c ,c,1);
    mpz_gcd(g ,c, n);
    mpz_add_ui(c , c, 1);

    //if g is a factor return else increment c
    if((mpz_cmp_si(g,1)) > 0 && (mpz_cmp(g,n)) < 0)
        return;
    else if (x > b)
        break;
    else
        x++;
    }

}


int main()
{
    mpz_t x;
    mpz_t g;

    //intialize g and x
    mpz_init(g);
    mpz_init_set_str(x,"167698757698757868925234234253423534235342655234234235342353423546435347",10);

    //p-1 will factor x as long as it has a factor p where p - 1 is b-smooth(has all prime factors less than bound b)
    factor(g , x, 1000);

    //output the factor, it will output 1 if algorithm fails
    mpz_out_str(NULL, 10, g);

    return 0;
}

输出 - 7465647 执行时间 - 0.003 秒

J.Pollard 创造的另一个分解算法是 Pollards Rho 算法,它不是很快但需要非常少的空间。还有并行化的方法。它的复杂度为 O(n^1/4)。

//Pollard rho factoring algorithm
void rho(mpz_t g, mpz_t n)
{
    mpz_t x;
    mpz_t y;
    mpz_init_set_ui(x ,2);
    mpz_init_set_ui(y ,2);//initialize x and y as 2
    mpz_set_ui(g , 1);
    mpz_t temp;
    mpz_init(temp);

    if(mpz_probab_prime_p(n,25) != 0)
        return;//test if n is prime with miller rabin test

    int count;
    int t1 = 0;
    int t2 = 1;
    int nextTerm = t1 + t2;
    while(mpz_cmp_ui(g,1) < 1)
    {
        f(x,n);//x is changed
        f(y,n);//y is going through the sequence twice as fast
        f(y,n);

        if(count == nextTerm)//calculate gcd every fibonacci number
        {
            mpz_sub(temp,x,y);
            mpz_gcd(g , temp, n);

            t1 = t2;
            t2 = nextTerm;
            nextTerm = t1 + t2;//calculate next fibonacci number
        }

        count ++;
    }

    return;
}

int main()
{
    mpz_t x;
    mpz_t g;

    //intialize g and x
    mpz_init(g);
    mpz_init_set_str(x,"167698757698757868925234234253423",10);


    rho(g , x);

    //output the factor, it will output 1 if algorithm fails
    mpz_out_str(NULL, 10, g);

    return 0;
}

输出 - 353 执行时间 - 0.003秒


“f(x,n);”中的函数 f() 在哪里? - raspiduino

2

2

你提供的任务很有趣!谢谢!

对我来说,用了两天时间编写高级解决方案是一件非常愉快的事情。我从头开始实现了三个分解算法:试除法Pollard RhoLenstra椭圆曲线方法(ECM)

众所周知,ECM 方法(基于椭圆曲线)是中等范围因素分解的最快方法之一。虽然试除法适合处理非常小的因子(每秒可达到2^20个因子),而Pollard Rho则适合更大但仍然比较小的因子(每秒可达到2^40个因子),而 ECM 则适合于中等范围因子(每10秒可达到2^60个因子)。

还有一些非常先进的方法,如通用数域筛(GNFS)(每月可实现2^700的因子分解),但它们很难实现。还有二次筛法方法也很先进(可能每月可实现2^400的因子分解),我也从头开始实现了这个算法,但代码太长了,虽然可以理解,但由于其体积,我没有在此附上。ECM 方法是高级方法中唯一相对容易实现的方法。

除了我实现的上述3种分解方法之外,我还在代码中使用了以下算法:模指数运算费马素性测试Barrett减少算法欧几里得算法扩展欧几里得算法模乘逆算法椭圆曲线点加法和乘法

实际上,您的任务非常容易快速解决:对于2048位大小的比特,在大约每个ln(2048) = 1420次数中会出现一个质数,因此您只需快速生成大约1500个数字,同时检查它们是否为素数,例如使用非常快的费马素性测试。如果一个数字是素数,那么根据定义它已经被分解了。
我在我的思维中进一步扩展了您的任务,使其更有趣。我不寻找质数,而是试图找到这样的2048位随机生成的数字,使其至少具有多个大的质因数。我称之为“有趣”的数字。当然,如果一个数字有几个微小的因子和一个大的质数,那么它就不那么有趣。但如果它有60位质因数,那么捕捉这样的数字就很有趣,这就是我在代码中所做的。
您可以在我的代码中看到,我将其适用于两种库:Boost MultiprecisionGMP。两者都包含在我的代码中(请参见#include <boost/multiprecision/cpp_int.hpp>#include <gmpxx.h>),因此您应该安装并链接两者。在Linux下,通过sudo apt install libboost-all-dev libgmp-dev非常容易安装两者。但是在Windows上有点棘手,请先安装Chocolatey Packet Manager,然后执行命令choco install boost-msvc-14.3。对于GMP,请按这里描述的方式安装VCPKG,然后vcpkg install gmp。如果您想要的话,也可以通过VCPKG安装Boost:vcpkg install boost
ECM(椭圆曲线)方法非常有趣和简单:
  1. 您生成许多随机曲线,带有随机的X、Y、A、B、N参数,其中N是需要因式分解的输入数字,其他参数是适合曲线方程的随机数Y^2 = X^3 + A * x + B (mod N)
  2. 将每个曲线乘以所有增长的质数(带一些小幂)。
  3. 在某个点上,您将获得第一个曲线的曲线阶的倍数,并且由于这种情况下曲线阶的属性,您将获得所谓的曲线上的无限点。
  4. 如果您查看点加法公式,则可以看到公式中有一个分母lambda = (y_q - y_p) / (x_q - x_p)。该分母是模N的模乘法逆元计算得出的。对于无限点,它变成不可逆,而根据逆数不可逆的属性,只有当GCD(N,x_q-x_p)!=1时才可能不可逆,此时GCD会给出一些非平凡因子(有时也是N),因此我们成功地通过给出第一个因子来分解了N。
  5. 如果我们没有得到无限点,则继续生成更多的曲线,并除以更多(越来越大)的质数。我们生成的曲线越多,乘以的质数越多,因式分解的成功率就越高。

在线尝试!

源代码在此处。由于StackOverflow每篇文章的符号限制为30K,而我的代码本身约为44K字节,因此我无法将其内联在此处,而是通过Github Gist共享它(下面链接)。同样的代码也可以通过Godbolt服务器上的在线尝试!链接获得。

GitHub Gist 代码

示例控制台输出:

TrialDiv time 8.230 sec
Num: 4343663370925180057127849780941698665126534031938076094687921681578209757551374613160773985436765755919464255163981465381983273353052491 (2^453.90)
Factored: 13 (2^3.70, prime), 167 (2^7.38, prime), 3853 (2^11.91, prime), 53831 (2^15.72, prime), 916471 (2^19.81, prime), 9255383 (2^23.14, prime), 
UnFactored: 11372390351822722497588418782148940973499109818654526670537593527638523385195910987808859992169924704037636069779 (2^372.24, composite), 
PollardRho time 8.51 sec
Num: 11372390351822722497588418782148940973499109818654526670537593527638523385195910987808859992169924704037636069779 (2^372.24)
Factored: 189379811 (2^27.50, prime), 2315962907 (2^31.11, prime), 50213994043 (2^35.55, prime), 
UnFactored: 5163708449171395447719565208770850251589387410889704005960043195676697732937073689 (2^278.09, composite), 
Curves   1, Ops   12.965 Ki, ExtraPrimes 512, Primes  0.500 Ki, IterTime   0.410 sec
Curves   2, Ops   50.912 Ki, ExtraPrimes 512, Primes  1.000 Ki, IterTime   8.062 sec
Curves   3, Ops  112.586 Ki, ExtraPrimes 464, Primes  1.453 Ki, IterTime  15.093 sec
Curves   4, Ops  162.931 Ki, ExtraPrimes 120, Primes  1.570 Ki, IterTime   4.853 sec
Curves   5, Ops  193.699 Ki, ExtraPrimes  80, Primes  1.648 Ki, IterTime   4.201 sec
ECM time 32.624 sec
Num: 5163708449171395447719565208770850251589387410889704005960043195676697732937073689 (2^278.09)
Factored: 955928964443 (2^39.80, prime), 
UnFactored: 540177004907359979630305557131905764121354872876442621652476639261690523 (2^238.29, composite), 
Final time 49.385 sec
Num: 4343663370925180057127849780941698665126534031938076094687921681578209757551374613160773985436765755919464255163981465381983273353052491 (2^453.90)
Factored: 13 (2^3.70, prime), 167 (2^7.38, prime), 3853 (2^11.91, prime), 53831 (2^15.72, prime), 916471 (2^19.81, prime), 9255383 (2^23.14, prime), 189379811 (2^27.50, prime), 2315962907 (2^31.11, prime), 50213994043 (2^35.55, prime), 955928964443 (2^39.80, prime), 
UnFactored: 540177004907359979630305557131905764121354872876442621652476639261690523 (2^238.29, composite), 

网页内容由stack overflow 提供, 点击上面的
可以查看英文原文,
原文链接