Designing a Tetration Library

[bo198214]

Andrew, I am looking for ways for efficiently working together.
I would like to provide the regular iteration things, as I have the most experience with them here on the forum.

I think there is a lot to discuss, perhaps mainly about naming, dividing into packages and how to handle bugs of sage, etc. So
1. we need a communication channel. Perhaps skype, googletalk, ICQ, YIM, MSN, etc? and roughly times when we are available for discussing the tetration library. I am per default available via skype and google talk.

Then we have to keep track of a common code
2. we need a common repository. I worked a lot with subversion, however it is said that the repository git developed by Torwalds for the Linux kernel is more sophisticated. Has anyone experiences with it?
We also need a server to run it.

So then I would like to start the naming discussion, with how the tetration library shall be called? I mean it surely shall not be just about tetration. It should include all the hyperoperations of different kind and also the methods for non-integer iteration of functions.

So I would propose the name HyperOp. However I am not familiar with typical naming in Sage (to be honest i am an absolute novice to sage).

[andydude]

Well, I've included an algorithm for hyperbolic iteration (the more general "complement" of parabolic iteration, where \( f(0) = 0 \) and \( f'(0) \ne 1 \)). However, I was hoping to make it general enough to include an implementation of matrix_exp and matrix_log at the same time, but I do not have working alorithms for these yet (but thanks to Gottfried we do, just in a different language). The matrix power code for hyperbolic iteration is very specific to triangular matrices, so I'm not sure how well it could be generalized to arbitrary matrices.

I found that the Sage built-in exp will only work on numeric matrices, and that the .left_eigenvectors() and .right_eigenvectors() methods are only defined on numeric matrices as well. The .eigenvalues() method, however, IS defined on symbolic matrices, but we don't need the eigenvalues (which are \( (1, f_1, f_1^2, f_1^3, \cdots) \)) so this is of no help. One of the purposes of this library is to investigate the symbolic case, which means I had to implement my own eigensystem decomposition algorithm. I checked that it produces the same result as it should in Mathematica, so at least we are being consistent.

Here are some examples of some new functions:

Code:

sage: from hyper.all import *
sage: t = var('t')
sage: hyperbolic_flow(h_poly(x), t, x, 0, 2)
x^2*C1^(t - 1)*(C1^t - 1)*C2/(C1 - 1) + x*C1^t
sage: hyperbolic_flow_coeffs(h_poly(x), t, x, 0, 2)

[0,
C1^t,
C1^(t - 1)*(C1^t - 1)*C2/(C1 - 1)]
sage: c = var('c')
sage: hyperbolic_flow_coeffs(e^(c*x) - 1, t, x, 0, 3)

[0,
c^t,
c^(t + 1)*(c^t - 1)/(2*(c - 1)),
c^(t + 2)*(c^t - 1)*(c^(t + 1) + 2*c^t - 2*c - 1)/(6*(c - 1)^2*(c + 1))]

which is what is to be expected.

I would very much like to work together! It would also help the project be more robust with more eyes making sure the code is right

I have a crappy server with a dynamic IP, but I think we would need one with a static IP to be more reliable. "HyperOp" is a good name for a library, but I was thinking "HyperSage". What do you think?

Now that both basic parabolic and hyperbolic methods are fleshed out, I'm not going to work much more on them. There is also a regular iteration sub-module, but the only stuff that would go here is Daniel Geisler's summation method, which I don't understand very well. Since most of regular iteration is working (parabolic/hyperbolic), I'm going to start working on the natural iteration sub-module.

Andrew Robbins

[Gottfried]

I would very much like to work together! It would also help the project be more robust with more eyes making sure the code is right

I have a crappy server with a dynamic IP, but I think we would need one with a static IP to be more reliable.

I could provide my static IP at my office-pc. I have Windows running, can implement html,ftp via xitami-server and news/mail via hamster. Also I could implement some specific cgi-program. (If there are other free, reliable components running on Win98/128 MB I could install them, too)

[bo198214]

Well, I've included an algorithm for hyperbolic iteration (the more general "complement" of parabolic iteration,

At least one person works at the tetration library ... I didnt even found the time to look through your code, but it will happen
But (!) I now have Sage actually working on my new Linux machine!

Though I found out that Sage has only support for truncated powerseries.
In Maple you could define a powerseries by a function that returns for each index the coefficient, and then apply various operations on powerseries like exponentation, composition, multiplication, addition, etc (though they were sometimes very slow ...).
Cool would also be a lazy evaluation ala Haskell, which unsurprisingly has a power series module, lazy evaluation seems just to be made for power series operations.

The matrix power code for hyperbolic iteration is very specific to triangular matrices, so I'm not sure how well it could be generalized to arbitrary matrices.

The method is generally known as Matrix function and better described at Jordan matrix. In our case the function is just the power \( x^t \).
In the matrix function article there is also mentioned a method via the Cauchy integral. I wonder whether this is a link (equality?) between Kouznetsov's method and the diagonalization method.

