Williams p+1 整数分解

Question

我被要求向我的同学介绍 Williams 的 p+1 算法来分解整数，但我认为我没有做对。据我了解，该算法采用整数 N 将其分解为素数 p,q (N=pq)，其中 p+1 是 B 平滑的。我理解为什么从这些前提开始，算法有效（我已经写过证明），但我不知道如何正确实现和使用它。我认为必须按如下方式实现：

1. 我在区间 [1,N-1] 中随机取一个
2. 我计算 x=gcd(a,N)。如果 x !=1，那么我返回 x （我不明白为什么我们不首先检查 x 是否为素数，因为我们实际上不知道 N 是否真的等于 p*q 并且 x 可以组合，对?)
3. 通常，x == 1，所以我必须计算 y = gcd(V_M-2,N) 其中 V_0 = 2, V_1 = a, V_n= aV_(n-1) - V_(n-2) 。我找到了一种计算 V_n 的方法来计算矩阵幂模数 N，但我不知道应该使用哪个 M（我复制了 Pollard 的方法，但我不知道这是否有效以及为什么）。
4. 如果 y!=1 且 y!=N，我返回 y（同样，与 x 一样，我认为我们应该检查 y 是否为素数，对吗？）。否则，只需尝试另一个随机 a 并重新开始。

所以，这主要是我关于实现的问题，总体上是关于 M 的构建，我猜这可能与 p+1 B-smoothness 的事实有关。

关于用法，我真的不知道在什么情况下我应该使用这种方法，我应该采取哪个B。我将在此处留下 Python3 中的代码以及一个真正让我抓狂的案例示例，看看您能否帮帮我。

import random
from math import floor, log, gcd

def is_prime(n): #funcion que determina si un numero es primo
    for d in range(2,n):
        if n%d == 0:
            return False
    return True

def primes_leq(B): #funcion para obtener los primos que son menor o igual que B
    l=[]
    for i in range(2,B+1):
        if is_prime(i):
            l.append(i)
    return l

def matrix_square(A, mod):
    return mat_mult(A,A,mod)


def mat_mult(A,B, mod):
  if mod is not None:
    return [[(A[0][0]*B[0][0] + A[0][1]*B[1][0])%mod, (A[0][0]*B[0][1] + A[0][1]*B[1][1])%mod],
            [(A[1][0]*B[0][0] + A[1][1]*B[1][0])%mod, (A[1][0]*B[0][1] + A[1][1]*B[1][1])%mod]]


def matrix_pow(M, power, mod):
    #Special definition for power=0:
    if power <= 0:
      return [[1,0],[0,1]]

    powers =  list(reversed([True if i=="1" else False for i in bin(power)[2:]])) #Order is 1,2,4,8,16,...

    matrices = [None for _ in powers]
    matrices[0] = M

    for i in range(1,len(powers)):
        matrices[i] = matrix_square(matrices[i-1], mod)

    result = None

    for matrix, power in zip(matrices, powers):
        if power:
            if result is None:
                result = matrix
            else:
                result = mat_mult(result, matrix, mod)

    return result

def williams(N, B):
    flag = False
    while not flag :
        a = random.randint(1,N-1)
        print("a : " + str(a))
        x = gcd(a,N)
        print("x : " + str(x))
        if x != 1:
            return x
        else :
            M = 1
            A = [[0,1],[-1,a]]
            
            for p in primes_leq(B):
                M *= p **(floor(log(N,p)))
            print("voy por aquí")
            
            C = matrix_pow(A,M,N)
            V = 2*C[0][0]+ a*C[0][1]
            y = gcd(V-2,N)
            print("y : " + str(y))
            if y != 1 and y != N:
                flag = True
    return y

为了测试我的实现，我尝试按照一些示例来检查我的因式分解是否工作正常。例如，我查看了https://members.loria.fr/PZimmermann/records/Pplus1.html 并尝试了williams(2**439-1,10**5)，我得到了 104110607，但我知道我应该得到 122551752733003055543（如网页中所示）。据我了解，两者都是因数 N=2**439-1 的素数，但这不是与 N 是两个素数 p*q 的乘积的假设相矛盾吗？

感谢您的帮助，我们将不胜感激

Answer 1

我认为你错过了这个算法的重点......

当 M 是 p+1 的倍数时，您会找到 N 的因数 p（或微除数）。

如果 p+1 是平滑的 -- 假设 p+1 的所有因素都是 M 是 all 这样可能的因素的倍数，像这样：

M=1
for x in all primes <= B:
    let y = largest power of x such that y < N
    M = M*y

您应该检查此循环产生的 M 的连续值。或者，您可以只检查所有连续的阶乘。关键是，在每次迭代中，你给M添加了新的因子，当p+1的所有因子都在M中时，那么M当然会是p+1的倍数。

棘手的部分是M会变得很大，而且你不能把M mod N。但是，您可以做的是计算所有 V_M mod N，而不是实际计算每个 M，您只需使用加法公式将 V_M 的下标乘以适当的因子：V_a+b = V_aV_b - V_a-b

Answer 2

2⁴³⁹-1 有两个以上的质因数。如果你没有得到 你想要的，你应该除以你得到的并保留 继续商数。