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(。我已经找到了一种通过矩阵幂模n计算V_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+1s的所有因子都是<=B——那么就有可能构建一个M，它是所有这些可能因素的倍数，比如：

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会变得很大，并且你不能接受MmodN。不过，你可以做的是计算所有V_MmodN，而不是实际计算每个M。你只需使用加法公式将V_M的下标乘以适当的因子：V-_a+b=V_a-b

Answer 2

2⁴³⁹−1有两个以上的素因子。如果你没有你想要的，你应该除以你得到的，然后保留进行商运算。