Subject: Re: I can has SHA-1 hashes for RFC 2409/3526 MODP groups?

From: Peter Gutmann <[email-masked]>
Date: Fri, 7 Mar 2014 12:09:14 +0000

Geoffrey Keating <[email-masked]> writes:

>I'd encourage you to do the derivation again: compute
>
>2^2048 - 2^1984 - 1 + 2^64 * { [2^1918 pi] + 124476 }
>
>and verify that it's prime.

That assumes that (a) my calculation of that (what on earth is "{ [2^1918 pi]
+ 124476 }", for example?) will be correct, and (b) my overall calculations
will also be correct, which is more or less the thing I'm trying to avoid: I'd
like an independent check on the values so that if I've messed up anywhere, I
can detect it.  Getting a hash of the byte string seems to be the easiest way
to do this.

Speaking of doing the derivation again, do we know if anyone's actually tried
to reproduce the values given in the RFC?  I'm assuming it came from something
like Mathematica which I don't have directly available, and Mathics gives me,
for '2^2048 - 2^1984 - 1 + 2^64 * ((2^1918 * pi) + 124476)' the not terribly
helpful:

32317006071311007298962968279132820423223593024363796434440041316671543565873

which if fed to bc as:

32317006071311007298962968279132820423223593024363796434440041316671543565873
440916351244591771371847086149164662580591029693953927426555290878070964232927
695481445976665844136502242593753417365574610202650424305468123288740527986069
345182600707793335155055488870009358517077286832777252562865057729531709228943
077561263607903104272439999956067085492978564938102680241320640398099496994465
906209396510474094085543397979368739086794572799533369234200999216018730602022
705374745273152831923492968543993770580612325275982255795197317429444934474851
545487458964741061916562950289834595688065052813922041488166690320547839+18446
744073709551616*(124476+237427862951972853513674659239712424976731176586712791
857002390826462915716740084687764468520802762540432372120969949748452696941370
146798185382596204798325192468656471080494445867437239094388362141523332248179
455193652520455298404233189464786144777417873508682027722667738642546640317659
478006398477769954880726243534178790472636787575530378644324814366518352670268
626390091728133613800395438606915682767011151451496901064489828517647053407532
449143217533161990782555667852270573728212149993563349925018485149102146341244
69728944541006176128457062647724487151147320327110918144*pi)

is reported to be:

323170060713110072989629682791328204232235930243637964344400413166715435658734
409163512445917713718470861491646625805910296939539274265552908780709642329276
954814459766658441365022425937534173655746102026504243054681232887405279860693
451826007077933351550554888700093585170772868327772525628650577295317092289430
775612636079031042724399999560670854929785649381026802413206403980994969944659
062093965104740940855433979793687390867945727995333692342009992160187306020227
053747452731528319234929685439937705806123252759822557951973174294449344748515
45487458964741061916562950289834595688065052816218218403485760467501055

which is:

0xffffffffffffffff000000000000000000000000000000000000000000000000000000000000
000000000000000000000000000000000000000000000000000000000000000000000000000000
000000000000000000000000000000000000000000000000000000000000000000000000000000
000000000000000000000000000000000000000000000000000000000000000000000000000000
000000000000000000000000000000000000000000000000000000000000000000000000000000
000000000000000000000000000000000000000000000000000000000000000000000000000000
00000000000000000000000001e63bffffffffffffffff

which isn't right.

Peter.

---

From: Mouse <[email-masked]>
Date: Fri, 7 Mar 2014 18:18:11 -0500 (EST)

>> and verify that it's prime.

Or, at least, that it's probably prime.

> That assumes that (a) my calculation of that (what on earth is "{
> [2^1918 pi] + 124476 }", for example?)

Take pi, multiply by the 1918th power of two, discard the fractional
part, and add 124476.  Or, at least, that's what I'd read it as.

This means computing pi to nearly two thousand bits of precision (or
finding a reference value you trust enough), but that's really the only
hard part.

> will be correct, and (b) my overall calculations will also be
> correct, which is more or less the thing I'm trying to avoid: I'd
> like an independent check on the values so that if I've messed up
> anywhere, I can detect it.

Well, if you get the same thing the RFC specifies, that's reasonable
confirmation you've done the calculation correctly.

> Getting a hash of the byte string seems to be the easiest way to do
> this.

Not just computing it and comparing it against the value in RFC3526?

> Speaking of doing the derivation again, do we know if anyone's
> actually tried to reproduce the values given in the RFC?  I'm
> assuming it came from something like Mathematica which I don't have
> directly available, and Mathics gives me, for '2^2048 - 2^1984 - 1 +
> 2^64 * ((2^1918 * pi) + 124476)' the not terribly helpful:

Few general-purpose calculator programs will actually compute pi to the
necessary number of places.  I know of someone who tried such a thing
and got a wrong result which turned out to match the value resulting
from using the irrational number in question (pi, sqrt(2), whatever it
was) truncated to IEEE double precision.

