On 03/08/2014 01:12 PM, Mark D. Baushke wrote:
> You raise a good point. Folks with 'bc -l' should be able to verify the
> MODP groups.
The k-bit MODP discrete-log moduli is formulated as the smallest
safe-prime that meets the form:
2^k - 2^(k-64) - 1 + 2^64 * { [2^(k-130) pi] + N }
where N is a positive integer.
More colloquially, this is: 64 bits of 0xFF, followed by (k-128) bits of
pi, followed by 64 more bits of 0xFF. I don't know why this sequence
was selected. Does anyone have any pointers to reasons you might want
the modulus structured this way?
Given a bit-length k, The attached python program tests all integers N
in sequence, stopping at the first safe-prime. For example:
0 dkg@alice:~$ python3 ./gen-modp.py 768
searching for 768-bit safe-prime modulus, starting at:
FFFFFFFF FFFFFFFF C90FDAA2 2168C234 C4C6628B 80DC1CD1
29024E08 8A67CC74 020BBEA6 3B139B22 514A0879 8E3404DD
EF9519B3 CD3A431B 302B0A6D F25F1437 4FE1356D 6D51C245
E485B576 625E7EC6 F44C42E9 A637ED6B FFFFFFFF FFFFFFFF
added 0x248b6 * 2^64 (decimal: 149686 * 2^64)
The prime is: 2^768 - 2^704 - 1 + 2^64 * { [2^638 pi] + 149686 }
FFFFFFFF FFFFFFFF C90FDAA2 2168C234 C4C6628B 80DC1CD1
29024E08 8A67CC74 020BBEA6 3B139B22 514A0879 8E3404DD
EF9519B3 CD3A431B 302B0A6D F25F1437 4FE1356D 6D51C245
E485B576 625E7EC6 F44C42E9 A63A3620 FFFFFFFF FFFFFFFF
0 dkg@alice:~$
Thus far, I've used it to generate the moduli for these groups:
* 768-bit (MODP 1)
* 1024-bit (MODP 2)
* 1536-bit (MODP 5)
* 2048-bit (MODP 14)
* 3072-bit (MODP 15)
i haven't yet generated these (mainly due to time):
* 4096-bit (MODP 16)
* 6144-bit (MODP 17)
* 8192-bit (MODP 18)
Regards,
--dkg
#!/usr/bin/python
# Author: Daniel Kahn Gillmor <dkg%fifthhorseman.net@localhost>
# Date: 2014-03-11
# derive p for the default groups in RFC 3526 and RFC 5996
# these RFCs use a particular approach to derive the prime modulus,
# based only on the size of the prime desired, and on the digits of pi.
# this code performs the same search pattern and verifies the result.
# Usage:
# $ python3 gen-modp.py 1024
# this will generate the 1024-bit modulus found in RFC 5996.
import sys
import gmpy2
def breakup(s,k):
'''return string s broken into chunks of size k'''
return (s[0+n:k+n] for n in range(0,len(s),k))
def displayhex(mpz,stream):
'''write an mpz to stream in the form preferred by the MODP RFC's'''
h = hex(mpz)[2:].upper() # strip leading 0x
linelen = 48
if len(h) == 6144//4:
linelen += 8 # section 6 of RFC 3526 is unusually-formatted
for x in breakup(h,linelen):
for y in breakup(x,8):
stream.write(" " + y)
stream.write("\n")
def t_pi(n):
'''return an mpz containing the first n bits of pi'''
gmpy2.get_context().precision=n+1
return gmpy2.mpz(gmpy2.floor((gmpy2.mpz(1) << (n-2)) * gmpy2.const_pi(n+1)))
def checkanswer(bitsize,inc):
'''Given a number of bits and an selected offset, verify that the
expected result of the MODP modulus selection scheme is a safe prime.
We're using probabilistic primality testing here, but leaning hard on
the test, making sure it's expensive and much more thorough than the
default tests used during selection.
'''
one = gmpy2.mpz(1)
val = (one << bitsize) - (one << (bitsize-64)) - one + \
((one<<64)* (t_pi(bitsize - 128) + inc))
if not (gmpy2.is_prime(val, bitsize) and \
gmpy2.is_prime((val-1)//2,bitsize)):
raise Exception("expensive primality tests contradicted "+\
"the cheap ones!")
print("The prime is: 2^%d - 2^%d - 1 + 2^64 * { [2^%d pi] + %d }"%\
(bitsize, bitsize-64, bitsize-(64+66), inc))
displayhex(val, sys.stdout)
def search(bitsize):
'''search for a prime with bitsize bits, using the pattern established
by the MODP RFC's.
'''
bits = gmpy2.mpz(0xffffffffffffffff)
# 64 bits of 0xff
# the remaining bits of pi
# 64 bits of 0xff
# then count up until we find a prime
base = (bits << (bitsize-64)) + (t_pi(bitsize-128) << 64) + bits
print("searching for %d-bit safe-prime modulus, starting at:"%(bitsize,))
displayhex(base,sys.stdout)
inc = 1 # initial offset started to clear the middle bits for pi
while not gmpy2.is_prime(base) \
or not gmpy2.is_prime((base-1)>>1):
if (inc % 17) == 1:
sys.stdout.write("\radding 0x%x * 2^64"%(inc))
inc += 1
base += 2**64
sys.stdout.write("\radded 0x%x * 2^64 (decimal: %d * 2^64)\n"%\
(inc,inc))
return inc
if sys.argv.__len__() > 1:
target = int(sys.argv[1])
else:
target = 768
inc = search(target)
checkanswer(target,inc)
Attachment:
signature.asc
Description: OpenPGP digital signature