Noc*_*wer 4 algorithm math primes
有人告诉我,Frobenius pseudoprime算法运行时间比Miller-Rabin素数测试长三倍,但分辨率是其七倍.那么,如果一个人运行前十次和后三十次运行,两者都需要相同的时间运行,但前者将提供大约233%的分析能力.在试图找出如何进行测试时,最后用算法发现了以下文章:
尝试实现下面的算法,但程序从不打印出数字.那些更熟悉数学符号或算法的人可以验证发生了什么事吗？
编辑1:下面的代码添加了更正,但compute_wm_wm1缺少实施.有人可以从算法的角度解释递归定义吗？它不是为我"点击".
编辑2:已删除错误代码,并在compute_wm_wm1下面添加了该功能的实现.它似乎有效,但可能需要进一步优化才能实用.
from random import SystemRandom
from fractions import gcd
random = SystemRandom().randrange
def find_prime_number(bits, test):
number = random((1 << bits - 1) + 1, 1 << bits, 2)
while True:
for _ in range(test):
if not frobenius_pseudoprime(number):
break
else:
return number
number += 2
def frobenius_pseudoprime(integer):
assert integer & 1 and integer >= 3
a, b, d = choose_ab(integer)
w1 = (a ** 2 * extended_gcd(b, integer)[0] - 2) % integer
m = (integer - jacobi_symbol(d, integer)) >> 1
wm, wm1 = compute_wm_wm1(w1, m, integer)
if w1 * wm != 2 * wm1 % integer:
return False
b = pow(b, (integer - 1) >> 1, integer)
return b * wm % integer == 2
def choose_ab(integer):
a, b = random(1, integer), random(1, integer)
d = a ** 2 - 4 * b
while is_square(d) or gcd(2 * d * a * b, integer) != 1:
a, b = random(1, integer), random(1, integer)
d = a ** 2 - 4 * b
return a, b, d
def is_square(integer):
if integer < 0:
return False
if integer < 2:
return True
x = integer >> 1
seen = set([x])
while x * x != integer:
x = (x + integer // x) >> 1
if x in seen:
return False
seen.add(x)
return True
def extended_gcd(n, d):
x1, x2, y1, y2 = 0, 1, 1, 0
while d:
n, (q, d) = d, divmod(n, d)
x1, x2, y1, y2 = x2 - q * x1, x1, y2 - q * y1, y1
return x2, y2
def jacobi_symbol(n, d):
j = 1
while n:
while not n & 1:
n >>= 1
if d & 7 in {3, 5}:
j = -j
n, d = d, n
if n & 3 == 3 == d & 3:
j = -j
n %= d
return j if d == 1 else 0
def compute_wm_wm1(w1, m, n):
a, b = 2, w1
for shift in range(m.bit_length() - 1, -1, -1):
if m >> shift & 1:
a, b = (a * b - w1) % n, (b * b - 2) % n
else:
a, b = (a * a - 2) % n, (a * b - w1) % n
return a, b
print('Probably prime:\n', find_prime_number(300, 10))
Dan*_*her 14
由于不熟悉符号,您似乎完全误解了算法.
def frobenius_pseudoprime(integer):
assert integer & 1 and integer >= 3
a, b, d = choose_ab(integer)
w1 = (a ** 2 // b - 2) % integer
这来自于这条线
w ^ 0 ≡2(模N)和W 1 ≡一个2 b -1 - 2(MOD N)
但B -1并不意味着1/b这里,但模逆的b模数n,即整数c带b·c ? 1 (mod n).通过扩展的欧几里德算法,您可以c通过连续的分数展开b/n或等效地但稍微更多的计算来找到这样的a .既然你可能不熟悉连续分数,我推荐后者.
m = (integer - d // integer) // 2
来自
n - (Δ/ n)= 2m
并且误解了Jacobi符号作为分数/除法(诚然,我在这里显示它更像是一个分数,但由于该站点不支持LaTeX渲染,我们必须做到).雅可比符号是勒让德符号的一般化-相同地标识-其指示号码是否是二次剩余模奇素数(如果n是一个二次剩余模p,即,存在一个k与k^2 ? n (mod p)和n不的倍数p,然后(n/p) = 1,如果n是p,然后(n/p) = 0,否则(n/p) = -1)的倍数.雅可比符号解除了"分母"为奇数素数并允许任意奇数作为"分母"的限制.它的值是勒让德符号与所有素数除以n(根据多重性)相同的"分子"的乘积.更多相关内容,以及如何在链接文章中有效地计算雅可比符号.该行应该正确读取
m = (integer - jacobi_symbol(d,integer)) // 2
以下几行我完全无法理解,从逻辑上讲,这里应该使用递归计算W m和W m + 1
W¯¯ 2J ≡W¯¯ Ĵ 2 - 2(MOD N)
w ^ 2J + 1 ≡W¯¯ Ĵ w ^ J + 1 - w ^ 1(MOD N)
在PDF的公式(11)周围给出了使用该递归来计算所需值的有效方法.
w_m0 = w1 * 2 // m % integer
w_m1 = w1 * 2 // (m + 1) % integer
w_m2 = (w_m0 * w_m1 - w1) % integer
函数的其余部分几乎是正确的,当然,由于早期的误解,它现在得到了错误的数据.
if w1 * w_m0 != 2 * w_m2:
这里的(in)相等应该是模数integer,即if (w1*w_m0 - 2*w_m2) % integer != 0.
return False
b = pow(b, (integer - 1) // 2, integer)
return b * w_m0 % integer == 2
但请注意,如果n是素数,那么
b^((n-1)/2) ? (b/n) (mod n)
其中(b/n)是勒让德(或雅可比)符号(对于素数'分母',雅可比符号是勒让德符号),因此b^((n-1)/2) ? ±1 (mod n).所以你可以使用它作为额外的检查,如果W m不是2或者n-2,n不能是素数,也不能b^((n-1)/2) (mod n)是1或者不是n-1.
可能b^((n-1)/2) (mod n)首先计算并检查是否为1或者n-1是一个好主意,因为如果检查失败(顺便说一下是Euler pseudoprime测试),你不再需要另一个,也不需要更便宜的计算,如果成功的话,无论如何,你很可能需要计算它.
关于更正,他们似乎是正确的,除了一个故障,我以前忽略可能更糟糕:
if w1 * wm != 2 * wm1 % integer:
这仅适用于模数2 * wm1.
关于W j的递归,我认为最好用一个有效的实现来解释,首先是toto,以便于复制和粘贴:
def compute_wm_wm1(w1,m,n):
a, b = 2, w1
bits = int(log(m,2)) - 2
if bits < 0:
bits = 0
mask = 1 << bits
while mask <= m:
mask <<= 1
mask >>= 1
while mask > 0:
if (mask & m) != 0:
a, b = (a*b-w1)%n, (b*b-2)%n
else:
a, b = (a*a-2)%n, (a*b-w1)%n
mask >>= 1
return a, b
然后解释:
def compute_wm_wm1(w1,m,n):
我们需要W 1的值,所需数字的索引以及将模数作为输入的数量.值W 0始终为2,因此我们不需要将其作为参数.
称之为
wm, wm1 = compute_wm_wm1(w1,m,integer)
in frobenius_pseudoprime(旁边:不是一个好名字,大部分返回的数字True都是真正的素数).
a, b = 2, w1
我们初始化a并b与W 0和W 1分别.在每个点处,a保持W的值Ĵ和bW的值J + 1,其中,j是的比特的值m到目前为止消耗.例如,使用m = 13,值j,a并b发展如下:
consumed remaining j a b
1101 0 w_0 w_1
1 101 1 w_1 w_2
11 01 3 w_3 w_4
110 1 6 w_6 w_7
1101 13 w_13 w_14
这些位是从左到右消耗的,因此我们必须找到第一个设置位m并将我们的"指针"放在它之前
bits = int(log(m,2)) - 2
if bits < 0:
bits = 0
mask = 1 << bits
我从计算的对数中减去一点只是为了完全确定我们不会被浮点错误所欺骗(顺便说一下,使用log限制你最多1024位的数字,大约308个十进制数字;如果你想对待更大的数字,你必须以m不同的方式找到base-2对数,使用log是最简单的方法,它只是一个概念证明,所以我在这里使用它).
while mask <= m:
mask <<= 1
移动掩码直到它大于m,所以设置位指向m第一个设置位之前.然后向后移动一个位置,所以我们指向该位.
mask >>= 1
while mask > 0:
if (mask & m) != 0:
a, b = (a*b-w1)%n, (b*b-2)%n
如果设置了下一位,则消耗位的初始部分的值m从0 j开始2*j+1,因此我们需要的W序列的下一个值是W 2j + 1 for a和W 2j + 2 for b.通过上面的递归公式,
W_{2j+1} = W_j * W_{j+1} - W_1 (mod n)
W_{2j+2} = W_{j+1}^2 - 2 (mod n)
由于a式为W Ĵ和b式为W j + 1,a变为(a*b - W_1) % n和b变(b * b - 2) % n.
else:
a, b = (a*a-2)%n, (a*b-w1)%n
如果未设置下一位,所消耗的比特的初始部分的值m从进入j到2*j,所以a变得W¯¯ 2J =(W Ĵ 2 - 2)(模n),并且b变成w ^ 2J + 1 =(W Ĵ*W j + 1 - W 1)(mod n).
mask >>= 1
将指针移动到下一位.当我们移过最后一位时,mask变为0并且循环结束.消耗位的初始部分m现在是所有m的位,因此值当然是m.然后我们可以
return a, b
一些补充说明:
def find_prime_number(bits, test):
while True:
number = random(3, 1 << bits, 2)
for _ in range(test):
if not frobenius_pseudoprime(number):
break
else:
return number
大数字中的素数并不太频繁,因此只需选择随机数就可能需要大量尝试.如果您选择一个随机数并按顺序检查候选人,您可能会更快地找到素数(或可能的素数).
另一点是,如Frobenius测试这样的测试不成比例地发现例如3的倍数是复合的.在使用这样的测试(或Miller-Rabin测试,或Lucas测试,或Euler测试,......)之前,你应该做一些试验分工来清除大部分复合材料,并且只在那里工作它有战斗机会值得.
哦,is_square函数不准备处理小于2的参数,潜在的零除错误,
def is_square(integer):
if integer < 0:
return False
if integer < 2:
return True
x = integer // 2
应该有所帮助