Williams p+1整数分解算法实现、参数选择及测试问题咨询

我需要向同班同学介绍用于整数分解的Williams' p+1算法，但目前对该算法的掌握仍有不足。

算法前提与实现步骤疑问 #

据我理解，该算法的适用前提为待分解整数N可分解为素数乘积N=pq，且p+1是B-光滑的，我已经完成该前提下算法成立的证明，但不清楚如何正确实现和使用该算法，我设想的实现步骤如下：

• 在区间[1,N-1]内随机选取a
• 计算x = gcd(a,N)，若x≠1则返回x。我不理解为什么这一步不需要先校验x是否为素数：我们实际并不能确定N真的是两个素数的乘积，x有可能是合数，对吗？
• 正常情况下x=1，此时需要计算y = gcd(V_M-2,N)，其中V序列的递推规则为V₀ = 2, V₁ = a, Vₙ= aVₙ₋₁ - Vₙ₋₂。我已经找到了用模N矩阵幂计算Vₙ的方法，但不知道应该如何选择M的取值：我直接复用了Pollard算法的M计算逻辑，但不确定是否可行，也不明白背后的原理。
• 若y≠1且y≠N则返回y，此处我同样认为应该校验y是否为素数，我的想法是否正确？如果y不符合要求，就重新选取随机a重试。

核心疑问 #

我目前在实现上的核心疑问是M的构造逻辑，猜测和p+1的B-光滑性相关。在使用层面，我也不清楚该算法的适用场景和B参数的选取规则。

现有实现代码 #

下方是我编写的Python3实现代码：

import random
from math import floor, log, gcd

def is_prime(n): # 判断输入数是否为素数
    for d in range(2,n):
        if n%d == 0:
            return False
    return True

def primes_leq(B): # 获取所有小于等于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):
    # 幂为0时返回单位矩阵
    if power <= 0:
      return [[1,0],[0,1]]

    powers =  list(reversed([True if i=="1" else False for i in bin(power)[2:]])) # 按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("当前处理至M计算完成")
            
            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

测试用例疑问 #

我尝试用示例测试实现是否正确，调用williams(2**439-1,10**5)得到的因子是104110607，而参考记录中对应的因子应为122551752733003055543。据我所知这两个数都是N=2**439-1的素因子，这是不是和“N是两个素数的乘积”的前提矛盾？

问题解答 #

1. 关于因子是否需要校验素数的疑问
   Williams' p+1算法是通用合数分解算法，你学习证明时用到的“N=pq两素数乘积”只是简化场景，实际算法可以对任意合数输出非平凡因子，不管该因子是素数还是合数。如果你的目标是得到N的完整素因子分解，返回因子后额外做素数校验、对合数因子递归分解是合理的；如果只需要获取任意非平凡因子，不需要额外校验素数。

2. 关于M的构造逻辑
   你的猜测是正确的，M的构造和p+1的B-光滑性直接相关，复用Pollard p-1算法的M构造逻辑完全正确。M的本质是所有≤B的素数的足够高次幂的乘积，保证所有B-光滑数都能整除M。当p+1是B-光滑数时，p+1|M，结合Lucas序列的性质，此时V_M ≡ 2 mod p，因此gcd(V_M-2, N)就能将p这个因子提取出来，这正是算法的核心原理。

3. 关于适用场景与B参数选择

• 适用场景：Williams' p+1是Pollard p-1算法的补充，当Pollard p-1无法找到因子（即待分解素因子p的p-1有大素因子、不光滑）时，就可以尝试p+1算法，只要p+1是光滑的就能快速分解。
• B参数选择：通常从较小的值开始逐步提升，比如先试10³、再试10⁴、10⁵，B越大能覆盖的p+1素因子上限越高，但算法运行时间也越长，根据你的分解需求权衡即可。

4. 关于测试用例的疑问
   你的测试结果完全正常，不存在实现错误。2⁴³⁹-1本身就包含多个素因子，你得到的104110607的p+1光滑度更高，在当前B=1e5的参数下先被提取出来，和你学习的两素数简化证明场景并不冲突。

内容的提问来源于stack exchange，提问作者Nuria Gallego Ariño