ic.sage

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)