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.