A function in my factor package
return zeros instead of the correct value.
It did this for several times, and then returned the correct value.
I inserted print statements in the function and in the calling routine
which document this behavior.
The calling routine is function function fermat.
The function being called is function strongfac.
To see the test results , run the function testrange
with parameters
testrange(10**14+37,10**14+37).
Since it takes quite a lot of work to make a copy that survives the
email
removal of leading blanks
I attach the package of function routines here.
Results as run on my machine:
Python 2.4.2 (#67, Sep 28 2005, 12:41:11) [MSC v.1310 32 bit (Intel)] on
win32
Type "copyright", "credits" or "license()" for more information.
IDLE 1.1.2
>>> import factor34
>>> from factor34 import testrange
>>> testrange(10**14+37,10**14+37)
Begin calculating constants to be used in factoring.
Finished calculating constants to be used in factoring.
Search for a factor of 100000000000037
Could not find factors by strong probable prime test using witnesses <
10000.
Try to find factors by using power of 2 mod z as a witness.
Retrieved from strongfac, face = [0, 0, 0]
In fermat: k = 0 x = 0 y = 0 face = [0, 0, 0]
Retrieved from strongfac, face = [0, 0, 0]
In fermat: k = 1 x = 0 y = 0 face = [0, 0, 0]
Retrieved from strongfac, face = [0, 0, 0]
In fermat: k = 2 x = 0 y = 0 face = [0, 0, 0]
Retrieved from strongfac, face = [0, 0, 0]
In fermat: k = 3 x = 0 y = 0 face = [0, 0, 0]
Retrieved from strongfac, face = [0, 0, 0]
In fermat: k = 4 x = 0 y = 0 face = [0, 0, 0]
Retrieved from strongfac, face = [0, 0, 0]
In fermat: k = 5 x = 0 y = 0 face = [0, 0, 0]
Retrieved from strongfac, face = [0, 0, 0]
In fermat: k = 6 x = 0 y = 0 face = [0, 0, 0]
Retrieved from strongfac, face = [0, 0, 0]
In fermat: k = 7 x = 0 y = 0 face = [0, 0, 0]
Retrieved from strongfac, face = [0, 0, 0]
In fermat: k = 8 x = 0 y = 0 face = [0, 0, 0]
Retrieved from strongfac, face = [0, 0, 0]
In fermat: k = 9 x = 0 y = 0 face = [0, 0, 0]
Retrieved from strongfac, face = [0, 0, 0]
In fermat: k = 10 x = 0 y = 0 face = [0, 0, 0]
Retrieved from strongfac, face = [0, 0, 0]
In fermat: k = 11 x = 0 y = 0 face = [0, 0, 0]
Retrieved from strongfac, face = [0, 0, 0]
In fermat: k = 12 x = 0 y = 0 face = [0, 0, 0]
In strongfac
Found1: x = 53799857
Found1: y = 1858741
face = [53799857L, 1858741L, 0]
Found1: factor by strong probable prime test, using witnes 15306543515214
.
Found1: factors are 53799857 1858741
Found1: returning [53799857L, 1858741L, 0] to fermat routine.
Retrieved from strongfac, face = [0, 0, 0]
In fermat: k = 13 x = 0 y = 0 face = [0, 0, 0]
In strongfac
Found1: x = 53799857
Found1: y = 1858741
face = [53799857L, 1858741L, 0]
Found1: factor by strong probable prime test, using witnes 12123044576953
.
Found1: factors are 53799857 1858741
Found1: returning [53799857L, 1858741L, 0] to fermat routine.
Retrieved from strongfac, face = [0, 0, 0]
In fermat: k = 14 x = 0 y = 0 face = [0, 0, 0]
In strongfac
Found1: x = 53799857
Found1: y = 1858741
face = [53799857L, 1858741L, 0]
Found1: factor by strong probable prime test, using witnes 45391315949900
.
Found1: factors are 53799857 1858741
Found1: returning [53799857L, 1858741L, 0] to fermat routine.
Retrieved from strongfac, face = [0, 0, 0]
In fermat: k = 15 x = 0 y = 0 face = [0, 0, 0]
In strongfac
Found1: x = 53799857
Found1: y = 1858741
face = [53799857L, 1858741L, 0]
Found1: factor by strong probable prime test, using witnes 59571344259390
.
Found1: factors are 53799857 1858741
Found1: returning [53799857L, 1858741L, 0] to fermat routine.
Retrieved from strongfac, face = [0, 0, 0]
In fermat: k = 16 x = 0 y = 0 face = [0, 0, 0]
In strongfac
Found1: x = 53799857
Found1: y = 1858741
face = [53799857L, 1858741L, 0]
Found1: factor by strong probable prime test, using witnes 78029752396948
.
Found1: factors are 53799857 1858741
Found1: returning [53799857L, 1858741L, 0] to fermat routine.
Retrieved from strongfac, face = [0, 0, 0]
In fermat: k = 17 x = 0 y = 0 face = [0, 0, 0]
In strongfac
Found1: x = 53799857
Found1: y = 1858741
face = [53799857L, 1858741L, 0]
Found1: factor by strong probable prime test, using witnes 35863146075772
.
Found1: factors are 53799857 1858741
Found1: returning [53799857L, 1858741L, 0] to fermat routine.
Retrieved from strongfac, face = [0, 0, 0]
In fermat: k = 18 x = 0 y = 0 face = [0, 0, 0]
In strongfac
Found1: x = 53799857
Found1: y = 1858741
face = [53799857L, 1858741L, 0]
Found1: factor by strong probable prime test, using witnes 19712913203085
.
Found1: factors are 53799857 1858741
Found1: returning [53799857L, 1858741L, 0] to fermat routine.
Retrieved from strongfac, face = [0, 0, 0]
In fermat: k = 19 x = 0 y = 0 face = [0, 0, 0]
In strongfac
Found1: x = 53799857
Found1: y = 1858741
face = [53799857L, 1858741L, 0]
Found1: factor by strong probable prime test, using witnes 14607414373499
.
Found1: factors are 53799857 1858741
Found1: returning [53799857L, 1858741L, 0] to fermat routine.
Retrieved from strongfac, face = [0, 0, 0]
In fermat: k = 20 x = 0 y = 0 face = [0, 0, 0]
In strongfac
Found1: x = 53799857
Found1: y = 1858741
face = [53799857L, 1858741L, 0]
Found1: factor by strong probable prime test, using witnes 57761947468766
.
Found1: factors are 53799857 1858741
Found1: returning [53799857L, 1858741L, 0] to fermat routine.
Retrieved from strongfac, face = [0, 0, 0]
In fermat: k = 21 x = 0 y = 0 face = [0, 0, 0]
In strongfac
Found1: x = 53799857
Found1: y = 1858741
face = [53799857L, 1858741L, 0]
Found1: factor by strong probable prime test, using witnes 75604347243674
.
Found1: factors are 53799857 1858741
Found1: returning [53799857L, 1858741L, 0] to fermat routine.
Retrieved from strongfac, face = [0, 0, 0]
In fermat: k = 22 x = 0 y = 0 face = [0, 0, 0]
In strongfac
Found1: x = 53799857
Found1: y = 1858741
face = [53799857L, 1858741L, 0]
Found1: factor by strong probable prime test, using witnes 3561873332543
Found1: factors are 53799857 1858741
Found1: returning [53799857L, 1858741L, 0] to fermat routine.
In strongfac
Found1: x = 1858741
Found1: y = 53799857
face = [1858741L, 53799857L, 0]
Found1: factor by strong probable prime test, using witnes 96058894729069
.
Found1: factors are 1858741 53799857
Found1: returning [1858741L, 53799857L, 0] to fermat routine.
Retrieved from strongfac, face = [1858741L, 53799857L, 0]
In fermat: k = 23 x = 1858741 y = 53799857 face = [1858741L,
53799857L, 0]
index of torial is 23
z = 100000000000037 x = 1858741 y = 53799857 difficulty = 0
>>>
Kermit Rose < [EMAIL PROTECTED] >
def testrange(za,zb):
print " Begin calculating constants to be used in factoring."
pmap = setpmap()
plist = setplist(pmap)
sqlist = setsqlist(plist)
sqresidue = setsqresidue(sqlist)
sqroot = setsqroot(sqlist)
torials = range(100)
for k in range(100):
mpy = 1
for j in range(1,101):
mpy = mpy * (100 * k + j)
torials[k] = mpy
print " Finished calculating constants to be used in factoring."
for z in range(za,zb+1,2):
testz(z,pmap,plist,sqlist,sqresidue,torials)
def ABCD(z):
D = (z-1)/2
if D%2 == 0:
stlist = [ [0,0],[1,1]]
cst = 2
else:
stlist = [ [0,1],[1,0] ]
cst = 2
qlist = [ [0,2,2,1,1,D,cst,stlist ] ]
cabcd = 0
maxcabcd = 0
stop = 0
level = 0
stop = 0
while stop == 0:
# print " "
# print " cabcd = ",cabcd
# print " qlist[cabcd] = ",qlist[cabcd]
stop1 = 0
while stop1 == 0:
# print " "
# print " cabcd = ",cabcd
if cabcd < 0:
break
cst = qlist[cabcd][6]
# print " cst = ",cst
qlist[cabcd][6] = cst = cst-1
# print " new cst = ",cst
A = qlist[cabcd][2]
D = qlist[cabcd][5]
# print " A = ",A," D = ",D," cst = ",cst
flag = 0
if cst < 0:
flag = 1
if (A - z) > 0:
flag = 1
if (D+z) < 0:
flag = 1
if flag ==0:
break
if flag == 1:
cabcd = cabcd - 1
# print " Backtrack to previous qlist entry. "
# print " cabcd = ",cabcd
# if cabcd > -1:
# print " qlist[cabcd] = ",qlist[cabcd]
if cabcd < 0:
break
level = level + 1
# print " "
# print " cabcd = ",cabcd
# print " qlist[cabcd] = ",qlist[cabcd]
# print " cabcd = ",cabcd," cst = ",cst
# print " qlist[cabcd] = ",qlist[cabcd]
s = qlist[cabcd][7][cst][0]
t = qlist[cabcd][7][cst][1]
A = qlist[cabcd][2]
B = qlist[cabcd][3]
C = qlist[cabcd][4]
D = qlist[cabcd][5]
p = qlist[cabcd][1]
# print " cabcd = ",cabcd, " p = ",p," A = ",A," B = ",B," C = ",C," D =
",D," cst = ",cst," s = ",s," t = ",t
# print " qlist[cabcd] = ",qlist[cabcd]
testz = A * D + B * C
if testz != z:
print " Error. z should have be = A * D + B * C ", " A * D + B * C
= ",testz," z = ",z
if D == 0:
fac = testfac(z,B)
if fac[0] != 0:
fac[2] = [p,0,level]
return fac
fac = testfac(z,B)
if fac[0] != 0:
fac[2] = [p,1,level]
return fac
fac = testfac(z,C)
if fac[0] != 0:
fac[2] = [p,2,level]
return fac
fac = testfac(z,kabs(A - B) )
if fac[0] != 0:
fac[2] = [p,3,level]
return fac
fac = testfac(z,kabs(A+B))
if fac[0] != 0:
fac[2] = [p,4,level]
return fac
fac = testfac(z,kabs(A-C))
if fac[0] != 0:
fac[2] = [p,5,level]
return fac
fac = testfac(z,kabs(D - B))
if fac[0] != 0:
fac[2] = [p,6,level]
return fac
fac = testfac(z,kabs(D + B))
if fac[0] != 0:
fac[2] = [p,7,level]
return fac
fac = testfac(z,kabs(D - C))
if fac[0] != 0:
fac[2] = [p,8,level]
return fac
fac = testfac(z,kabs(D + C))
if fac[0] != 0:
fac[2] = [p,9,level]
return fac
fac = testfac(z,kabs(A + B + C - D))
if fac[0] != 0:
fac[2] = [p,10,level]
return fac
fac = testfac(z,kabs(A + B - C + D))
if fac[0] != 0:
fac[2] = [p,11,level]
return fac
fac = testfac(z,kabs(A - B + C + D))
if fac[0] != 0:
fac[2] = [p,12,level]
return fac
fac = testfac(z,kabs(-A + B + C + D))
if fac[0] != 0:
fac[2] = [p,13,level]
return fac
# print " "
# print " calculate next set of A,B,C,D,s,t"
# print " s = ",s," t = ",t," A = ",A," B = ",B, " C = ",C," D = ",D
A,B,C,D = A*p,B + A * t,C + A*s,(D - A * s * t - B * s - C * t)/p
# print " A = ",A," B = ",B," C = ",C," D = ",D
cst = 0
stlist = []
Dtest = D%p
for s in range(p):
for t in range(p):
test = (A * s * t + B * s + C * t)%p
if test == Dtest:
stlist.append([s,t])
cst = cst + 1
# print " Adding solution ", " cst = ",cst, " [s,t] = ",[s,t]
cabcd = cabcd + 1
# print " Added ",cst," solutions."," cabcd = ",cabcd
if cabcd > maxcabcd:
qlist.append( [cabcd,p,A,B,C,D,cst,stlist] )
qlist[cabcd][6] = cst
maxcabcd = maxcabcd + 1
else:
qlist[cabcd] = [cabcd,p,A,B,C,D,cst,stlist]
qlist[cabcd][6] = cst
# print " qlist[cabcd] = ",qlist[cabcd]
# print " "
# print " cabcd = ",cabcd
# print " qlist[cabcd] = ",qlist[cabcd]
return [1,z,[p,14,level]]
def testfac(z,w):
if w == 0:
fac = [0,0,0]
return fac
x = gcd(z,w)
if x > 1:
y = z/x
fac = [x,y,0]
return fac
fac = [0,0,0]
return fac
def kabs(w):
if w < 0:
return -w
return w
def testz(z,pmap,plist,sqlist,sqresidue,torials):
fac = fermat(z,1,pmap,plist,sqlist,sqresidue,torials)
print " z = ",z," x = ",fac[0]," y = ",fac[1]," difficulty = ",fac[2]
def setsqlist(plist):
sqlist = plist[0:]
sqlist[1] = 8
return sqlist
def setsqresidue(sqlist):
sqresidue = range(1230)
for j in range(1,1230):
p = sqlist[j]
sqresidue[j] = range(p)
for k in range(p):
sqresidue[j][k] = -1
# print " j = ",j," p = ",p," sqresidue[j] = ",sqresidue[j]
half = (p+1)/2
for k in range(half):
k2 = k*k%p
sqresidue[j][k2] = 1
# print " j = ",j," p = ",p," sqresidue[j] = ",sqresidue[j]
# print " "
return sqresidue
def setsqroot(sqlist):
sqroot = range(1230)
sqroot[1] = [0,1,-1,-1,2,-1,-1,-1]
for j in range(2,1230):
p = sqlist[j]
sqroot[j] = range(p)
for k in range(p):
sqroot[j][k] = -1
half = (p+1)/2
for k in range(half):
k2 = k*k%p
sqroot[j][k2] = k
return sqroot
def ifprobsq(z,sqlist,sqresidue):
for k in range(1,1230):
p = sqlist[k]
r = z%p
q = sqresidue[k][r]
if q == -1:
return -1
return 1
def fermat(z,maxdiff,pmap,plist,sqlist,sqresidue,torials):
# maxdiff should be set to one less than number of expected factors.
print " "
print " Search for a factor of ",z
if maxdiff < 1:
maxdiff = 1
if z < 0:
fac = range(3)
fac[0] = -1
fac[1] = z
fac[2] = 0
print " Found factors from z is negative."
return fac
if z == 0:
fac = range(3)
fac[0] = 0
fac[1] = 0
fac[2] = 0
print " Found factors from z == 0."
return fac
if z == 1:
fac = range(3)
fac[0] = 1
fac[1] = 1
fac[2] = 0
print " Found factors from z == 1."
return fac
if z%2 == 0:
fac = range(3)
fac[0] = 2
fac[1] = z/2
fac[2] = 0
print " Found factors from z is even."
return fac
if z < 10000:
if pmap[z] == z:
fac = range(3)
fac[0] = 1
fac[1] = z
fac[2] = 0
print " Found factors by z is prime < 10000."
return fac
q = ifprime(z,plist)
if q == 1:
fac = range(3)
fac[0] = 1
fac[1] = z
fac[2] = 1
print " Found factors by z is probabble prime > 10000."
return fac
for k in range(2,1230):
p = plist[k]
if z%p == 0:
fac = range(3)
fac[0] = p
fac[1] = z/p
fac[2] = k
print " Found factors by z is divisible by prime < 10000."
return fac
print " "
fac = miller(z)
if fac[0] != 0:
return fac
for j in range(2,6):
print " "
print " Try to find factors by using power of ",j," mod z as a witness."
w2 = j
for k in range(100):
w = pow(j,torials[k],z)
face = strongfac(z,w)
print " Retrieved from strongfac, face = ",face
x = face[0]
y = face[1]
print " In fermat: k = ",k," x = ",x," y = ",y," face = ",face
if face[0] != 0:
print " index of torial is ",k
return face
if k == 0:
w2 = w
else:
w2 = pow(w2,torials[k],z)
fac = strongfac(z,w2)
print " "
print " Could not find factors by strong probable prime test."
r = ksqrt(z)
if r * r == z:
x = r
y = r
fac = range(3)
fac[0] = x
fac[1] = y
fac[2] = 0
print " Found factors by z is square."
return fac
r = ksqrt(z)
for j in range(1,maxdiff+1):
t = r + j
d = t * t - z
ifsq = ifprobsq(d,sqlist,sqresidue)
if ifsq == 1:
r2 = ksqrt(d)
if d == r2 * r2:
x = t - r2
y = t + r2
fac = range(3)
fac[0] = x
fac[1] = y
fac[2] = 1
if j%10 == 1:
print " Found factors by ",j,"st square > z is z + a
square."
elif j%10 == 2:
print " Found factors by ",j,"nd sqare > z is z + a square."
elif j%10 == 3:
print " Found factors by ",j,"rd square > z is z + a
square."
else:
print " found factors by ",j,"th, square > z is z + a
square."
return fac
# call def diffsq(mpy,z,maxdiff,difficulty,sqlist,sqresidue) for mpy in
range(2,10000)
print " "
print " Begin to look for simple multiple of z to be difference of squares."
difficulty = [0,0]
for mpy in range(3,10000,4):
difficulty[0] = difficulty[0] + 1
fac = diffsq(mpy,z,maxdiff,difficulty,sqlist,sqresidue)
if fac[0] != 0:
return fac
difficulty[0] = difficulty[0] + 1
fac = diffsq(mpy+1,z,maxdiff,difficulty,sqlist,sqresidue)
if fac[0] != 0:
return fac
difficulty[0] = difficulty[0] + 1
fac = diffsq(mpy+2,z,maxdiff,difficulty,sqlist,sqresidue)
if fac[0] != 0:
return fac
print " "
print " Could not factor z by finding a simple multiple of z to be a
difference of squares."
print " "
print " Try to find factors by AD + BC = z method. "
fac = ABCD(z)
if fac[0] != 0:
print " Found factor by method AD + BC = z."
return fac
print " Could not factor z by AD + BC = z method. "
print " Begin final resort routine."
mpylist = setmpylist(z,plist)
# print " mpylist = ",mpylist
lenmult = len(mpylist)
mult = range(lenmult)
for k in range(lenmult):
mult[k] = 0
# zlim = log_2(2 * z)
zlim = 0
temp = 2 * z
while temp > 1:
temp = temp/2
zlim = zlim + 1
# print " zlim = ",zlim
flag = 0
difficulty = range(2)
difficulty[0] = 1
difficulty[1] = 0
while flag == 0:
difficulty[0] = difficulty[0] + 1
mult = incr(mult,zlim,mpylist)
# print " mult = ",mult," difficulty = ",difficulty
if mult[0] == -1:
print "Fermat algorithm Could not factor z "
fac = range(3)
fac[0] = 1
fac[1] = z
fac[2] = difficulty
return fac
mpy = 1
lenmult = len(mult)
for k in range(lenmult):
j = lenmult-1 - k
if mult[j] != 0:
mpy = mpy * mpylist [j] [0] **mult[j]
fac = diffsq(mpy,z,maxdiff,difficulty,sqlist,sqresidue)
if fac[0] != 0:
return fac
# mpylist contains the vector of multipliers.
# each element of mpylist is a vector of size 2.
# mpylist[k][0] is kth multiplier.
# mpylist[k][1] is approximately log_2 of kth multiplier
def setmpylist(z,plist):
flag = 0
mpylist = [[2,1]]
curr = 1
lenmpylist = 1
dex = 1
while flag == 0:
m = 2 * curr + 1
dex = dex + 1
actdex = dex
if m > z:
break
q = m
for k in range(lenmpylist):
d = mpylist [k][0]
f = mpylist [k][1]
rflag = 0
while rflag == 0:
quo = q/d
r = q - quo * d
# print " z = ",z," m = ",m," q = ",q," k = ",k," d = ",d," f =
",f," quo = ",quo," r = ",r," actdex = ",actdex
if r == 0:
q = quo
actdex = actdex - f
else:
rflag = 1
if q > 1:
mpylist.append([q,actdex])
lenmpylist = lenmpylist + 1
curr = m
return mpylist
def diffsq(mpy,z,maxdiff,difficulty,sqlist,sqresidue):
za = mpy * z
r = ksqrt(za)
if za == r * r:
x = gcd(z,r)
if x > 1:
if x < z:
y = z/x
fac = range(3)
fac[0] = x
fac[1] = y
fac[2] = difficulty
return fac
for j in range(1,maxdiff+1):
t = r + j
d = t * t - za
ifsq = ifprobsq(d,sqlist,sqresidue)
if ifsq == 1:
difficulty[1] = difficulty[1] + 1
rd = ksqrt(d)
if rd * rd == d:
x = gcd(z,t - rd)
if x > 1:
if x < z:
y = z/x
fac = range(3)
fac[0] = x
fac[1] = y
fac[2] = difficulty
if j%10 == 1:
print " Found factors by ",j,"st square > some
multiple of z, za, is za + a square."
elif j%10 == 2:
print " Found factors by ",j,"nd square > some
multiple of z, za, is za + a square."
elif j%10 == 3:
print " Found factors by ",j,"rd square > some
multiple of z, za, is za + a square."
else:
print " found factors by ",j,"th, square > some
multiple of z, za, is za + a square."
return fac
y = gcd(z,t+rd)
if y > 1:
if y < z:
x = z/y
fac = range(3)
fac[0] = x
fac[1] = y
fac[2] = difficulty
if j%10 == 1:
print " Found factors by ",j,"st square > z, za, is
za + a square."
elif j%10 == 2:
print " Found factors by ",j,"nd square > z, za, is
za + a square."
elif j%10 == 3:
print " Found factors by ",j,"rd square > z, za, is
za + a square."
else:
print " found factors by ",j,"th, square > z, za,
is za + a square."
return fac
fac = range(3)
fac[0] = 0
fac[1] = 0
fac[2] = difficulty
return fac
def incr(mult,zlim,mpylist):
# mult is a vector of exponents for the multipliers in mpylist.
# z is the integer to be factored
# zlim is the critical value for the sum of multiplier * index
# where kth multiplier is mpylist [k][0]
# and kth index is mpylist [k][1]
# function incr returns the next value of vector mult.
# mult may be thought as a number written in a variable base.
# mult[0] is the least significant and matches multiplier mpylist[0][0]
# adding one to mult would mean adding 1 to mult[0]
# unless doing so would make sum of mpylist[k][1] * mult[k] exceed zlim.
# in that case mult[0] is set to zero, and 1 is added to mult[1]
# unless doing so would make sum of multilier * index exceed zlim
# in that case, mult[0] and mult[1] is set to zero,
# and 1 is added to mult[2]
# unless . . .
# mult[0] is set to -1 to indicate that the largest possible value of mult has
been exceeded.
# mpylist[0] = 2 and all other values in mpylist are odd.
# mult[0], the index on multiplier 2, must not equal 1. It may equal zero, or
any integer > 1,
# provided it meets the zlim constraint.
smalllim = 13
lenmult = len(mult)
lim = lenmult+0
for k in range(lim):
if mult[k] == 0:
lenmult = lenmult - 1
else:
break
mult[0] = mult[0]+1
if mult[0] ==1:
mult[0] = 2
testsum = 0
for j in range(0,lenmult):
testsum = testsum + mpylist[j][1] * mult[j]
if testsum < smalllim:
mult[0] = mult[0] + (smalllim - testsum)/2
if testsum < zlim:
return mult
for k in range(1,lenmult):
mult[k] = mult[k] + 1
testsum = 0
if k > 0:
mult[k-1] = 0
for j in range(k,lenmult):
testsum = testsum + mpylist[j][1] * mult[j]
if testsum < smalllim:
mult[k] = mult[k] + (smalllim - testsum)/mpylist[k][0]
if testsum < zlim:
return mult
mult[0] = -1
return mult
def ksqrt(z):
if z == 0:
return 0
j = 1
k = z
flag = 0
while flag == 0:
j = k
k = (j + z/j)/2
if k == j + 1:
return j
if k == j:
return j
def gcd(a,b):
r1 = a
r2 = b
flag = 0
while flag == 0:
r3 = r1%r2
if r3== 0:
return r2
r1 = r2
r2 = r3
# print "Testing gcd: r1 = ",r1," r2 = ",r2
### write alternative strong(z,w) subroutine that
### also tests for factors in the case that
### w is witness for z being a weak probable prime.
def strong(z,w):
trace = 0
t = z - 1
a = 0
while t%2 == 0:
t = t/2
a = a + 1
test = pow(w,t,z)
if test ==1:
q = 1
return q
if test + 1 == z:
q = 1
return q
for j in range (1,a):
test = (test * test)%z
if test + 1 == z:
q = 1
return q
q = 0
if trace > 0:
print " found ",z," is composite, with witness = ",w
return q
def miller(z):
for p in range(2,10000):
fac = strongfac(z,p)
if fac[0] != 0:
return fac
print " Could not find factors by strong probable prime test using
witnesses < 10000."
fac[0] = 0
fac[1] = 0
fac[2] = 0
return fac
def strongfac(z,w):
trace = 0
t = z - 1
a = 0
while t%2 == 0:
t = t/2
a = a + 1
test = pow(w,t,z)
if test ==1:
q = 1
fac = [1,z,0]
return fac
else:
x = gcd(test-1,z)
if x > 1:
print " "
print " In strongfac "
print " Found1: x = ",x
if x < z:
y = z/x
print " Found1: y = ",y
face = [x,y,0]
print " face = ",face
face[0] = x
face[1] = y
face[2] = 0
print " Found1: factor by strong probable prime test, using
witnes ",w," ."
print " Found1: factors are ",face[0],face[1]
print " Found1: returning ",face," to fermat routine."
return face
if test + 1 == z:
q = 1
fac = [1,z,0]
return fac
for j in range (1,a):
test = (test * test)%z
test1 = test + 1
if test1 == z:
q = 1
fac = [1,z,j]
return fac
else:
x = gcd(test1,z)
if x > 1:
if x < z:
print " "
print " In Strongfac "
y = z/x
fac = [x,y,j]
print " Foundj factor by strong probable prime test, using
witness,",w," ."
print " Returning ",fac," to fermat routine."
return fac
q = 0
fac = [0,0,0]
return fac
def ifprime(z,plist):
for j in range(1,1230):
w = plist[j]
if w == z:
q = 1
return q
q = strong(z,w)
if q == 0:
return q
return 1
def setplist(pmap):
plist = range(1230)
n = 0
if pmap[3] != 3:
print " pmap has not been calculated! "
return
for j in range(2,10000):
if pmap[j] != -1:
n = n + 1
plist[n] = j
plist[0] = n
return(plist)
def setpmap():
pmap=range(10001)
pmap[0] = 0
pmap[1] = 0
pmap[2] = 2
for j in range(3,10000,2):
pmap[j] = j
pmap[j+1] = -1
for k in range(3,100,2):
if pmap[k] != -1:
jstart = k*k
jend = 10000
jinc = 2 * k
for j in range(jstart,jend,jinc):
pmap[j] = -1
return(pmap)
_______________________________________________
Tutor maillist - Tutor@python.org
http://mail.python.org/mailman/listinfo/tutor
