Williams p+1 整数分解

Question

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

1. 我在区间 [1,N-1] 内随机取 a
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 平滑度的事实有关。

关于用法，我实在不明白什么情况下该用这个方法，该用哪个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 的所有因子都是 <= B -- then it becomes becomes possible to construct an 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 具有两个以上素因数。如果你没有得到 你想要的，你应该除以你得到的并保留 继续商。

Answer 3

我们的同伴和朋友们以 p+1 的方式呈现威廉姆斯的方法，并以类似的方式进行了类似的操作，我们将继续采取行动，以解决数字问题或计算数字问题。没有任何人参与实施算法和数字的任务。 Llegaste a resolver tus dudas conrespecto a esto？ Nos podrías echar una mano？