More useful for implementing the matrix power - without needing to perform an actual diagonalization - is Aldrovandi's explanation featured by Gottfried here.

In Maple everything is already there. You just can apply any function to matrices and you have even a specialized function Matrixpower.

I found that the Sage built-in exp will only work on numeric matrices, and that the .left_eigenvectors() and .right_eigenvectors() methods are only defined on numeric matrices as well.

Usually you have symbolic eigenvalues only up to 3 or 4. Because they are the solutions of the characteristic polynomial. I guess there are only a few occasions where symbolic eigenvalues make sense.

I would very much like to work together! It would also help the project be more robust with more eyes making sure the code is right

I have a look at it soon. Btw I would really be keen on having a public git project here or here or here. How about registering one?

[bo198214]

Ok, I now created a public git repository, having Andrew's sources checked in.
Its on github, you can find it here.

The public git url is:
git://github.com/bo198214/hyperops.git

So what do you have to do:

1. Install git if not done already,
   
2. Read about git basics here if not done already. And do those
   
   Code:

   $ git config --global user.name "Your Name Comes Here"
$ git config --global user.email you@yourdomain.example.com
3. go to the directory where you store your projects and then do
   
   Code:

   git clone git://github.com/bo198214/hyperops.git
   This creates the directory 'hyperops' with all the sources that Andrew just provided contained.
   
4. You can locally work on it with 'git commit', 'git branch' , 'git merge' etc.
   
5. Give me your public ssh-key if you want to make changes/contribute to the code. It is found in: .ssh/id_rsa.pub or .ssh/id_dsa.pub. If unsure perhaps first read this
   
6. I ll add it to the repository and then you can write back your local changes via 'git push'.
   
7. thats it, have fun
   

Edit: changed hyperop to hyperops as this obeys the Sage recommendations

[bo198214]

Ok, here are my results of learning python and sage via implementing a proper powerseries module "hyperops.powerseries". I didnt dare to already incorporate it into our git repository before consilidating and synchronizing us about it. So I just attach the single file.

I also have a Sage question: How does one incorporate own packages into the sage hierarchy so that they are recognized from whatever path one starts sage?

Ok, now the module description:
The main difference to the existing powerseries package is that we deal with untruncated powerseries, which are evaluated as soon as we request it.
A powerseries is mainly given as a function that returns the coefficient for the index, i.e. what we write often as \( f_n \), the \( n \)-th coefficient of the powerseries \( f \).
Maple has such a package however its awfully slow if one want to compute for example hyperbolic iteration, because it lacks caching. This is a simple mechanism that just stores a once computed value for subsequent access. It is life-saving in the case of recursive computations as we need to do.
There are also other ugly matters with the Maple powerseries, for example that they always need a name. When created anonymously one get errors.

So thatswhy I was burning to create a proper powerseries package without all those drawbacks. And found that python is an absolute nice language to perform such tasks including functional programming style. And before long explanations I just give the example code, which is also included in the documentation section of the module.

Our beloved regular iteration is of course contained (hyperbolic and parabolic gathered in one routine deciding by the first coefficient which one to chose) and demonstrated here with the decremented exponential.

Code:

sage: from hyperops.powerseries import PowerSeriesI
sage: #Predefined PowerSeries                                                        
sage: expps = PowerSeriesI().Exp()
sage: expps.poly(10,x)
x^9/362880 + x^8/40320 + x^7/5040 + x^6/720 + x^5/120 + x^4/24 + x^3/6 + x^2/2 + x + 1
sage: expps
[1, 1, 1/2, 1/6, 1/24, 1/120, 1/720, 1/5040, 1/40320, 1/362880, 1/3628800, ...]
sage: PowerSeriesI().Zero()
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ...]
sage: PowerSeriesI().One()
[1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ...]
sage: PowerSeriesI().Identity()
[0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ...]
sage: #finite powerseries                                                            
sage: p = PowerSeriesI([1,2,3,4])
sage: p.poly(10,x)
4*x^3 + 3*x^2 + 2*x + 1
sage: p
[1, 2, 3, 4, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ...]
sage: one = PowerSeriesI([1])
sage: id = PowerSeriesI([0,1])

sage: #power series are just functions that map the index to the coefficient          
sage: expps(30)
1/265252859812191058636308480000000

sage: #power series operations                                                        
sage: p+p
[2, 4, 6, 8, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ...]
sage: p-p
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ...]
sage: p*p
[1, 4, 10, 20, 25, 24, 16, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ...]
sage: #/ not yet implemented                                                          
sage: p**2
[1, 4, 10, 20, 25, 24, 16, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ...]
sage: #composition only works for coefficient 0 being 0 in the second operand        
sage: dexp = expps - one
sage: expps.compose(dexp)
[1, 1, 1, 5/6, 5/8, 13/30, 203/720, 877/5040, 23/224, 1007/17280, ...]