> [...]
> which if fed to bc as:
> [...]+[...]*(124476+[...]*pi)
> is reported to be:
> [...]
> which is:
> 0xffffffffffffffff000000000000000000000000000000000000000
> 000

> which isn't right.

However, if you add 18 more 0s to the long string of them (which the
formatting makes it look likely got dropped in whatever produced the
weird line breaks above - if I add 18 more 0s and paste the first two
lines together, the result is exactly as long as the next five lines),
that _is_ 2^2048 - 2^1984 - 1 + (2^64 * 124476).  That is, it's what
you get if you take pi to be 0.  I suspect your bc didn't recognize the
character string "pi" and treated it as an unset variable or some such,
effectively replacing it with 0.

I have a calculator program capable of working with numbers of this
size.  I told it to compute pi in hex to 500 places; truncating the
result to 1918 fraction bits gives (still in hex, linebreaks added for
email consumption)

3.243f6a8885a308d313198a2e03707344a4093822299f31d0082efa98ec4e6c
89452821e638d01377be5466cf34e90c6cc0ac29b7c97c50dd3f84d5b5b54709
179216d5d98979fb1bd1310ba698dfb5ac2ffd72dbd01adfb7b8e1afed6a267e
96ba7c9045f12c7f9924a19947b3916cf70801f2e2858efc16636920d871574e
69a458fea3f4933d7e0d95748f728eb658718bcd5882154aee7b54a41dc25a59
b59c30d5392af26013c5d1b023286085f0ca417918b8db38ef8e79dcb0603a18
0e6c9e0e8bb01e8a3ed71577c1bd314b2778af2fda55605c60e65525f3aa55ab
945748986263e8144055ca396a2aab10b4

