Page 1 of 1

Python & Sage debugging

Posted: Fri Feb 06, 2009 2:09 am
by afbase
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.

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()
 
 
The only revisions from the tarball that I made are lines 1 and 15.
1)line 1 was modified from

Code: Select all

#!/usr/bin/env sage-python
to

Code: Select all

#!/usr/bin/env sage
2)line 15 was modified from

Code: Select all

sys.path.append("/usr/lib/python2.4/site-packages")
to

Code: Select all

sys.path.append("/usr/lib/python2.5/site-packages")
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:
Traceback (most recent call last):
File "./genroots.py", line 12, in <module>
from xcombine import xcombine
ImportError: cannot import name xcombine
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. :D

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.

Re: Python & Sage debugging

Posted: Wed Feb 11, 2009 3:26 pm
by visu
Hiya, best asking on a python forum tbh, try changing the import

Code: Select all

from xcombine import xcombine
to either

Code: Select all

from xcombine import * 
#or
import xcombine