sage: #we come into interesting regions ...                                          
sage: dexp.compose(dexp)
[0, 1, 1, 5/6, 5/8, 13/30, 203/720, 877/5040, 23/224, 1007/17280, ...]
sage: dexp.iterate(2)
[0, 1, 1, 5/6, 5/8, 13/30, 203/720, 877/5040, 23/224, 1007/17280, ...]
sage: hdexp = dexp.iterate(1/2)
sage: hdexp
[0, 1, 1/4, 1/48, 0, 1/3840, -7/92160, 1/645120, 53/3440640, -281/30965760, ...]
sage: #verifying that shoudl be Zero                                                  
sage: hdexp.iterate(2) - dexp
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ...]

sage: #symbolic (parabolic) iteration                                                
sage: dexp.iterate(x)
[0, 1, x/2, 5*(x - 1)*x/12 - (x - 2)*x/6, ...]
sage: q = dexp.iterate(1/x).iterate(x)
sage: q(3)
(5*(1/x - 1)/(6*x) - (1/x - 2)/(3*x) + 1/(2*x^2))*(x - 1)*x/2 - (5*(1/x - 1)/(12*x) - (1/x - 2)/(6*x))*(x - 2)*x
sage: #simiplify and compare                                                          
sage: expand(q(3))
1/6
sage: dexp(3)
1/6

sage: #you can initialize power series with arbitrary functions on natural numbers    
sage: #for example the power series of sqrt(2)^x can be given as                      
sage: bsrt = PowerSeriesI(lambda n: diff(sqrt(2)^x,x,n)(x=0)/factorial(n))

sage: #making the first coefficient 0 to get the decremented exponential  
sage: def coeff(n):
sage:     if n == 0:
sage:         return 0
sage:     else:
sage:         return bsrt(n)

sage: dbsrt = PowerSeriesI(coeff)

sage: #and now starting hyperbolic iteration                                          
sage: dbsrt2 = dbsrt.iterate(x).iterate(1/x)

sage: #Sage is not able to simplify                                                  
sage: simplify(dbsrt2(3))
...

sage: #but numerically we can verify equality                                        
sage: RR(dbsrt2(3)(x=0.73)-dbsrt(3))
-8.67361737988404e-19

The pure command file (without the outputs) is appended. You can just paste into sage.

[bo198214]

Wah I found this site sage notebook. There you can online interactively work with sage on the server and publish your notebooks so that everyone can reproduce your results.

I just took the chance, invested some work and made the powerseries implementation into an online notebook accessible for everyone:

Try it out here!

PS: Naturally the site is sometimes somewhat slow so dont be impatient

EDIT: Oh I just see that you cant test the script as visitor, you have to create an own account and I have to add you into the collaborators list ... hm

[andydude]

Wow! Nice power series package! I just looked at the .py file, and I noticed a few things:
(1) Several methods are missing underscores (for example _repr_ should be __repr__)
(2) Your inverse() method is empty, but I have implemented series_inverse() already, so we can use this.

I'll let you know when I have "git" installed.

Andrew Robbins

[andydude]

@Henryk
I just read your implementation for parabolic and hyperbolic iteration, and I recognize that you are using Jabotinsky's formula for the parabolic case, but for the hyperbolic case, you seem to be using a method I've never seen before. Embedded in your implementation is the implicit recurrence equation:

\(
f^{\circ t}(x) = x f_1^t + \sum_{n=2}^{\infty} \frac{x^n}{n!} \sum_{k=1}^{n-1}
\frac{f^{(k)} ((f^{\circ t})^k)^{(n)} - (f^{\circ t})^{(k)}(f^k)^{(n)}}{(f_1^n - f_1)k!}
\)

Where did you find this? Is this your own formula? Have you posted this before?

Andrew Robbins

[bo198214]

but for the hyperbolic case, you seem to be using a method I've never seen before. Embedded in your implementation is the implicit recurrence equation:

I described it here. The principle \( f^{\circ t} \circ f=f\circ f^{\circ t} \) with \( {f^{\circ t}}_1={f_1}^t \) is common in the regular iteration literature, however the formula is perhaps too dirty that I have seen it explicitely mentioned.

(1) Several methods are missing underscores (for example _repr_ should be __repr__)

You can read it in the sage package guidelines somewhere that you should use the sage specific _repr_ so that the user has the possibility to redefine __repr__

(2) Your inverse() method is empty, but I have implemented series_inverse() already, so we can use this.

Yes, the nice thing about the (hyp,par) iterate method is, that it is universally also applicable for exponent (-1). So the inverse function is just a wrapper for iterate(-1).

However perhaps in the future I will also include the direct formulas, for example to take the compositional root, you can easily derive an iterative formula solving \( f^{\circ n} = g \) for \( f \). Though the result should be the same as if you just do \( g^{\circ 1/n} \) vie the regular iteration. I also plan to implement the power \( f^{x} \) for arbitrary \( x \) and not just for integers as it is now.

Originally I wanted to use "^" via __and__ for the iteration, like I used "**" via __pow__ for the powers. Do you have an idea why this is not recognized? And also the question about how to put the hyperops package into sage that it is automatically recognized is still open for me.