SAGE code for computing flow matrix for exp(z)-1

[jaydfox]

Alas, it seems that Maxima and/or Lisp seems to crash for a large size, or after repeated calculations, due to small, hard-coded memory limit in the Maxima installation. According to discussions on the sage-support google group, it might require a recompilation of Maxima!

Well, I know ways to work around this, as there are other ways to compute the taylor series of composed functions, such as through an alternative library like gp/pari, but I don't have time to work on a workaround today. Sorry, folks!

Okay, so I decided to just calculate the taylor series of the iterations "by hand", using matrix math. It turns out to be way faster. I am now getting times comparable to what Andrew was getting in Mathematica. Hooray!

Code:

# Define the custom IDM:
def parabolic_idm(field, expr, x, size):
    mulm = matrix(field, size+1, size+1);
    abel = matrix(field, size+1, size+1);
    ret = matrix(field, size+1, size+1);
    cof = taylor(expr, x, 0, size).coeffs();
    for rrow in xrange(size+1):
        for ccol in xrange(len(cof)):
            col2 = cof[ccol][1]+rrow;
            if col2 < size+1:
                mulm[rrow, col2] = field(cof[ccol][0]);
    vec = vector(field, size+1);
    vec[0] = field(1);
    abel[0] = vec;
    for ii in xrange(size):
        vec = vec * mulm;
        abel[ii+1] = vec;
    vec = vector(field, size+1);
    vec[1] = field(1);
    ret[0] = vec;
    for ii in xrange(size):
        vec = vec * abel;
        ret[ii+1] = vec;
    return ret;

# Create a Lagrange Interpolation polynomial for each
# coefficient of the iteration function. Note that
# multiplying on the left by a row-vector of powers
# of n (where n is the desired iteration) will give
# a row vector which represents the coefficients of
# a power series in z. Alternatively, multiplying
# on the right by a column vector of powers of z
# will give the coefficients of a power series in n,
# which allows direct calculation of the nth iterate,
# centered at z.
def LagrangePolynomialMatrix(pidm, size):
    XI = matrix(ZZ, size+1);
    XI[0,0] = factorial(size);
    for ccol in xrange(1, size+1):
        XI[0, ccol] = (1-ccol)/ccol * XI[0, ccol-1]
    for rrow in xrange(1,size+1):
        for ccol in xrange(rrow, size+1):
            XI[rrow, ccol] =  (XI[rrow-1, ccol-1] + (1-ccol) * XI[rrow, ccol-1])/ccol
    PI = matrix(ZZ, size+1);
    PI[0,0] = 1;
    for rrow in xrange(1, size+1):
        for ccol in xrange(rrow+1):
            PI[rrow, ccol] = PI[rrow-1, ccol-1] - PI[rrow-1, ccol];
    return (XI*PI)*matrix(QQ, pidm)/factorial(size);

# test data: change or delete as needed
size = 100
time mtx = parabolic_idm(QQ, exp(x)-1, x, size)
time fma = LagrangePolynomialMatrix(mtx, size)

Update: It occurs to me that memory can be saved by eliminating the "mulm" variable and repurposing the "ret" variable. Thus:

Code:

# Define the custom IDM:
def parabolic_idm(field, expr, x, size):
    abel = matrix(field, size+1, size+1);
    ret = matrix(field, size+1, size+1);
    cof = taylor(expr, x, 0, size).coeffs();
    for rrow in xrange(size+1):
        for ccol in xrange(len(cof)):
            col2 = cof[ccol][1]+rrow;
            if col2 < size+1:
                ret[rrow, col2] = field(cof[ccol][0]);
    vec = vector(field, size+1);
    vec[0] = field(1);
    abel[0] = vec;
    for ii in xrange(size):
        vec = vec * ret;
        abel[ii+1] = vec;
    vec = vector(field, size+1);
    vec[1] = field(1);
    ret[0] = vec;
    for ii in xrange(size):
        vec = vec * abel;
        ret[ii+1] = vec;
    return ret;