from operator import mul
from sage.matrix.all import *
def factorbaseexponents(fb,f):
fbex = [0]*len(fb) # list of zeros
for t in list(f): # loop over tuples in factors
i = fb.index(t[0])
fbex[i] = t[1]
return fbex
def kraitchik(N,y,r):
fb = prime_range(y)
print 'Length of factorbase: ' + str(len(fb))
x = int(ceil(sqrt(N)))
expmat = []
xlist = []
# find the matrix
xsqr = x*x
for t in range(1,r+1):
d = xsqr % N
f = factor(d)
if f[-1][0]<= y:
fbex = factorbaseexponents(fb,f)
expmat.append(fbex)
xlist.append(x)
xsqr = xsqr + 2*x + 1
x = x+1
print 'Length of list of x: ' + str(len(xlist))
# get null basis
A = matrix(GF(2), len(expmat), len(fb) , expmat)
vecs = A.left_kernel().basis()
print 'Size of null space: ' + str(len(vecs))
# check for valid u,v
for vec in vecs:
u = 1
vsqr = 1
for i in range(len(vec)):
if vec[i] == 1:
x = xlist[i]
u = (u*x)
mulv = reduce(mul, [a**b for a,b in zip(fb,expmat[i])],1)
vsqr = vsqr*mulv
v = int(sqrt(vsqr))
if ((((u+v)%N) != 0) and (((u-v)%N) != 0)):
print gcd(u-v,N), gcd(u+v,N)
return
return 'Failed.'