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 %d" % 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)