Python & Sage debugging
Posted: Fri Feb 06, 2009 2:09 am
I downloaded some code, polynomial-roots.tar, from here to build some images from solving polynomials (those things you did in your high school/middle school math classes). The code is about 3 years old, and therefore a little old.
So far I've downloaded sage, psyco, and SciPy. and tweaked the 'genroots.sage' file from the tar ball.
The only revisions from the tarball that I made are lines 1 and 15.
1)line 1 was modified from
to
2)line 15 was modified from
to
I
The reason for the first change is to update sage-python to my current version, sage. The second change is to move from the old python to the newer python I have on my computer. I thought these were fairly reasonable to change, and shouldn't impact the code much if at all.
The only trouble trouble is when I execute './genroots.sage' I get the error is:
I'm doing this for fun primarily. I tried making some matlab code to build the image found on the first link on my own. I ran into some coding trouble as usual, and was left with rudimentary 'bottle cap' images. So now I'm trying to use this code to build the images.
So far I've downloaded sage, psyco, and SciPy. and tweaked the 'genroots.sage' file from the tar ball.
Code: Select all
#!/usr/bin/env sage
# A version which stores roots in files, organized
# by degree of minimal polynomial.
# ./genroots.sage
# Took about 10 minutes for d=5, rng=5 on jdc.
# 176622 polynomials, 14 megabytes of output.
from xcombine import xcombine
from sage.all import *
import sys
sys.path.append("/usr/lib/python2.5/site-packages")
from scipy import *
# psyco only helps a little. Need to install custom version
# for sage, since python version is different.
import psyco
psyco.full()
def sht(x):
'''Return a short string representation of a float.'''
# %.4f includes exactly four digits to the right of the
# decimal point, i.e. it specifies the absolute precision,
# which is what we care about for plotting. The rstrip
# removes unneeded trailing 0's and decimal point.
# (rstrip('0.') changes '0.0000' to ''.)
# (%.4g gives four significant digits, which is silly
# for really small numbers.)
return ('%.4f' % (x,)).rstrip('0').rstrip('.')
def shtcmpl(x):
'''Return a short string representation of a complex number.'''
return ('%.4f' % (x.real,)).rstrip('0').rstrip('.') + \
('%+.4f' % (x.imag,)).rstrip('0').rstrip('.') + 'j'
def rts(d=6, rng=5):
if d < 1:
return
fd = {}
for deg in range(1,d+1):
for height in range(1, rng+1):
fd[(deg, height)] = open('roots%d-%d.txt' % (deg, height), 'w')
# Put 0 at end, so low degree polynomials come after
# high degree polynomials.
r = range(-rng,0)+range(1,rng+1)+[0]
# Use just positive constant coefficient, since negating
# polynomial doesn't change roots, and if constant coeff is zero,
# polynomial factors.
rpos = xrange(1, rng+1)
# rnn used for linear factor; see comment below.
rnn = range(1, rng+1) + [0]
R=PolynomialRing(ZZ)
for coeffs in xcombine(*((r,)*(d-1)+(rnn,rpos))):
# Replacing x with -x just reflects the roots, so we try to
# save time by not generating both. We use rnn to ensure that
# linear coefficient coeffs[-2] is non-negative. If it is zero,
# we can force the cubic coefficient coeffs[-4] to be non-negative.
# If *that* is zero, we can force the fifth degree coefficient to
# be non-zero. We stop there, so we lose some efficiency for d >= 7.
if d >= 3 and coeffs[-2] == 0 and coeffs[-4] < 0:
continue
if d >= 5 and coeffs[-2] == 0 and coeffs[-4] == 0 and coeffs[-6] < 0:
continue
rcoeffs = list(coeffs)
rcoeffs.reverse()
order = d
while order >= 0 and rcoeffs[order] == 0:
order -= 1
# If the polynomial is constant, it has no roots.
# Constant polynomials give an error for .is_irreducible(),
# so that's why we eliminate them here.
if order <= 0:
continue
# If coeffs have a common factor, then this polynomial is redundant:
if reduce(gcd, coeffs) != 1:
continue
# Could also look at x |--> 1/x, plus multiply by x^degree,
# which is just the same as reversing the coefficients.
# This takes irreducibles to irreducibles, giving
# giving another factor of two symmetry. Using this would
# save time computing the roots, but still produce the same number
# of dots. Only deal with absolute values of coefficients here,
# since I haven't figured out how this interacts with the other
# two symmetries; so I suspect there are some redundant polynomials
# still generated.
abscoeffs = map(abs, coeffs)
absrcoeffs = map(abs, rcoeffs)
num = (order+1)//2
left = abscoeffs[-(order+1):-(order+1)+num]
right = absrcoeffs[:num]
if right < left:
continue
#print abscoeffs, ':', left, '+', right
# If we didn't exclude polynomials with zero constant term, we'd
# have to make sure we didn't try to find the roots of the zero polynomial.
# Scipy is roughly 5 times faster than sage at finding roots:
p = poly1d(coeffs)
#rts = p2.complex_roots()
# For deg=4, rng=4, all roots have abs <= 5.
# In my tests, since in most cases we need to compute the
# roots anyways, it turns out to be faster to compute all
# roots and look for integer roots than to test integers
# in range(-5,6). Surprising to me.
#if order > 1:
# skip = False
# for x in range(-5,6):
# if p(x) == 0:
# skip = True
# break
# if skip: continue
rts = p.r
#for rt in rts:
# maxabs = max(maxabs, abs(rt))
# This succeeded:
#assert(len(rts) == len(rts2))
#for r in rts:
# found = False
# for r2 in rts:
# if abs(r-r2) < 1e-9:
# found = True
# break
# assert(found)
# Search for linear factors; this gets rid of some duplicates.
# Overall, it takes more time than it saves, but I really want
# to associate each root with the polynomial of lowest degree
# that it is a root of.
# If we find an integer root, then polynomial factors over Z[x].
# But, how do I know that the factors have coeffs in -rng..rng????
# We don't know, but I still think it makes sense to eliminate
# all reducible polynomials. This changes the plot (in theory),
# but is "correct" in some sense.
if order > 2:
skip = False
for rt in rts:
if -1e-6 < rt.imag < 1e-6:
closeint = int(round(rt.real))
# This evaluates p(closeint) using integer arithmetic.
# The first condition just saves some time.
if abs(closeint-rt) < 1e-6 and p(closeint) == 0:
skip = True
break
if skip:
# This succeeds:
#assert(not p2.is_irreducible())
continue
# Look for quadratic factors?
# Eisenstein test?
# Since we only use this to check irreducibility, it
# doesn't actually matter whether we reverse the coeffs.
p2 = R(rcoeffs)
# Just calling is_irr on each poly doubles the run time.
# For deg=4, rng=4, 10450 of 12588 polys are irreducible:
if not p2.is_irreducible():
continue
height = max(abscoeffs)
fdoh = fd[(order, height)]
for c in coeffs[-(order+1):]:
fdoh.write(sht(c)+' ')
fdoh.write('; ')
for rt in rts:
fdoh.write(shtcmpl(rt)+' ')
fdoh.write('\n')
fdoh.flush()
# Include 0 here, since polynomials with zero constant term excluded:
fd[(1,1)].write('1 0 ; 0+0j \n')
if __name__ == '__main__':
rts()
1)line 1 was modified from
Code: Select all
#!/usr/bin/env sage-pythonCode: Select all
#!/usr/bin/env sageCode: Select all
sys.path.append("/usr/lib/python2.4/site-packages")Code: Select all
sys.path.append("/usr/lib/python2.5/site-packages")The reason for the first change is to update sage-python to my current version, sage. The second change is to move from the old python to the newer python I have on my computer. I thought these were fairly reasonable to change, and shouldn't impact the code much if at all.
The only trouble trouble is when I execute './genroots.sage' I get the error is:
I have xcombine.py in the same folder as genroots.sage. I have no clue why it is having an import error. Any help is appreciated.Traceback (most recent call last):
File "./genroots.py", line 12, in <module>
from xcombine import xcombine
ImportError: cannot import name xcombine
I'm doing this for fun primarily. I tried making some matlab code to build the image found on the first link on my own. I ran into some coding trouble as usual, and was left with rudimentary 'bottle cap' images. So now I'm trying to use this code to build the images.