Setting pi to this value and telling it to compute (still working in
hex, hence the #d prefixes all over and the cvt() call because it
doesn't like mixed-base arithmetic)

pow(2,#d2048)-pow(2,#d1984)-1+(pow(2,#d64)*((pi*pow(2,#d1918))+cvt(#d124476,#d16)))

I get (again, line breaks added for email)

which is the value in RFC3526 (after case conversion and whitespace
fixup, the two text strings are identical).

I haven't done anything about checking it for primality.

If anyone would like to look at the calculator program in question, git
clone git://git.rodents-montreal.org/Mouse/calc and look at calc.c.  It
might even build for you, though not with the Makefile from the git
repo unless you pick up
ftp.rodents-montreal.org:/mouse/local/src/makefiles/makefiles-20100906/local-prog
and maybe not even then.  (icalc.c is ignorable - known broken.)

/~\ The ASCII				  Mouse
\ / Ribbon Campaign
 X  Against HTML		[email-masked]
/ \ Email!	     7D C8 61 52 5D E7 2D 39  4E F1 31 3E E8 B3 27 4B

---

From: Kurt Roeckx <[email-masked]>
Date: Sat, 8 Mar 2014 11:21:42 +0100

> effectively replacing it with 0.

If you want pi in bc, you need to start bc with the -l option and do:
pi=4*a(1)

Kurt

---

From: "Mark D. Baushke" <[email-masked]>
Date: Sat, 8 Mar 2014 10:12:58 -0800

Hi Kurt,

You raise a good point. Folks with 'bc -l' should be able to verify the
MODP groups.

> If you want pi in bc, you need to start bc with the -l option and do:

Yes. Of course, to get the same values as is used in the RFCs more
fractional bits need to be calculated

Using bc to calculate the group14 and group15 primes:

$ bc -l
scale=1000 /* should be enough for group 14 and group 15 */
obase=16
/* A function to return the largest integral value not greater than x */
define floor(x) {
  auto s
  s = scale
  scale = 0
  x /= 1  /* round x down */
  scale = s
  return (x)
}

pi=4*a(1)
pi
3.243F6A8885A308D313198A2E03707344A4093822299F31D0082EFA98EC4E6C8945\
2821E638D01377BE5466CF34E90C6CC0AC29B7C97C50DD3F84D5B5B54709179216D5\
D98979FB1BD1310BA698DFB5AC2FFD72DBD01ADFB7B8E1AFED6A267E96BA7C9045F1\
2C7F9924A19947B3916CF70801F2E2858EFC16636920D871574E69A458FEA3F4933D\
7E0D95748F728EB658718BCD5882154AEE7B54A41DC25A59B59C30D5392AF26013C5\
D1B023286085F0CA417918B8DB38EF8E79DCB0603A180E6C9E0E8BB01E8A3ED71577\
C1BD314B2778AF2FDA55605C60E65525F3AA55AB945748986263E8144055CA396A2A\
AB10B6B4CC5C341141E8CEA15486AF7C72E993B3EE1411636FBC2A2BA9C55D741831\
F6CE5C3E169B87931EAFD6BA336C24CF5C7A325381289586773B8F48986B4BB9AFC4\
BFE81B6628219361D809CCFB21A991487CAC605DEC8032EF845D5DE98575B1DC2623\
02EB651B8823893E81D396ACC50F6D6FF383F442392E0B4482A484200469C8F04A9E\
1F9B5E21C66842F6E96C9A670C9C61ABD388F06A51A0D2D8542F68960FA728AB5133\
A36EEF0B6C137A3B8
2^2048 - 2^1984 - 1 + 2^64 * (floor(pi * 2^1918) + 124476) /* group 14 */
FFFFFFFFFFFFFFFFC90FDAA22168C234C4C6628B80DC1CD129024E088A67CC74020B\
BEA63B139B22514A08798E3404DDEF9519B3CD3A431B302B0A6DF25F14374FE1356D\
6D51C245E485B576625E7EC6F44C42E9A637ED6B0BFF5CB6F406B7EDEE386BFB5A89\
9FA5AE9F24117C4B1FE649286651ECE45B3DC2007CB8A163BF0598DA48361C55D39A\
69163FA8FD24CF5F83655D23DCA3AD961C62F356208552BB9ED529077096966D670C\
354E4ABC9804F1746C08CA18217C32905E462E36CE3BE39E772C180E86039B2783A2\
EC07A28FB5C55DF06F4C52C9DE2BCBF6955817183995497CEA956AE515D2261898FA\
051015728E5A8AACAA68FFFFFFFFFFFFFFFF
2^3072 - 2^3008 - 1 + 2^64 * (floor(2^2942 * pi) + 1690314) /* group 15 */
FFFFFFFFFFFFFFFFC90FDAA22168C234C4C6628B80DC1CD129024E088A67CC74020B\
BEA63B139B22514A08798E3404DDEF9519B3CD3A431B302B0A6DF25F14374FE1356D\
6D51C245E485B576625E7EC6F44C42E9A637ED6B0BFF5CB6F406B7EDEE386BFB5A89\
9FA5AE9F24117C4B1FE649286651ECE45B3DC2007CB8A163BF0598DA48361C55D39A\
69163FA8FD24CF5F83655D23DCA3AD961C62F356208552BB9ED529077096966D670C\
354E4ABC9804F1746C08CA18217C32905E462E36CE3BE39E772C180E86039B2783A2\
EC07A28FB5C55DF06F4C52C9DE2BCBF6955817183995497CEA956AE515D2261898FA\
051015728E5A8AAAC42DAD33170D04507A33A85521ABDF1CBA64ECFB850458DBEF0A\
8AEA71575D060C7DB3970F85A6E1E4C7ABF5AE8CDB0933D71E8C94E04A25619DCEE3\
D2261AD2EE6BF12FFA06D98A0864D87602733EC86A64521F2B18177B200CBBE11757\
7A615D6C770988C0BAD946E208E24FA074E5AB3143DB5BFCE0FD108E4B82D120A93A\
D2CAFFFFFFFFFFFFFFFF

They both compare exactly to rfc3526.

To calculate the group 16 value, I need to add more digits to pi

scale=4000 /* need more digits for group 16 */
pi=4*a(1)
2^4096 - 2^4032 - 1 + 2^64 * (floor(pi * 2^3966) + 240904) /* group 16 */
FFFFFFFFFFFFFFFFC90FDAA22168C234C4C6628B80DC1CD129024E088A67CC74020B\
BEA63B139B22514A08798E3404DDEF9519B3CD3A431B302B0A6DF25F14374FE1356D\
6D51C245E485B576625E7EC6F44C42E9A637ED6B0BFF5CB6F406B7EDEE386BFB5A89\
9FA5AE9F24117C4B1FE649286651ECE45B3DC2007CB8A163BF0598DA48361C55D39A\
69163FA8FD24CF5F83655D23DCA3AD961C62F356208552BB9ED529077096966D670C\
354E4ABC9804F1746C08CA18217C32905E462E36CE3BE39E772C180E86039B2783A2\
EC07A28FB5C55DF06F4C52C9DE2BCBF6955817183995497CEA956AE515D2261898FA\
051015728E5A8AAAC42DAD33170D04507A33A85521ABDF1CBA64ECFB850458DBEF0A\
8AEA71575D060C7DB3970F85A6E1E4C7ABF5AE8CDB0933D71E8C94E04A25619DCEE3\
D2261AD2EE6BF12FFA06D98A0864D87602733EC86A64521F2B18177B200CBBE11757\
7A615D6C770988C0BAD946E208E24FA074E5AB3143DB5BFCE0FD108E4B82D120A921\
08011A723C12A787E6D788719A10BDBA5B2699C327186AF4E23C1A946834B6150BDA\
2583E9CA2AD44CE8DBBBC2DB04DE8EF92E8EFC141FBECAA6287C59474E6BC05D99B2\
964FA090C3A2233BA186515BE7ED1F612970CEE2D7AFB81BDD762170481CD0069127\
D5B05AA993B4EA988D8FDDC186FFB7DC90A6C08F4DF435C934063199FFFFFFFFFFFF\
FFFF

The above was calculated using GNU bc version 1.06.

	-- Mark

---

From: Daniel Kahn Gillmor <[email-masked]>
Date: Tue, 11 Mar 2014 18:14:55 -0400

> You raise a good point. Folks with 'bc -l' should be able to verify the

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 <[email-masked]>
# 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)

---

From: Henrick Hellström <[email-masked]>
Date: Wed, 12 Mar 2014 02:18:20 +0100

> the modulus structured this way?

- If the least significant word of the modulus equals 2^w-1, Montgomery 
reduction becomes more efficient.
- If the most significant word of the modulus equals 2^w-1, school-book 
reduction and Barrett reduction becomes more efficient.
- Pi usually plays the role of a "nothing-up-my-sleeve" value. If a 
random value had been chosen, there would have been no easy way to 
verify that a specially crafted value hadn't been chosen. For instance, 
it is possible to generate a prime, such that the discrete logarithm of 
selected smooth numbers becomes to known to the entity generating the 
prime, and that would make the discrete logarithm problem easier for 
that entity.

> i haven't yet generated these (mainly due to time):

Python is likely too slow for this. Generating the 8192 bit prime takes 
a couple of hours on a contemporary PC using reasonably optimized native 
code.

---

From: Watson Ladd <[email-masked]>
Date: Tue, 11 Mar 2014 20:56:14 -0700

On Tue, Mar 11, 2014 at 6:18 PM, Henrick Hellström <[email-masked]> wrote:
> - If the least significant word of the modulus equals 2^w-1, Montgomery reduction becomes more efficient.
> - If the most significant word of the modulus equals 2^w-1, school-book reduction and Barrett reduction becomes more efficient.
> - Pi usually plays the role of a "nothing-up-my-sleeve" value. If a random value had been chosen, there would have been no easy way to verify that a specially crafted value hadn't been chosen. For instance, it is possible to generate a prime, such that the discrete logarithm of selected smooth numbers becomes to known to the entity generating the prime, and that would make the discrete logarithm problem easier for that entity.
>
>
>>
> Python is likely too slow for this. Generating the 8192 bit prime takes a couple of hours on a contemporary PC using reasonably optimized native code.

A few hours? The trick is to sieve: the numbers for various N are in
an arithmetic progression with a=2^64. As a result, you can determine
which N are possible because of small prime obstructions, thus
decreasing the amount of work to be done. In fact, I think your first
N is entirely determined by considerations modulo the primes less than
100.

Secondly look at (p-1)/2=q instead. If this is known prime determining
if 2*q+1 is prime is quick using Pocklington's Theorem: you don't need
to do the extensive primality testing.

Sincerely,
Watson Ladd
>
>
> _______________________________________________
> TLS mailing list
> [email-masked]
> https://www.ietf.org/mailman/listinfo/tls

-- 
"Those who would give up Essential Liberty to purchase a little
Temporary Safety deserve neither  Liberty nor Safety."
-- Benjamin Franklin

---

From: Henrick Hellström <[email-masked]>
Date: Wed, 12 Mar 2014 10:53:17 +0100

> to do the extensive primality testing.

Exactly, have you tried this?

For the 8192 bit prime, you are to eliminate 4,743,158 candidates. Suppose:

- you are using a double sieve on both primes that eliminates all but 
0.6-0.7% of the candidates,
- you are doing the Fermat test of the Pocklington criterion first, and 
this eliminates all but 1-2% of the remaining candidates,

then this means your search will still be dominated by some 300,000 
modular exponentiations with a 8192 bit modulus and with 8191 bit 
exponents. Let's say a reasonably optimized implementation performs such 
an exponentiation in 400ms on a contemporary PC. This leaves you with an 
expected running time of 20 min.

---

From: [email-masked] (Niels Möller)
Date: Wed, 12 Mar 2014 15:08:05 +0100

> Exactly, have you tried this?

I couldn't resist, so I have now tried that. Below is a lite program to
search for strong primes of the given form. It generates a 4096-bit
prime in about 30s on my office machine (an intel core i5). I haven't had
the patience to let it complete search for the 8192 prime.

Some comments:

1. This is a quick hack, not thoroughly tested or profiled, etc.

2. It may be advantageous to use an even larger prime table for the
   sieving. Current code overflows somewhere if it is made substantially
   larger, though. Generating the prime table is more or less
   instantaneous.

3. The 8192-bit prime definitely seems a bit expensive to generate.

4. I let the pari/gp compute the pi approximation for me. The program
   could be made more self-contained if one adds some suitable
   function for computing pi to arbitrary precision, but I'm not
   really familiar with that problem.

Regards,
/Niels

/*
    Copyright (C) 2014  Niels Möller

    This program is free software: you can redistribute it and/or modify
    it under the terms of the GNU General Public License as published by
    the Free Software Foundation, either version 3 of the License, or
    (at your option) any later version.

    This program is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    GNU General Public License for more details.

    You should have received a copy of the GNU General Public License
    along with this program.  If not, see <http://www.gnu.org/licenses/>.
*/

#include <assert.h>
#include <stdint.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>

#include <gmp.h>

static void *
xalloc (size_t size)
{
  void *p = malloc (size);
  if (!p)
    {
      fprintf (stderr, "Memory exhausted.\n");
      abort();
    }
  return p;
}

static void *
xrealloc (void *p, size_t size)
{
  p = realloc (p, size);
  if (!p)
    {
      fprintf (stderr, "Memory exhausted.\n");
      abort();
    }
  return p;
}

struct prime_info
{
  uint32_t p;	   /* An odd prime. */ 
  uint32_t inv;    /* p^{-1} mod 2^32 */
  uint32_t limit;  /* floor ((2^32-1) / p) */
};

static struct prime_info *primes;
static unsigned nprimes;

#define DIVISIBLE_P(x, i) ((x) * primes[i].inv <= primes[i].limit)

static uint32_t binvert(uint32_t a)
{
  unsigned char tab[8] = { 1, 11, 13, 7, 9, 3, 5, 15 };
  uint32_t x = tab[(a >> 1) & 7]; /* 4-bit inverse */
  x = 2*x - a*x*x; /* 8 bits */
  x = 2*x - a*x*x; /* 16 bits */
  x = 2*x - a*x*x; /* 32 bits */
  assert (a*x == 1);

  return x;
}

static void
init_table (unsigned limit)
{
  /* Number of odd numbers in the range 3 <= j < limit */
  unsigned sieve_size = limit/2-1;
  char *sieve = xalloc (sieve_size);
  unsigned table_size = 0;
  unsigned i;
  
  primes = NULL;
  nprimes = 0;

  memset (sieve, 1, sieve_size);
  for (i = 0; i < sieve_size; i++)
    {
      unsigned p = 3 + 2*i;
      unsigned j;
      if (sieve[i])
	{
	  /* fprintf (stderr, "%u\n", p); */
	  if (nprimes <= table_size)
	    {
	      table_size = 50 + 3*table_size/2;
	      primes = realloc (primes, table_size);
	    }
	  primes[nprimes].p = p;
	  primes[nprimes].inv = binvert(p);
	  primes[nprimes].limit = (~(uint32_t) 0) / p;
	  nprimes++;
	}
      /* Clear elements corrsponding to p^2, (p+2)^p, ... */
      for (j = (p*p-3)/2; j < sieve_size; j += p)
	sieve[j] = 0;
    }
  free (sieve);
}

/* Find the first strong prime in the sequence start + k step */
void
next_strong_prime(mpz_t r, const mpz_t start, const mpz_t step)
{
  mpz_t t, q, p, pm1, h;

  uint32_t *qmod = xalloc (nprimes * sizeof(*qmod));
  uint32_t *hmod = xalloc (nprimes * sizeof(*hmod));
  unsigned i;

  mpz_inits (t, q, p, pm1, h, NULL);
  /* Some other variants are possible, but exclude them, for
     simplicity */
  assert (mpz_odd_p (start));
  assert (mpz_even_p (step));

  mpz_gcd (t, start, step);
  assert (mpz_cmp_ui (t, 1) == 0);

  mpz_sub_ui (q, start, 1);
  mpz_fdiv_q_2exp (q, q, 1);
  mpz_fdiv_q_2exp (h, step, 1);

  for (i = 0; i < nprimes; i++)
    hmod[i] = mpz_fdiv_ui (h, primes[i].p);

  for (;;)
    {
      /* New q, recompute qmod */
      for (i = 0; i < nprimes; i++)
	qmod[i] = mpz_fdiv_ui (q, primes[i].p);

      /* Limit to avoid overflow in i * hmod[j] */
      for (i = 0; i < 100000; i++)
	{
	  unsigned j;
	  for (j = 0; j < nprimes; j++)
	    {
	      if (DIVISIBLE_P(qmod[j] + i * hmod[j], j))
		/* q not prime */
		goto not_prime;
	      if (DIVISIBLE_P(2*qmod[j] + 1 + 2*i * hmod[j], j))
		/* p not prime */
		goto not_prime;
	    }
	  fprintf (stderr, ".");
	  mpz_addmul_ui (q, h, i);
	  mpz_mul_2exp (pm1, q, 1);
	  mpz_add_ui (p, pm1, 1);

	  /* Since q is updated, we need to reset i and update qmod */
	  for (j = 0; j < nprimes; j++)
	    qmod[j] = (qmod[j] + i * hmod[j]) % primes[j].p;
	  i = 0;

	  /* First check if q prime implies p prime. By pocklington,
	     check that gcd (a^2-1, p) = 1 and a^{2q} = 1 (mod p).
	     Need use only a == 2, and then we already know that a^2-1
	     = 3 doesn't have any common factor with p. */

	  mpz_set_ui (t, 2);
	  mpz_powm (t, t, q, p);
	  /* If t = +/- 1 (mod p), then t^2 = 1, and the conditions of
	     Pocklingtons's theorem (except q prime, which remain to
	     be checked) are satisfied.

	     On the other hand, if p and q are indeed prime, then 1 or
	     -1 are the only possibilities, since t^2 = a^{p-1} = 1
	     (mod p), and those values are the only square roots of 1.
	  */
	  if (mpz_cmp_ui (t, 1) != 0 && mpz_cmp (t, pm1) != 0)
	    goto not_prime;

	  fprintf (stderr, "p");

	  if (!mpz_millerrabin (q, 25))
	    goto not_prime;	    

	  assert (mpz_probab_prime_p (p, 25));
	  fprintf (stderr, "pq\n");

	  mpz_swap (r, p);
	  mpz_clears (t, q, p, pm1, h, NULL);
	  free (hmod);
	  free (qmod);
	  return;
	  
	not_prime:
	  /* Go on with the same i and q */
	  ;
	}
      /* Can this ever happen? */
      mpz_addmul_ui (q, h, i);
      fprintf (stderr, "100000 candidates passed trial division!\n");
    }
}

int
main (int argc, char **argv)
{
  mpz_t pi8192, start, step, t, p;
  int bits;
  int limit = 1000000;

  if (argc < 2)
    {
      fprintf (stderr,"%s BITS [SIEVE-LIMIT]\n", argv[0]);
      return EXIT_FAILURE;
    }
  bits = atoi(argv[1]);
  if (bits < 256 || bits > 8192)
    {
      fprintf (stderr,"Supported range is 256 <= bits < 8192.\n");
      return EXIT_FAILURE;
    }

  if (argc > 2)
    {
      limit = atoi(argv[2]);
      if (limit < 4)
	{
	  fprintf (stderr,"Invalid prime limit.\n");
	  return EXIT_FAILURE;
	}
    }
  init_table(limit);

  fprintf (stderr, "Using %u primes for sieving, largest p = %u\n",
	   nprimes, primes[nprimes-1].p);
  mpz_inits (pi8192, start, step, t, p, NULL);

  /* floor (pi * 2^8192), computed with pari/gp. This is an 8194 bit number */
  mpz_set_str (pi8192,
	       "342668632977872056340614340318593589633845496407221537975716"
	       "613737152546230132865696034849115540419392405119782364448996"
	       "505597009061206839355038690511597649633340755295662641609317"
	       "535642924345851483170828622195425552862989984357380546470618"
	       "654558998494810401062770045070958319118968772694396124916431"
	       "894580767686315880694138861428875902873600542294826287748558"
	       "777680794805812784463450507747144885813404423833320387798321"
	       "609280015466841325521574345612804300925224268961852838014727"
	       "344949392988773343665477596160680355250188118765396205271825"
	       "594335340262310201378513530164474118131690724342410799910757"
	       "546982847641291056024097727491901066776286401345290418249625"
	       "628084349896086834732113313412019501787434146151912994099903"
	       "795562822901556466065108423731653397934665628418225503906500"
	       "252356756723136478156580557983552276992132926102607206902926"
	       "217226186459755888047793134565096032547470753739485581485056"
	       "095310312774412642091527392486496373093667048791992749203017"
	       "039714338536895459352056225848086680604116795009774640833217"
	       "634173864865403767984673414411303165787683355096108212088449"
	       "911140581907950288312790601048394648396324285541310288154248"
	       "299687810226823038493790254875059185044994880417764356982840"
	       "836688145748032677215052979431605921845872542983666091618301"
	       "641191123769382721257794537529591295620218822973870974507625"
	       "710916011007688849347790863162497296922421716445647812915537"
	       "174656134053117238096786067011739714153581530618764175406523"
	       "316433003785370745957635127014688173775323182910793894295533"
	       "869774970638119508426058911918833814428165195037324835873478"
	       "264560759508465057903722354035637978838385867006178927628223"
	       "655893247618439468688098957702197415152483332250099523049205"
	       "947113206854728848346101557957368387613667773778399604898871"
	       "571620512445391665708436224753527427390770571921354058146861"
	       "534703545990059875592521055483476922407124171344793745585445"
	       "612078729092110740281827121644569208704341570355808361384144"
	       "320068653876547846047724572021904389025572441127188782597462"
	       "975206221252835136986819772656964303920606775754918430878535"
	       "831653872047705625068521848678547523189492202180161565293546"
	       "004893947202674258029340791722950503743099899657901284673761"
	       "082932857148614112547296837932346608390996266386438042203034"
	       "975635130970578410146576920597437945752613183601468824949930"
	       "014430184646230304464605951116975119907446615541881480270560"
	       "054980165635722725085125496591712178571547448110006296573138"
	       "030681065709805770745255837670620283490202933751846093959476"
	       "8886090", 10);

  mpz_setbit (step, 64);
  mpz_sub_ui (t, step, 1);
  mpz_mul_2exp (start, t, bits - 64);
  mpz_add (start, start, t);

  mpz_fdiv_q_2exp (t, pi8192, (8194 + 128 - bits));
  assert (mpz_sizeinbase (t, 2) == bits - 128);
  mpz_mul_2exp (t, t, 64);
  mpz_add (start, start, t);

  next_strong_prime (p, start, step);

  mpz_sub (t, p, start);
  mpz_fdiv_q_2exp (t, t, 64);
  gmp_printf ("p = 0x%Zx\nt = %Zd\n", p, t);

  mpz_clears (pi8192, start, step, t, p, NULL);

  free (primes);
  return EXIT_SUCCESS;
}

-- 
Niels Möller. PGP-encrypted email is preferred. Keyid C0B98E26.
Internet email is subject to wholesale government surveillance.

---

From: [email-masked] (Niels Möller)
Date: Wed, 12 Mar 2014 20:29:18 +0100

> 1. This is a quick hack, not thoroughly tested or profiled, etc.

And memory allocation was totally broken. I've put an updated version at
http://www.lysator.liu.se/~nisse/misc/next-strong-prime.c

Regards,
/Niels

---

From: [email-masked] (Niels Möller)
Date: Wed, 12 Mar 2014 22:24:35 +0100

> 3. The 8192-bit prime definitely seems a bit expensive to generate.

I've now ran the updated version of the program on a more high-end pc.
For 8192 bits, it completed in about 52 minutes. Output:

p = 0xffffffffffffffffc90fdaa22168c234c4c6628b80dc1cd129024e088a67cc74020bbea63b139b22514a08798e3404ddef9519b3cd3a431b302b0a6df25f14374fe1356d6d51c245e485b576625e7ec6f44c42e9a637ed6b0bff5cb6f406b7edee386bfb5a899fa5ae9f24117c4b1fe649286651ece45b3dc2007cb8a163bf0598da48361c55d39a69163fa8fd24cf5f83655d23dca3ad961c62f356208552bb9ed529077096966d670c354e4abc9804f1746c08ca18217c32905e462e36ce3be39e772c180e86039b2783a2ec07a28fb5c55df06f4c52c9de2bcbf6955817183995497cea956ae515d2261898fa051015728e5a8aaac42dad33170d04507a33a85521abdf1cba64ecfb850458dbef0a8aea71575d060c7db3970f85a6e1e4c7abf5ae8cdb0933d71e8c94e04a25619dcee3d2261ad2ee6bf12ffa06d98a0864d87602733ec86a64521f2b18177b200cbbe117577a615d6c770988c0bad946e208e24fa074e5ab3143db5bfce0fd108e4b82d120a92108011a723c12a787e6d788719a10bdba5b2699c327186af4e23c1a946834b6150bda2583e9ca2ad44ce8dbbbc2db04de8ef92e8efc141fbecaa6287c59474e6bc05d99b2964fa090c3a2233ba186515be7ed1f612970cee2d7afb81bdd762170481cd0069127d5b05aa993b4ea988d8fddc186ffb7d!
 c90a6c08f4df435c93402849236c3fab4d27c7026c1d4dcb2602646dec9751e763dba37bdf8ff9406ad9e530ee5db382f413001aeb06a53ed9027d831179727b0865a8918da3edbebcf9b14ed44ce6cbaced4bb1bdb7f1447e6cc254b332051512bd7af426fb8f401378cd2bf5983ca01c64b92ecf032ea15d1721d03f482d7ce6e74fef6d55e702f46980c82b5a84031900b1c9e59e7c97fbec7e8f323a97a7e36cc88be0f1d45b7ff585ac54bd407b22b4154aacc8f6d7ebf48e1d814cc5ed20f8037e0a79715eef29be32806a1d58bb7c5da76f550aa3d8a1fbff0eb19ccb1a313d55cda56c9ec2ef29632387fe8d76e3c0468043e8f663f4860ee12bf2d5b0b7474d6e694f91e6dbe115974a3926f12fee5e438777cb6a932df8cd8bec4d073b931ba3bc832b68d9dd300741fa7bf8afc47ed2576f6936ba424663aab639c5ae4f5683423b4742bf1c978238f16cbe39d652de3fdb8befc848ad922222e04a4037c0713eb57a81a23f0c73473fc646cea306b4bcbc8862f8385ddfa9d4b7fa2c087e879683303ed5bdd3a062b3cf5b3a278a66d2a13f83f44f82ddf310ee074ab6a364597e899a0255dc164f31cc50846851df9ab48195ded7ea1b1d510bd7ee74d73faf36bc31ecfa268359046f4eb879f924009438b481c6cd7889a002ed5ee382bc9190da6fc026e47955!
 8e4475677e9aa9e3050e2765694dfc81f56e880b96e7160c980dd98edd3dff!
 fffffffffffffff
k = 4743157
elapsed time 3133.28s

Which after a cursory look appears to agree with
https://www.rfc-editor.org/rfc/rfc3526.txt (I guess the 47473157 vs
4743158 is an off-by-one error in my code. Does the spec's brackets in
in "[2^8062 pi]" denote floor or round-to-nearest?

Regards,
/Niels

---

From: Daniel Kahn Gillmor <[email-masked]>
Date: Wed, 12 Mar 2014 17:30:42 -0400

> in "[2^8062 pi]" denote floor or round-to-nearest?

i was using floor for my calculations, not round-to-nearest.

	--dkg

---

From: Jeffrey Hutzelman <[email-masked]>
Date: Wed, 12 Mar 2014 17:50:48 -0400

On Wed, 2014-03-12 at 22:24 +0100, Niels Möller wrote:
> p = 0xffffffffffffffffc90fdaa22168c234c4c6628b80dc1cd129024e088a67cc74020bbea63b139b22514a08798e3404ddef9519b3cd3a431b302b0a6df25f14374fe1356d6d51c245e485b576625e7ec6f44c42e9a637ed6b0bff5cb6f406b7edee386bfb5a899fa5ae9f24117c4b1fe649286651ece45b3dc2007cb8a163bf0598da48361c55d39a69163fa8fd24cf5f83655d23dca3ad961c62f356208552bb9ed529077096966d670c354e4abc9804f1746c08ca18217c32905e462e36ce3be39e772c180e86039b2783a2ec07a28fb5c55df06f4c52c9de2bcbf6955817183995497cea956ae515d2261898fa051015728e5a8aaac42dad33170d04507a33a85521abdf1cba64ecfb850458dbef0a8aea71575d060c7db3970f85a6e1e4c7abf5ae8cdb0933d71e8c94e04a25619dcee3d2261ad2ee6bf12ffa06d98a0864d87602733ec86a64521f2b18177b200cbbe117577a615d6c770988c0bad946e208e24fa074e5ab3143db5bfce0fd108e4b82d120a92108011a723c12a787e6d788719a10bdba5b2699c327186af4e23c1a946834b6150bda2583e9ca2ad44ce8dbbbc2db04de8ef92e8efc141fbecaa6287c59474e6bc05d99b2964fa090c3a2233ba186515be7ed1f612970cee2d7afb81bdd762170481cd0069127d5b05aa993b4ea988d8fddc186ffb!
 7d!
>  c90a6c08f4df435c93402849236c3fab4d27c7026c1d4dcb2602646dec9751e763dba37bdf8ff9406ad9e530ee5db382f413001aeb06a53ed9027d831179727b0865a8918da3edbebcf9b14ed44ce6cbaced4bb1bdb7f1447e6cc254b332051512bd7af426fb8f401378cd2bf5983ca01c64b92ecf032ea15d1721d03f482d7ce6e74fef6d55e702f46980c82b5a84031900b1c9e59e7c97fbec7e8f323a97a7e36cc88be0f1d45b7ff585ac54bd407b22b4154aacc8f6d7ebf48e1d814cc5ed20f8037e0a79715eef29be32806a1d58bb7c5da76f550aa3d8a1fbff0eb19ccb1a313d55cda56c9ec2ef29632387fe8d76e3c0468043e8f663f4860ee12bf2d5b0b7474d6e694f91e6dbe115974a3926f12fee5e438777cb6a932df8cd8bec4d073b931ba3bc832b68d9dd300741fa7bf8afc47ed2576f6936ba424663aab639c5ae4f5683423b4742bf1c978238f16cbe39d652de3fdb8befc848ad922222e04a4037c0713eb57a81a23f0c73473fc646cea306b4bcbc8862f8385ddfa9d4b7fa2c087e879683303ed5bdd3a062b3cf5b3a278a66d2a13f83f44f82ddf310ee074ab6a364597e899a0255dc164f31cc50846851df9ab48195ded7ea1b1d510bd7ee74d73faf36bc31ecfa268359046f4eb879f924009438b481c6cd7889a002ed5ee382bc9190da6fc026e479!
 55!
> in "[2^8062 pi]" denote floor or round-to-nearest?

Neither 2412 nor 3526 says what they mean by that notation, but I'm
assuming floor.  However, in that case your result agrees with the spec.

The value given in the spec is

2^8192 - 2^8128 - 1 + 2^64 * { [2^8062 pi] + 4743158 }

This can be rewritten as:

2^8192 - 2^8128 + 2^64 * { [2^8062 pi] + 4743157 } + 2^64 - 1

The first two terms give the high-order 64 bits (all 1).
The next term gives the middle bits, and agrees with your result.
The last two terms give the low-order 64 bits.

-- Jeff