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