(Update: Look two posts down, where I have updated code that is much faster, probably on par with the Mathematica code that Andrew posted.)
I don't recall what the exact term is for the matrix which this code computes. I call it the FMA, which I believe stands roughly for "flow matrix". You see, "fma" was the name of a variable in some code which Andrew provided two years ago:
http://math.eretrandre.org/tetrationforu...147#pid147
I don't remember what the matrix is called, but I know what it represents. The columns (or rows, depending on transposition) are the coefficients of Lagrange interpolation polynomials in n, and each polynomial interpolates a coefficient of the power series of the nth iterate of the desired function, provided it's centered at a parabolic fixed point. For example, I have been using it to calculate for exp(z)-1.
Well, the code apparently runs much slower in SAGE than it does in Mathematica. Andrew converted the code directly, it seems, so efficiencies in Mathematica are apparently not conserved in going to SAGE. Fear not, this merely means we must find an equivalent form that is efficient in SAGE, or at least relatively so.
If you look at the code, there are three subfunctions that Andrew wrote. The first calculates the coefficients of the integer iterates of a function centered at its parabolic fixed point (can't recall if the fixed point needs to be 0, but the change is trivial if so). The second two functions create the interpolation polynomials.
Well, the first change I made was in the parabolic_idm function. Rather than computing the taylor series of ser(expr), I changed it to compute expr(ser). It turns out that internal optimizations make this latter form much faster, about four times faster for size=100, five times faster for size=150, and generally with the speed increase getting slightly better with increasing size.
Using the matrix approach that Gottfried helped me figure out, I was able to reduce the second and third function calls to a matrix multiplication, which takes only a matter of seconds for matrices this small. These latter two function calls were taking almost as much time as the first, so reducing this to a few seconds helps dramatically.
Without further ago, here is the code, which runs in SAGE 4.0.1 and 4.1, but is probably backwards compatible as well. (Depending on how you convert Andrew's fma list to a matrix, it might be a transposition of the one I calculate below. So don't freak out if you run into this when trying to compare.)
Note: the code gives a deprecation warning when calculating the Taylor series, but I'm not sure about the best way to avoid it. (I'll post an update if I can work around it without hurting performance.)
PS: If anyone tries this code and it doesn't work, let me know. If anyone can suggest further ways to optimize, I'd like to know, as learning to optimize SAGE code can be helpful for me in general.
I don't recall what the exact term is for the matrix which this code computes. I call it the FMA, which I believe stands roughly for "flow matrix". You see, "fma" was the name of a variable in some code which Andrew provided two years ago:
http://math.eretrandre.org/tetrationforu...147#pid147
I don't remember what the matrix is called, but I know what it represents. The columns (or rows, depending on transposition) are the coefficients of Lagrange interpolation polynomials in n, and each polynomial interpolates a coefficient of the power series of the nth iterate of the desired function, provided it's centered at a parabolic fixed point. For example, I have been using it to calculate for exp(z)-1.
Well, the code apparently runs much slower in SAGE than it does in Mathematica. Andrew converted the code directly, it seems, so efficiencies in Mathematica are apparently not conserved in going to SAGE. Fear not, this merely means we must find an equivalent form that is efficient in SAGE, or at least relatively so.
If you look at the code, there are three subfunctions that Andrew wrote. The first calculates the coefficients of the integer iterates of a function centered at its parabolic fixed point (can't recall if the fixed point needs to be 0, but the change is trivial if so). The second two functions create the interpolation polynomials.
Well, the first change I made was in the parabolic_idm function. Rather than computing the taylor series of ser(expr), I changed it to compute expr(ser). It turns out that internal optimizations make this latter form much faster, about four times faster for size=100, five times faster for size=150, and generally with the speed increase getting slightly better with increasing size.
Using the matrix approach that Gottfried helped me figure out, I was able to reduce the second and third function calls to a matrix multiplication, which takes only a matter of seconds for matrices this small. These latter two function calls were taking almost as much time as the first, so reducing this to a few seconds helps dramatically.
Without further ago, here is the code, which runs in SAGE 4.0.1 and 4.1, but is probably backwards compatible as well. (Depending on how you convert Andrew's fma list to a matrix, it might be a transposition of the one I calculate below. So don't freak out if you run into this when trying to compare.)
Code:
# Define the custom IDM:
def parabolic_idm(expr, x, size):
ret = [];
ser = x;
for i in xrange(size):
ser = taylor(expr(ser), x, 0, size);
cof = ser.coeffs(x);
ret.append([0] + [cof[int(k)][0] for k in xrange(len(cof))]);
top = [0, 1] + [0 for i in xrange(size-1)];
return [top] + 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 = 10
time mtx = parabolic_idm(exp(x)-1, x, size)
time fma = LagrangePolynomialMatrix(mtx, size)Note: the code gives a deprecation warning when calculating the Taylor series, but I'm not sure about the best way to avoid it. (I'll post an update if I can work around it without hurting performance.)
PS: If anyone tries this code and it doesn't work, let me know. If anyone can suggest further ways to optimize, I'd like to know, as learning to optimize SAGE code can be helpful for me in general.
~ Jay Daniel Fox
