使用 Python 高效地找到原根模 n?
Efficient finding primitive roots modulo n using Python?
我正在使用以下代码在 Python 中查找 primitive roots 模 n
:
代码:
def gcd(a,b):
while b != 0:
a, b = b, a % b
return a
def primRoots(modulo):
roots = []
required_set = set(num for num in range (1, modulo) if gcd(num, modulo) == 1)
for g in range(1, modulo):
actual_set = set(pow(g, powers) % modulo for powers in range (1, modulo))
if required_set == actual_set:
roots.append(g)
return roots
if __name__ == "__main__":
p = 17
primitive_roots = primRoots(p)
print(primitive_roots)
输出:
[3, 5, 6, 7, 10, 11, 12, 14]
代码片段摘自: Diffie-Hellman (Github)
primRoots
方法能否在内存使用和性能/效率方面进行简化或优化?
您可以在此处进行的一个快速更改(尚未有效优化)是使用列表和集合理解:
def primRoots(modulo):
coprime_set = {num for num in range(1, modulo) if gcd(num, modulo) == 1}
return [g for g in range(1, modulo) if coprime_set == {pow(g, powers, modulo)
for powers in range(1, modulo)}]
现在,您可以在此处进行的一项强大而有趣的算法更改是优化您的 gcd
function using memoization。或者更好的是,您可以简单地使用 Python-3.5+ 中的内置 gcd
函数形式 math
模块或以前版本中的 fractions
模块:
from functools import wraps
def cache_gcd(f):
cache = {}
@wraps(f)
def wrapped(a, b):
key = (a, b)
try:
result = cache[key]
except KeyError:
result = cache[key] = f(a, b)
return result
return wrapped
@cache_gcd
def gcd(a,b):
while b != 0:
a, b = b, a % b
return a
# or just do the following (recommended)
# from math import gcd
然后:
def primRoots(modulo):
coprime_set = {num for num in range(1, modulo) if gcd(num, modulo) == 1}
return [g for g in range(1, modulo) if coprime_set == {pow(g, powers, modulo)
for powers in range(1, modulo)}]
如评论中所述,作为一种更 pythoinc 优化器的方式,您可以使用 fractions.gcd
(或 Python-3.5+ math.gcd
)。
根据 Pete 的评论和 Kasramvd 的回答,我可以这样建议:
from math import gcd as bltin_gcd
def primRoots(modulo):
required_set = {num for num in range(1, modulo) if bltin_gcd(num, modulo) }
return [g for g in range(1, modulo) if required_set == {pow(g, powers, modulo)
for powers in range(1, modulo)}]
print(primRoots(17))
输出:
[3, 5, 6, 7, 10, 11, 12, 14]
变化:
有关内置 gcd 的其他信息位于此处:
在 p 是素数的特殊情况下,以下速度要快一些:
import sys
# translated to Python from http://www.bluetulip.org/2014/programs/primitive.js
# (some rights may remain with the author of the above javascript code)
def isNotPrime(possible):
# We only test this here to protect people who copy and paste
# the code without reading the first sentence of the answer.
# In an application where you know the numbers are prime you
# will remove this function (and the call). If you need to
# test for primality, look for a more efficient algorithm, see
# for example Joseph F's answer on this page.
i = 2
while i*i <= possible:
if (possible % i) == 0:
return True
i = i + 1
return False
def primRoots(theNum):
if isNotPrime(theNum):
raise ValueError("Sorry, the number must be prime.")
o = 1
roots = []
r = 2
while r < theNum:
k = pow(r, o, theNum)
while (k > 1):
o = o + 1
k = (k * r) % theNum
if o == (theNum - 1):
roots.append(r)
o = 1
r = r + 1
return roots
print(primRoots(int(sys.argv[1])))
您可以通过使用更高效的算法来大大改进您的 isNotPrime 函数。您可以通过对偶数进行特殊测试然后仅测试奇数直到平方根来使速度加倍,但与 Miller Rabin 测试等算法相比,这仍然非常低效。 Rosetta Code 站点中的这个版本将始终为任何少于 25 位数字的数字给出正确答案。对于大素数,这将 运行 花费的时间只是使用试除法的一小部分。
此外,在处理本例中的整数时,您应该避免使用浮点求幂运算符 **(即使我刚刚链接到的 Rosetta 代码做同样的事情!)。在特定情况下,事情可能会正常工作,但当 Python 必须从浮点数转换为整数时,或者当整数太大而无法精确表示为浮点数时,它可能是一个微妙的错误来源。您可以改用有效的整数平方根算法。这是一个简单的:
def int_sqrt(n):
if n == 0:
return 0
x = n
y = (x + n//x)//2
while (y<x):
x=y
y = (x + n//x)//2
return x
这些代码在很多方面都是低效的,首先你不需要迭代所有 n 的互素提醒,你只需要检查欧拉函数从 n 的除法幂.在 n 是质数的情况下,欧拉的函数是 n-1。如果 n 是素数,则需要对 n-1 进行因式分解并仅检查那些除法器,而不是全部。这背后有一个简单的数学。
其次。你需要更好的函数来为数字提供动力想象一下功率太大,我认为在 python 中你有函数 pow(g, powers, modulo) 它在每个步骤中进行除法并仅获得余数( _ % modulo ).
如果您要实现 Diff & Helman 算法,最好使用安全素数。它们是这样的质数,即p是质数,2p+1也是质数,所以2p+1被称为安全质数。如果你得到 n = 2*p+1,那么 n-1 的除数(n 是素数,n 的欧拉函数是 n-1)是 1、2、p 和 2p,你只需要检查数字是否2 次方的 g 和 p 次方的 g 如果其中一个给出 1,则该 g 不是原根,您可以丢弃该 g 并 select 另一个 g,下一个 g+1,如果 g^ 2 和 g^p 不等于 1 的模 n,那么 g 是原根,检查保证除 2p 之外的所有幂都会给出不同于 1 的模 n.
示例代码使用Sophie Germain素数p和对应的安全素数2p+1,计算该安全素数2p+1的原根。
您可以轻松地为任何素数或任何其他数字重新编写代码,方法是添加一个函数来计算欧拉函数并找到该值的所有约数。但这只是一个演示,不是完整的代码。或许还有更好的方法。
class SGPrime :
'''
This object expects a Sophie Germain prime p, it does not check that it accept that as input.
Euler function from any prime is n-1, and the order (see method get_order) of any co-prime
remainder of n could be only a divider of Euler function value.
'''
def __init__(self, pSophieGermain ):
self.n = 2*pSophieGermain+1
#TODO! check if pSophieGermain is prime
#TODO! check if n is also prime.
#They both have to be primes, elsewhere the code does not work!
# Euler's function is n-1, #TODO for any n, calculate Euler's function from n
self.elrfunc = self.n-1
# All divisors of Euler's function value, #TODO for any n, get all divisors of the Euler's function value.
self.elrfunc_divisors = [1, 2, pSophieGermain, self.elrfunc]
def get_order(self, r):
'''
Calculate the order of a number, the minimal power at which r would be congruent with 1 by modulo p.
'''
r = r % self.n
for d in self.elrfunc_divisors:
if ( pow( r, d, self.n) == 1 ):
return d
return 0 # no such order, not possible if n is prime, - see small Fermat's theorem
def is_primitive_root(self, r):
'''
Check if r is a primitive root by modulo p. Such always exists if p is prime.
'''
return ( self.get_order(r) == self.elrfunc )
def find_all_primitive_roots(self, max_num_of_roots = None):
'''
Find all primitive roots, only for demo if n is large the list is large for DH or any other such algorithm
better to stop at first primitive roots.
'''
primitive_roots = []
for g in range(1, self.n):
if ( self.is_primitive_root(g) ):
primitive_roots.append(g)
if (( max_num_of_roots != None ) and (len(primitive_roots) >= max_num_of_roots)):
break
return primitive_roots
#demo, Sophie Germain's prime
p = 20963
sggen = SGPrime(p)
print (f"Safe prime : {sggen.n}, and primitive roots of {sggen.n} are : " )
print(sggen.find_all_primitive_roots())
此致
我正在使用以下代码在 Python 中查找 primitive roots 模 n
:
代码:
def gcd(a,b):
while b != 0:
a, b = b, a % b
return a
def primRoots(modulo):
roots = []
required_set = set(num for num in range (1, modulo) if gcd(num, modulo) == 1)
for g in range(1, modulo):
actual_set = set(pow(g, powers) % modulo for powers in range (1, modulo))
if required_set == actual_set:
roots.append(g)
return roots
if __name__ == "__main__":
p = 17
primitive_roots = primRoots(p)
print(primitive_roots)
输出:
[3, 5, 6, 7, 10, 11, 12, 14]
代码片段摘自: Diffie-Hellman (Github)
primRoots
方法能否在内存使用和性能/效率方面进行简化或优化?
您可以在此处进行的一个快速更改(尚未有效优化)是使用列表和集合理解:
def primRoots(modulo):
coprime_set = {num for num in range(1, modulo) if gcd(num, modulo) == 1}
return [g for g in range(1, modulo) if coprime_set == {pow(g, powers, modulo)
for powers in range(1, modulo)}]
现在,您可以在此处进行的一项强大而有趣的算法更改是优化您的 gcd
function using memoization。或者更好的是,您可以简单地使用 Python-3.5+ 中的内置 gcd
函数形式 math
模块或以前版本中的 fractions
模块:
from functools import wraps
def cache_gcd(f):
cache = {}
@wraps(f)
def wrapped(a, b):
key = (a, b)
try:
result = cache[key]
except KeyError:
result = cache[key] = f(a, b)
return result
return wrapped
@cache_gcd
def gcd(a,b):
while b != 0:
a, b = b, a % b
return a
# or just do the following (recommended)
# from math import gcd
然后:
def primRoots(modulo):
coprime_set = {num for num in range(1, modulo) if gcd(num, modulo) == 1}
return [g for g in range(1, modulo) if coprime_set == {pow(g, powers, modulo)
for powers in range(1, modulo)}]
如评论中所述,作为一种更 pythoinc 优化器的方式,您可以使用 fractions.gcd
(或 Python-3.5+ math.gcd
)。
根据 Pete 的评论和 Kasramvd 的回答,我可以这样建议:
from math import gcd as bltin_gcd
def primRoots(modulo):
required_set = {num for num in range(1, modulo) if bltin_gcd(num, modulo) }
return [g for g in range(1, modulo) if required_set == {pow(g, powers, modulo)
for powers in range(1, modulo)}]
print(primRoots(17))
输出:
[3, 5, 6, 7, 10, 11, 12, 14]
变化:
有关内置 gcd 的其他信息位于此处:
在 p 是素数的特殊情况下,以下速度要快一些:
import sys
# translated to Python from http://www.bluetulip.org/2014/programs/primitive.js
# (some rights may remain with the author of the above javascript code)
def isNotPrime(possible):
# We only test this here to protect people who copy and paste
# the code without reading the first sentence of the answer.
# In an application where you know the numbers are prime you
# will remove this function (and the call). If you need to
# test for primality, look for a more efficient algorithm, see
# for example Joseph F's answer on this page.
i = 2
while i*i <= possible:
if (possible % i) == 0:
return True
i = i + 1
return False
def primRoots(theNum):
if isNotPrime(theNum):
raise ValueError("Sorry, the number must be prime.")
o = 1
roots = []
r = 2
while r < theNum:
k = pow(r, o, theNum)
while (k > 1):
o = o + 1
k = (k * r) % theNum
if o == (theNum - 1):
roots.append(r)
o = 1
r = r + 1
return roots
print(primRoots(int(sys.argv[1])))
您可以通过使用更高效的算法来大大改进您的 isNotPrime 函数。您可以通过对偶数进行特殊测试然后仅测试奇数直到平方根来使速度加倍,但与 Miller Rabin 测试等算法相比,这仍然非常低效。 Rosetta Code 站点中的这个版本将始终为任何少于 25 位数字的数字给出正确答案。对于大素数,这将 运行 花费的时间只是使用试除法的一小部分。
此外,在处理本例中的整数时,您应该避免使用浮点求幂运算符 **(即使我刚刚链接到的 Rosetta 代码做同样的事情!)。在特定情况下,事情可能会正常工作,但当 Python 必须从浮点数转换为整数时,或者当整数太大而无法精确表示为浮点数时,它可能是一个微妙的错误来源。您可以改用有效的整数平方根算法。这是一个简单的:
def int_sqrt(n):
if n == 0:
return 0
x = n
y = (x + n//x)//2
while (y<x):
x=y
y = (x + n//x)//2
return x
这些代码在很多方面都是低效的,首先你不需要迭代所有 n 的互素提醒,你只需要检查欧拉函数从 n 的除法幂.在 n 是质数的情况下,欧拉的函数是 n-1。如果 n 是素数,则需要对 n-1 进行因式分解并仅检查那些除法器,而不是全部。这背后有一个简单的数学。
其次。你需要更好的函数来为数字提供动力想象一下功率太大,我认为在 python 中你有函数 pow(g, powers, modulo) 它在每个步骤中进行除法并仅获得余数( _ % modulo ).
如果您要实现 Diff & Helman 算法,最好使用安全素数。它们是这样的质数,即p是质数,2p+1也是质数,所以2p+1被称为安全质数。如果你得到 n = 2*p+1,那么 n-1 的除数(n 是素数,n 的欧拉函数是 n-1)是 1、2、p 和 2p,你只需要检查数字是否2 次方的 g 和 p 次方的 g 如果其中一个给出 1,则该 g 不是原根,您可以丢弃该 g 并 select 另一个 g,下一个 g+1,如果 g^ 2 和 g^p 不等于 1 的模 n,那么 g 是原根,检查保证除 2p 之外的所有幂都会给出不同于 1 的模 n.
示例代码使用Sophie Germain素数p和对应的安全素数2p+1,计算该安全素数2p+1的原根。
您可以轻松地为任何素数或任何其他数字重新编写代码,方法是添加一个函数来计算欧拉函数并找到该值的所有约数。但这只是一个演示,不是完整的代码。或许还有更好的方法。
class SGPrime :
'''
This object expects a Sophie Germain prime p, it does not check that it accept that as input.
Euler function from any prime is n-1, and the order (see method get_order) of any co-prime
remainder of n could be only a divider of Euler function value.
'''
def __init__(self, pSophieGermain ):
self.n = 2*pSophieGermain+1
#TODO! check if pSophieGermain is prime
#TODO! check if n is also prime.
#They both have to be primes, elsewhere the code does not work!
# Euler's function is n-1, #TODO for any n, calculate Euler's function from n
self.elrfunc = self.n-1
# All divisors of Euler's function value, #TODO for any n, get all divisors of the Euler's function value.
self.elrfunc_divisors = [1, 2, pSophieGermain, self.elrfunc]
def get_order(self, r):
'''
Calculate the order of a number, the minimal power at which r would be congruent with 1 by modulo p.
'''
r = r % self.n
for d in self.elrfunc_divisors:
if ( pow( r, d, self.n) == 1 ):
return d
return 0 # no such order, not possible if n is prime, - see small Fermat's theorem
def is_primitive_root(self, r):
'''
Check if r is a primitive root by modulo p. Such always exists if p is prime.
'''
return ( self.get_order(r) == self.elrfunc )
def find_all_primitive_roots(self, max_num_of_roots = None):
'''
Find all primitive roots, only for demo if n is large the list is large for DH or any other such algorithm
better to stop at first primitive roots.
'''
primitive_roots = []
for g in range(1, self.n):
if ( self.is_primitive_root(g) ):
primitive_roots.append(g)
if (( max_num_of_roots != None ) and (len(primitive_roots) >= max_num_of_roots)):
break
return primitive_roots
#demo, Sophie Germain's prime
p = 20963
sggen = SGPrime(p)
print (f"Safe prime : {sggen.n}, and primitive roots of {sggen.n} are : " )
print(sggen.find_all_primitive_roots())
此致