18.783 Lecture 11 Index calculus.sagews

# generates a random Sophie-Germaine prime (p=2*q+1 with q prime)
# we use these to make Fp^times as close to a cyclic group of prime order as possible
# this means a Pohlig-Hellman attack on Fp^* will have very little benefit
def random_sg_prime(bits):
    while true:
        p=random_prime(2^bits,2^(bits-1))
        if is_prime(2*p+1): break
    return p

def L(a,c,N):
    z = ceil(exp(c*log(N)^a*log(log(N))^(1-a)))
    if z%2: z+=1
    return ceil(z)

def trial_factor(x,B):
    f=[]
    lastp=0
    while x > 1:
        p=x.trial_division(B)
        if p > B: return []
        if p == lastp: f[-1][1] += 1
        else: f.append([p,1])
        x = x.divide_knowing_divisible_by(p)
        lastp = p
    return f
    
def dlog(a,b,p):
    B=ceil(2*L(1/2,1/2,p))
    print "smoothness bound is ", B
    pi=[q for q in primes(B)]
    pimap=[0 for i in range(0,B)]
    i=0
    for q in primes(B):
        pimap[q] = i
        i += 1
    n=len(pi)
    print "factor base has %d elements, searching for relations..."%n
    
    M=matrix(ZZ,n+2,n+2,sparse=true)
    t = cputime()
    bi = 1/b
    i = 0
    while i < M.nrows():
        e = randint(1,p-1)
        f=trial_factor(ZZ(a^e*bi),B)
        if f and f[-1][0] < B:
            for q in f:
                M[i,pimap[q[0]]]=q[1]
            M[i,n]=1; M[i,n+1]=e
            i += 1
    print "Found %d relations in %.2f s, attempting to solve system"%(i,cputime()-t)
    q = ZZ((p-1)/2)
    print M
    Mq=M.change_ring(GF(q)).echelon_form()
    M2=M.change_ring(GF(2)).echelon_form()
    print Mq
    print M2
    i = n; j = n
    while not Mq[i,n]: i-=1
    while not M2[j,n]: j-=1
    print i,j,n
    return ((q+1)*ZZ(Mq[i,n+1])+q*ZZ(M2[j,n+1]))%(2*q)
    

%time
p=2*random_sg_prime(50)+1
print p
F=GF(p)
a=F.random_element()
while a.multiplicative_order() != p-1: a=F.random_element()
x=randint(1,p-1)
b=a^x
print b, "=", a, "^", x
y=dlog(a,b,p)
print "dlog is", y
assert x==y

1700823059257403 926235546846520 = 474138947680118 ^ 1328456286604230 smoothness bound is 536 factor base has 99 elements, searching for relations... Found 101 relations in 13.94 s, attempting to solve system 101 x 101 sparse matrix over Integer Ring 101 x 101 sparse matrix over Finite Field of size 850411529628701 101 x 101 sparse matrix over Finite Field of size 2 97 97 99 dlog is 1328456286604230 CPU time: 14.59 s, Wall time: 14.75 s