Skip to content

Instantly share code, notes, and snippets.

@grocid 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
You can’t perform that action at this time.
You signed in with another tab or window. Reload to refresh your session. You signed out in another tab or window. Reload to refresh your session.