(08/19/2009, 09:17 PM)jaydfox Wrote: 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!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!

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!

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;

~ Jay Daniel Fox