Instantly share code, notes, and snippets.

@grocid /ic.sage
Last active Sep 5, 2018

Embed
What would you like to do?
p = 3122882963 # 32-bit safe prime
#p = 13019428505045596307
R = IntegerModRing(p - 1)
Field = GF(p)
g = Field(2*3*5*7*11*17*19)
exp = 1234562364 # just some bang-on-the-keyboard exponent
h = g**exp
B = 199 # there are 45 primes < B
# find basis using approximation
def search_for_basis(x, bound):
L = list(factor(int(x)))
z = 1
for q, e in L:
if q < B:
z *= q**e
else:
break
if x/z < int(sqrt(sqrt(p))): # set to 2 to get normal ic, sqrt(p) to approx
v = vector(R, B)
for q, e in L:
if q < B:
v[q] = e
return v, x / z
return None, x / z
def append_row(M, v):
return matrix(Field, M.rows() + [v])
def solve(A, v, M, N):
LL = []
for q, e in list(factor(p - 1)):
AA = matrix(GF(q), A.rows())
vv = matrix(GF(q), list(v))
LL.append((AA.solve_left(vv), q, e))
UU = []
for xx, yy in zip(LL[0][0].list(), LL[1][0].list()):
UU.append(crt(int(xx), int(yy), int(LL[0][1]), int(LL[1][1])))
return vector(R, UU)
def check(D, J):
prod = 1
for q, e in zip(D.list(), J.list()):
prod *= pow(q, e, p)
return prod % p
# first stage: generate basis
print "First stage..."
gg = 1
# we set this to gf(p) although it should be IntegerModRing(p-1)
A = matrix(Field, 0, B)
D = matrix(Field, 0, 1)
J = matrix(IntegerModRing(p - 1), 0, 1)
jj = 0
while True:
v, d = search_for_basis(gg, sqrt(p))
if v is not None:
AA = append_row(A, v)
if AA.rank() > A.rank():
A = AA
D = append_row(D, [d])
J = append_row(J, [jj])
print "Current rank", AA.rank(), jj
if A.rank() == 45:
print "no linear indep"
break
gg = gg * g
jj += 1
#print jj
print "Second stage..."
# second stage: generate target
j = 0
while True:
while True:
u, d = search_for_basis(h, sqrt(p))
h = h * g
j += 1
if u is not None:
break
try:
xx = vector(R, J.transpose().list()) * \
vector(R, list(solve(A, u, J, D)))# % (p - 1)
except Exception as e:
print "nope...", j, h, e
continue
break
print "Looking for {}, found {}".format(exp % (p - 1), (xx - j + 1) % (p-1))
print factor(int((xx - j + 1) ))
print J.transpose()
print D.transpose()
print check( D, J)
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment