Creating colliding x509 certificates: Crafting RSA moduli step

Question

I am trying to generate 2 x509 certificates with the same signature but different values in the common name field, based on md5 collisions, as it was specified in this paper (page 7).

Now I have successfully created 2 files with the same md5 value, that are valid certificates, but the public key field doesn't contain valid RSA moduli yet. According to the paper, they can be generated using the method described in this paper (page 2, step 4), with sizes for p and q of 1976 and 6216 bits in one hour on a regular laptop. It doesn't say how to adjust the values in the second paper exactly (which was for 2048 bit moduli). Moreover, hashclash produced 9 md5 blocks instead of 8 in the paper, so it might not work with these sizes.

Here's the algorithm from the second paper for a quick reference: www.imgur.com/Z35J8Af

The difference is that $b_1$ and $b_2$ are 4704 bits each (9 * 512 + 96), so $b$ should be 3488 bits.

I simply tried different values that might make sense (1024, 1976, 1740, 1730, ...), for several hours, but the script never finished. It worked for 1024 bit moduli though.

Now, what I need to know is:

1. What size to pick for p, if it's not one of the values that I tried
2. Or is there another parameter I might have to change in the script below ?

Here is the script:

from mymath import *
import random
import sys

e = 65537

def read_file_as_number(name, skip):

    f = open(name)
    n = 0
    byte = f.read(1)

    while byte != "":
        if skip == 0:
            n = n << 8 | ord(byte)
        else:
            skip -= 1
        byte = f.read(1)

    f.close()

    return n

b1 = read_file_as_number("file1_9.bin", 244)
b2 = read_file_as_number("file2_9.bin", 244)

print hex(b1)
print hex(b2)

b_size = 3488

print "b length in bits:", b_size

b1 <<= b_size
b2 <<= b_size

def random_prime(bitsize):
    x = random.randint(0, 1 << (bitsize - 1))

    return next_prime(x)


# random prime x, where x - 1 is coprime to n
def random_prime_coprime(bitsize, n):

    x = random_prime(bitsize)

    if gcd(x - 1, n) == 1:
        return x
    else:
        return random_prime_coprime(bitsize, n)

# chinese remainder theorem
def crt(a1, a2, n1, n2):
    N = n1 * n2

    inv1 = mod_inv(n2, n1)
    inv2 = mod_inv(n1, n2)

    return -(a1 * inv1 * n2 + a2 * inv2 * n1) % N

p_size = 1740

print "p size:", p_size

maxb = 1 << b_size

while True:
    b = 0
    k = 0
    p1 = random_prime_coprime(p_size, e)
    p2 = random_prime_coprime(p_size, e)

    b0 = crt(b1, b2, p1, p2)

    assert gcd(p1, b1 + b0) == p1
    assert gcd(p2, b2 + b0) == p2
    assert b0 < p1 * p2

    while b < maxb:

        b = b0 + k * p1 * p2

        q1 = (b1 + b) / p1
        q2 = (b2 + b) / p2

        if gcd(q1 - 1, e) == 1 and gcd(q2 - 1, e) == 1 and is_prime(q1) and is_prime(q2):

            print b0

            print "DONE"
            print "p1 =", p1
            print "p2 =", p2
            print "q1 =", q1
            print "q2 =", q2

            print "n1 =", b1 + b
            print "n2 =", b2 + b

            exit()

        k += 1
        sys.stdout.write('-')
        sys.stdout.flush()

    sys.stdout.write('.')
    sys.stdout.flush()