Thread Rating:
  • 0 Vote(s) - 0 Average
  • 1
  • 2
  • 3
  • 4
  • 5
Designing a Tetration Library
#1
I aim to address many things in these posts (I have 3 more planned):
  1. Definitions of accuracy and precision
  2. Keeping track of accuracy and precision
  3. Algorithms for natural super-logarithms
  4. Algorithms for regular tetration
  5. General algorithms for tetration
  6. General algorithms for hyper-operations

To do this, I am going to bring together several of Jay Fox's posts on accuracy and precision, since I think I understand them enough to shed some light on the issue.

Jay Fox has posted several threads about improving natural tetration (as we call it), which investigate acceleration (overview (slog/sexp), details (slog)) and accuracy & precision of the natural super-logarithm. I'm currently building a library based on these techniques, which will allow one to ask: "What is e^^pi to 14 significant figures?" (which Jay has already found to 25 sig.fig.) and since Jay Fox has found an upper bound for the possible error, we will know for sure that the result will be accurate to that many decimal places.

How do we do this? First of all, we need a general algorithm for computing tetrate(b, h), superroot(h, z), and the superlog(b, z). The natural method applies to superlog() first, and tetrate() second (you must invert the function somehow), and the regular method applies to tetrate() first, then superlog() second (again, because you must invert it). I will focus on the general algorithms first, in other words, where this inversion takes place, rather than focusing on the method and saying "then invert it". This will be the context in which all other functions will fit in.

Some of the algorithms involve accuracy_of() and precision_of(), which imply a language with arbitrary-precision. These can be, or are, implemented in Sage, Maple, Mathematica, and GMP. But first, I should define these terms a little before continuing. When I say "precision" I mean the number of significant digits. When I say "accuracy" I mean the negative of the least significant digit (base 10). So a number like 3.14 will have "accuracy" 2 and "precision" 3. I will define these in more depth in the next post.

Here is the pseudo-code which is meant to be read, but not executed:

First, a basic framework for tetrate()
Code:
tetrate(base, h) :=
  # try to match the accuracy of 'h'
  # accuracy could also be d=15 by default
  tetrate(base, h, accuracy_of(h))

tetrate(base, h, d) :=

  if h is integer
  
  then: if base is in {-infinity, -1, 0, 1}
        then:
          # use "Limit" method (mainly for b=0)
          limit_tetrate(base, h)
        else:
          # use "Naive" method
          naive_tetrate(base, h)
        
  else: case base == -1: (undefined) or -1 ??
             base == 0: (undefined)
             base == 1: if h > -1
                     then: return 1
                     else: (undefined)
  
             base == infinity:
                case h < -2: (undefined)
                     -2<h<0: return 0
                     h > 0: return infinity

             # prefer "Regular" method for tetrate()
             # this could be Re(base) <= 1 if desired
             # but the series may not converge
             |log(H(base))| <= 1:
                # use "Regular" method
                n = regular_tetrate_size(d)
                regular_tetrate(base, h, n)
  
             Re(base) > 1:
                # use "Natural" method
                [n_sol, n_est] := natural_superlog_size(d)
                natural_tetrate(base, h, n_sol, n_est)
                
             otherwise: (undefined)
end

That was easy! but how are all the sub-functions like natural_tetrate() defined? I will cover these in a future post. Also, the reason why the domains of regular and natural tetration are limited above can be explained in this post.

Second, a basic framework for superlog()
Code:
superlog(base, z) :=
  # try to match the accuracy of 'z'
  # accuracy could also be d=15 by default
  superlog(base, z, accuracy_of(z))

superlog(base, z, d) :=

  [h, t] := ceiling_superlog(base, z)
  
  if t == 0
  
  then: return h

  else: case base == -1: (undefined)
             base == 0: (undefined)
             base == 1: case z > 1:  return infinity
                             z == 1: return 0
                             z < 1:  return -1
  
             base == -infinity: (undefined)
             base == infinity: case z > 0:  return 0
                                    z == 0: return -1
                                    z < 0:  return -2
  
             # prefer "Natural" method for superlog
             Re(base) > 1:
                  # use "Natural" method
                  [n_sol, n_est] := natural_superlog_size(d)
                  natural_superlog(base, z, n_sol, n_est)

             # this could be Re(base) <= 1 if desired
             # but the series may not converge
             |log(H(base))| <= 1:
                  # use "Regular" method
                  n := regular_superlog_size(d)
                  regular_superlog(base, z, n)
  
             otherwise: (undefined)
end

What are all these functions? like ceiling_superlog()?
These will be explained in a future post.

Third, a basic framework for superroot()
Code:
superroot(z, h) :=

  if h is integer

  then: case h == 0: (undefined)
             h == 1: z
             h == 2: ssqrt(z)
             h == infinity: z^(1/z)
             otherwise:
                # use a power series method
                Lagrange_inverse_series(
                  tetrate(base, h) == z, base )

  else:
    # use a numerical method
    Newton_method( tetrate(base, h) - z == 0, base )

Unfortunately I have noticed very little progress in the calculation of superroots.

Andrew Robbins
Reply
#2
Before we move on, lets overview the natural method for a bit.

There are actually 4-5 "natural" methods that I can think of for obtaining the coefficients of superlog(). Each has certain constants (indicated by 0 or the like), and variables (indicated by #) which can be set to almost anything.
  • (b=E, z_0=0, n_sol=#, n_est=0) -- Peter Walker's original method
  • (b=#, z_0=0, n_sol=#, n_est=0) -- An improved method I described in my first paper
  • (b=#, z_0=#, n_sol=#, n_est=0) -- An improved method I described in this post
  • (b=#, z_0=0, n_sol=#, n_est=#) -- Jay Fox's accelerated method described here, here, and here.
  • (b=#, z_0=#, n_sol=#, n_est=#) -- Unknown

Since the global behavior of the coefficients of the super-logarithm are still not very well known, the only expansion point which can be estimated so far seems to be . Maybe Jay can correct me on this. So in order to take advantage of the accelerated superlog, we must use the expansion about zero. All input values must be appropriately exponentiated or logarithmized, before applying the series.

In general, the natural method uses the a matrix that represents the Abel function, the solution to the matrix equation is a vector which represents the coefficients of the Abel function (superlog in this case). To summarize what Jay Fox found, instead of solving for all of the coefficients, we estimate all coefficients from 50..100, and we only solve for coefficients 1..50 for example. By estimating the later coefficients, the accuracy of the early coefficients is improved. I will call this the natural_superlog() function. Since it depends on an approximation number (called n), this will be another argument to the function.

I am very impressed with Jay Fox's error bound. This would allow us to ask questions like "What is e^^pi to 14 significant figures?" since we would know what is the matrix size that would allow this precision. According to his error bound, 14 digits would require a matrix using my method. With Jay's method however, there seems to be a variable that does not exist with my method, the size of the estimated/accelerated coefficients. Jay's method can reduce the matrix size down to 2000 or 200 depending on the ratio between solved/natural coefficients and estimated/accelerated coefficients. It seems that the more estimated coefficients there are, the greater the precision. Tuning this number to increase the speed of the calculations will be most valuable in calculating the fixed points of tetration and pentation as Ivars and GFR have mentioned.

Here is an example matrix for the accelerated natural super-log approach:

Code:
[ # # # # # # # ] [ X_1 ]   [ 1 ]
[ # # # # # # # ] [ X_2 ]   [ 0 ]
[ # # # # # # # ] [ X_3 ]   [ 0 ]
[ # # # # # # # ] [ X_4 ] = [ 0 ]
[ 0 0 0 0 1 0 0 ] [ X_5 ]   [ # ]
[ 0 0 0 0 0 1 0 ] [ X_6 ]   [ # ]
[ 0 0 0 0 0 0 1 ] [ X_7 ]   [ # ]

what Jay's method is doing is replacing the bottom part of the Abel matrix with the identity matrix, and replacing the 1-vector with the estimated coefficients. While this speeds up the solving process, it also seems to increase the accuracy and precision of the final solution. This matrix shows that there are two variables that dictate the matrix (aside from the base b): the number of solved coefficients and the number of estimated coefficients . In the example above these numbers would be and giving a matrix.

Interestingly, Jay has talked about solving the matrix given above, then simply putting more estimated coefficients in the series. This doesn't actually increase the accuracy, but decreases it, because instead of obeying the Abel functional equation exactly, it obeys it to 0.000001 accuracy, for example. This method of first solving, THEN using estimated coefficients is equivalent to the matrix:

Code:
[ # # # # 0 0 0 ] [ X_1 ]   [ 1 ]
[ # # # # 0 0 0 ] [ X_2 ]   [ 0 ]
[ # # # # 0 0 0 ] [ X_3 ]   [ 0 ]
[ # # # # 0 0 0 ] [ X_4 ] = [ 0 ]
[ 0 0 0 0 1 0 0 ] [ X_5 ]   [ # ]
[ 0 0 0 0 0 1 0 ] [ X_6 ]   [ # ]
[ 0 0 0 0 0 0 1 ] [ X_7 ]   [ # ]

which does not increase accuracy of the solution, but makes it worse! This is because the very nice estimates that we have for the outer coefficients are not being used in the solution at all. So when you are constructing an algorithm to do this, make sure that you implement the first matrix, NOT the second matrix. The first matrix will make it better, the second matrix will make it worse. And if you didn't want to use estimated coefficients at all, that would make the matrix look like this:

Code:
[ # # # # ] [ X_1 ]   [ 1 ]
[ # # # # ] [ X_2 ]   [ 0 ]
[ # # # # ] [ X_3 ]   [ 0 ]
[ # # # # ] [ X_4 ] = [ 0 ]

as in my original paper. The reason why it takes so long for these coefficients to converge using my method is because of the implicit assumption that when in fact or something.

This method can be used to increase the accuracy and precision of the coefficients of the natural super-logarithm, but how do they affect the accuracy and precision of the final answer? I found a paper that solves this problem exactly, called The Taylor Polynomial Error Formula, which allows you to ask what precision can be found when only using the first n terms of a Taylor series, i.e. a Taylor polynomial.

This formula can be used to derive the series-precision from the coefficient-precision that Jay has found.

Applying this formula to superlog will hopefully be the subject of my next post.

Andrew Robbins
Reply
#3
Some remarks from my side:

1. This is a great undertaking
2. Dont mix different methods into one tetration function. Each method should have its own dedicated function name(or one function which you can provide with the method as parameter). This is essential if you want to compare different methods.
3. The computation of the regular iteration is quite well-developed in the literature. There is the direct numeric method and there is the method via powerseries coefficients. There is a numeric method to directly compute the regular Abel function, as well as the regular iteration (without explicitly using function inversion). You can also directly compute the powerseries coefficients of the exponential of the Abel function.
4. Though I trust Jay's intuition, I think there was no proof that the accellerated method indeed yields the same limit as the unaccelareted/your original method.
5. Determining error-bounds would be a great thing, I hope you can make it. However is this possible if we have no proof of the convergence of the coefficients yet? Is it on the line that *if* they converge, then they would have these errorbounds? Perhaps we need extra research about the error-bounds of the other methods. Generally I think properly dealing with error-bounds is a quite messy thing. You need some kind of function how the precision of the output depends on the precision of the input, so that you can compute the overall error of a computation. Often then it is necessary to drastically increase the precision of intermediate steps.
6. I am not quite sure whether to develop it in pseudo code really makes sense. I mean if this is several pages long, you dont know whether there are errors in it, and who wants to read that?
The dilemma here is that we are quite differently equipped with mathematical software. I would suggest to implement it in a real language (and hence have the ability of debugging) or/and describe all methods in a scientific enough way that it should be easy to implement. For the implementation the language of choice would be Sage, as nobody is restrained from its use.
Reply
#4
I have learned alot about Sage by building this library. For example, in my implementation of parabolic iteration, I used "SymbolicVariable" to make a new variable, where I have since learned that Sage makes a simpler function available for this, namely "var" so you can do "t = var('t')" instead. I have also learned that every Sage module has a sub-module called "all" which allows you to import the most commonly used functions from that package at the same time. So where I used to do:
Code:
from sage.calculus.calculus import SymbolicVariable
now I use
Code:
from sage.all import *
which is a quick-and-dirty way to make all of Sage available when making a new Python module. I have also learned the difference between ".sage" and ".py" files, in that ".sage" and ".sobj" files are loaded with "load('filename.sage')" whereas ".py" files are loaded with "import filename". So instead of naming my files with ".sage" I will be naming them with ".py" so they can be used with the Python module system. So far I have finalized two files so far, and while I'm still working on the other files, I think there is enough to start using.

The 2 files I have available so far are "hyper/special/matrix.py" and "hyper/polynomial.py" which I have attached. I am very close to finishing matrix exp, log, and power, but I need help with the implementation. I see that the methods discussed in this thread talk about 2 methods of calculating matrix log and I'm still wondering which is better.

To use these, or try them out, just put the "hyper" directory in the directory you are going to be running sage. So if you run sage in "/home/me" then you should see "/home/me/hyper/" before you run Sage. You can type the following at the Sage prompt:
Code:
sage: from hyper.all import *
sage: p_poly(x)
x^5*C5 + x^4*C4 + x^3*C3 + x^2*C2 + x
You can test to see if its working by just making a parabolic function. The coefficients are pre-declared as symbolic variables in Sage, so you don't have to make your own. You can now test The Bell matrix with this:
Code:
sage: Bell_matrix(p_poly(x), x)

[  1   0              0              0     0  0]
[  0   1              0              0     0  0]
[  0  C2              1              0     0  0]
[  0  C3           2*C2              1     0  0]
[  0  C4    2*C3 + C2^2           3*C2     1  0]
[  0  C5 2*C4 + 2*C2*C3  3*C3 + 3*C2^2  4*C2  1]
sage: Bell_matrix(exp(x) - 1, x, n_row=6)

[     1      0      0      0      0      0      0]
[     0      1      0      0      0      0      0]
[     0    1/2      1      0      0      0      0]
[     0    1/6      1      1      0      0      0]
[     0   1/24   7/12    3/2      1      0      0]
[     0  1/120    1/4    5/4      2      1      0]
[     0  1/720 31/360    3/4   13/6    5/2      1]
sage: Carleman_matrix(x + 1, x, n_row=8)

[ 1  0  0  0  0  0  0  0  0]
[ 1  1  0  0  0  0  0  0  0]
[ 1  2  1  0  0  0  0  0  0]
[ 1  3  3  1  0  0  0  0  0]
[ 1  4  6  4  1  0  0  0  0]
[ 1  5 10 10  5  1  0  0  0]
[ 1  6 15 20 15  6  1  0  0]
[ 1  7 21 35 35 21  7  1  0]
[ 1  8 28 56 70 56 28  8  1]

Have Fun! Also, if you want to learn which functions are in this package, since it is obviously going to be growing (or learn about Sage!) is to do something like this:
Code:
sage: import hyper.all
sage: dir(hyper.all)
this will give you a listing of all of the functions in that module.

@Gottfried
I think it would also be nice to write two more functions for your use, namely V = Vandermonde_vector() and Vandermonde_matrix(), but I haven't done this yet...

@Henryk
I agree with all of your points. I will do it in Sage.

Andrew Robbins


Attached Files
.gz   hyperlib-20080423.tar.gz (Size: 2.2 KB / Downloads: 282)
Reply
#5
Hi Andrew -

concerning the matrix-exponential and matrixlogarithm...

Are you aware of the article "19 dubious ways to compute the matrix-exponential" of Moler/van Loan? They work out, that the "naive" exponential-series is just one of the poorest selections... (if it is difficult to find online, I cut put a copy on my webspace - I won't attach it here as usual, since it is about 3 MB)
So I didn't expect my own "naive" implementation in Pari/Gp as not really worth to publish, for instance. It should be made explicitely clear for which range of parameters/types of matrices the exponential-series is appropriate. The same is even more true for the matrix-logarithm.

On the other hand, for cases, where the "naive" methods are applicable, I can provide some enhancements - optimization for computation and acceleration of sums using Euler-summation-coefficients. I could send my Pari/Gp-routines for this (or display them here). However, the Euler-summation is not yet "dynamic", which means it doesn't itself evaluate the given matrix-parameter to find the appropriate order but the user must give a guess (and find out the best order by trial&error).

If this would be helpful, please let me know.

Gottfried
Gottfried Helms, Kassel
Reply
#6
So I wrote a few more files. To start, I implemented a simple version of the right hyper operations, which can be found in "hyper/hyper.py", and to test these you can try this:
Code:
sage: from hyper.all import *
sage: hyper(3)(3, 2)
9
sage: hyper(4)(3, 2)
27
sage: hyper(5)(3, 2)
7625597484987
sage: plot(lambda x: hyper(4)(e.n(20), x), -1.9, 2.1)
You get the idea. Smile

I have also implemented a wrapper around Jay Fox's infinitely iterated exponential (what he calls 'FindFixed', I renamed it to KnoebelH), and copied the Wikipedia implementation of LambertW, which does not work very well below -1/e, so I hope to improve upon this using the research of Corless et.al. so, for now, we can use it for simple things like:
Code:
sage: LambertW(1)
0.567143290409784
sage: KnoebelH(e)
0.318131505204764 + 1.33723570143069*I
sage: KnoebelH(e^(pi/2))
1.49809129296906e-13 + 1.00000000000024*I
Well, almost...


Attached Files
.gz   hyper-20080429.tar.gz (Size: 6.5 KB / Downloads: 268)
Reply
#7
Try it out with a very interesting example:
Code:
sage: from hyper.all import *
sage: t = var('t')
sage: p = parabolic_flow_matrix(exp(x) - 1, t, x, 0, 10)
sage: def c(k): return map(lambda x: x[-k], p[k:])
....:
(Just press enter again when it looks like '....').
The flow matrix 'p' will give the coefficients of x and t in .

The following will give the first diagonal of the coefficients of x and t in the continuous iteration of , which is known to be as Jay Fox noticed here and we can ensure that it fits this pattern by dividing by it.
Code:
sage: [c(1)[k] for k in xrange(9)]
[1, 1/2, 1/4, 1/8, 1/16, 1/32, 1/64, 1/128, 1/256]
sage: [c(1)[k]*2**k for k in xrange(9)]
[1, 1, 1, 1, 1, 1, 1, 1, 1]

Now for the fun part! Since we defined the helper function 'c', we can do the same thing for the second diagonal, which up until now has been a mystery.
Code:
sage: [c(2)[k] for k in xrange(1,8)]
[-1/12, -5/48, -13/144, -77/1152, -29/640, -223/7680, -481/26880]
sage: [c(2)[k]*2**k*(-6)/harmonic_number(k, 2) for k in xrange(1,8)]
[1, 1, 1, 1, 1, 1, 1]
as you can see, the second diagonal of the power series of iterated is related to harmonic numbers, see the implementation notes for more details. This means that the power series of iterated is of the form:

Note that this is not the normal generalized harmonic numbers, but the ones used by Conway & Guy, which are slightly different. These harmonic numbers are defined by:


The reason why I took the time to include such an advanced algorithm for computing harmonic numbers, is because I believe they hold the key to a new closed form, identity, or tautology for iterated-dxp.

Andrew Robbins
Reply
#8
andydude Wrote:
Code:
sage: LambertW(1)
0.567143290409784
sage: KnoebelH(e)
0.318131505204764 + 1.33723570143069*I
sage: KnoebelH(e^(pi/2))
1.49809129296906e-13 + 1.00000000000024*I
Well, almost...

Andrew - I came across the similar problem with the wikipedia Lambert-W and thus still use an own routine which replaces this for some parameters.
It is a sort of "lasso"-procedure, the mathematics is *not* proven.
It makes use of the observation, that if you start with a complex value in the near of the fixpoint, and iterate, then the sequence of values form a very nice spiral, with the fixpoint near the center.

So if you have one round, you just use the center as new initial value, which provides then a much more narrow new lasso.

The spiral converges, if it is an attracting fixpoint, but because you switch to the center after the first round, you possibly save lots of tidy rounds on the spiral.
The special effect is: you can even catch repelling fixpoints, if the spiral diverges not too strong - just because the catching of the "center" of one round can over-compensate the divergence.
For instance I found
Code:
findfixpt(exp(Pi/2))
%180 = -1.0451585870966735600158078879579 E-41 + 1.0000000000000000000000000000000*I

findfixpt(exp(Pi/2),1e-80)
%181 = -1.0250080565445841470968576986582 E-81 + 1.0000000000000000000000000000000*I
in a blink with adjustable precision.

If you want, I post the pari/gp routine here for re-implementation/improving. It's quick&dirty, fast only on some range of parameters (seems to cover the slow range of the wikipedia-lambert, for instance) and I didn't investigate the math below the surface.
Gottfried Helms, Kassel
Reply
#9
Gottfried Wrote:If you want, I post the pari/gp routine here for re-implementation/improving. It's quick&dirty, fast only on some range of parameters (seems to cover the slow range of the wikipedia-lambert, for instance) and I didn't investigate the math below the surface.

It seems like you have all the best code! Could you post your pari/gp code for the matrix-exp and LambertW? I also remember that you posted some other pari/gp code somewhere, but I can't find the post...

Andrew Robbins
Reply
#10
andydude Wrote:
Gottfried Wrote:If you want, I post the pari/gp routine here for re-implementation/improving. It's quick&dirty, fast only on some range of parameters (seems to cover the slow range of the wikipedia-lambert, for instance) and I didn't investigate the math below the surface.

It seems like you have all the best code! Could you post your pari/gp code for the matrix-exp and LambertW? I also remember that you posted some other pari/gp code somewhere, but I can't find the post...

Andrew Robbins

Here the complete findfixpt-code, where the lambert-w-wikipedia method is selected for some range of b
[update] I also should mention, that I work with default precision of 200 digits in pari/gp. So if you implement, for instance, the matrix-exponential using a lower precision the results may not be as expected...

To check the appropriateness of the parameter-range-separation you should examine the three alternative fixpoint-searcher again. I don't have my earlier ad-hoc-checks at hand and just left this for future investigation, while I'm using mostly "simple" bases b.
[/update]

Code:
\\ finds a fixpoint for (real) base-parameter b.
\\ Depending on its value the best of three methods are chosen
{ findfixpt(b, prec=1e-40, maxiter=1000) = local(r);
   if(b>=exp(-exp(1))&b<=exp(exp(-1)),return(LH(b))); \\ 1) use "LH" for some b
   if(b<4,return(findfp(b,maxiter,prec)));            \\ 2) use "findfp" for some other b
   return(findfpLog(b,prec,maxiter));                 \\ 3) use "findfplog" for all other b
}

\\ 3) --------------------------------------------------------------
{ findfpLog(b, prec=1e-40, maxiter=1000) = local(lb=1/log(b), x2, x1=b+0.1);
  \\ using hint by Henryk
    for(k=1,maxiter,
           x2=log(x1)*lb;
           if(abs(x2-x1)<prec,return(x2));
           x1=x2
        );
    return(0); }


\\ 2) LH(x) - uses wikipedia Lambert-W -----------------------------------

{ LW(x, prec=1E-80, maxiters=200) = local(w, we, w1e);
  \\ implementation by translation of wikipedia-python source
   w=0;
  for(i=1,maxiters
    ,we=w*exp(w);
     w1e=(w+1)*exp(w);
     if(prec>abs((x-we)/w1e),return(w),w=w-(we-x)/(w1e-(w+2)*(we-x)/(2*w+2)))
     );
  print("W doesn't converge fast enough for ",x,"( we=",we);
  return(0); }

LH(t11, prec=1E-80, maxiters=200) = exp(-LW(-log(t11),prec,maxiters))



\\ 1) findfp - "Lasso-method"  uses fdiff1 ------------------------------------------------
  \\ fdiff1 - uses vmean

vmean(v, n) = return(sum(k=1,n,v[k])/n);

\\ iterates a function tmp=(base^tmp + a*tmp)/(a+1) which circles around the fixpoint
\\ "a" is an empirical parameter to enhance convergence, here set to optimal(?) value 2
\\ performs one "round" of the spiral, returns its center
\\ workarrays vectors "list" and "diff" are provided by calling function, so they
\\ need not be created at each function-call
{ fdiff1(base=2.0, start=1+I, a=2.0, list, diff) = local(tmp1, tmp, findex);
   tmp1=start;findex=1;
   for(k=1,length(list), \\ compute iterations of function(see below) for some length
       list[k]=tmp1;
       \\ is the current iteration near the first value at k=1 (one round closed)?
       diff[k]=abs(tmp1-list[1]);
       if(k>1
         ,if(findex==1
             ,if(diff[k]<diff[k-1],findex=-1)
             ,if(diff[k]>diff[k-1],return([k,tmp1,vmean(list,k)])))); \\ return center of round
       tmp=tmp1;
       tmp1=(base^tmp+a*tmp)/(a+1); \\ this seems to spiral around the fixpoint
      );
   return([0,0,0]) ; }



{ findfp(s1, maxiter=60, prec=1e-44) = local(diff, list, res, val, oldval);
    diff=vectorv(100);list=vectorv(100);
    val=0.824678546142+1.56743212385*I; \\ an empirically useful init-point. May be enhanced in another version
    oldval=val;
    for(k=1,maxiter,
       res=fdiff1(s1,val,,list,diff);
       val=res[3];  \\ use the center-point as new guess
       if(abs(s1^val-val)<prec,return(val));
       oldval=val;
       );
     return(val); }

\\ ----------------------------------------------------------

The optimized matrixexponential simply uses a cache, but no improvement over the "naive" exponential-series is made. (to learn about problems see Moler/Van Loan: "19 dubious ways to compute the matrix-exponential")

Code:
\\ performs matrix-exponential using matrix-cache
\\ optimized "naive" implementation of exp-series for matrices
\\ uses a matrix-cache for keeping an array of "hmatrices"-count powers of the argument.
\\ then
\\ exp(X) = I + X/1! + X^2/2! + X^3/3! + ...
\\ is replaced by, for instance cache-size "hmatrices= 3"
\\ cache = [I,X,X^2], X3 = X^3
\\ exp(X) =   I   *(cache[1]/0! + cache[2]/1! + cache[3]/2!)
\\          + X3  *(cache[1]/3! + cache[2]/4! + cache[3]/5!)
\\          + X3^2*(cache[1]/6! + cache[2]/7! + cache[3]/8!)
\\          + ...
\\ where also the powers of X3 are sequentially computed "on the fly"
\\ Cache should have size at least >8 to benefit from optimization, while many terms are requested
\\ uses MAsum(X) = sum of absolute values of entries in X as terminating-criterion
\\ for nilpotent matrices provide the parameters terms=0, prec=0


{ MExpOpt(m, terms=0, prec=1e-32, hmatrices=8) = local(ff1, ff2, ff3, hm, h, s, k, rs, cs);
  rs=rows(m);cs=cols(m);if(rs<>cs,print("MExpOpt: no square matrix");return(m));
  hm=vector(hmatrices);
  ff1=matid(rs);
   for(h=0,hmatrices-1,
        hm[1+h]=ff1;ff1=ff1*m);
  if(terms==0,terms=rs);
  ff2=matid(rs);
  s=sum(h=0,hmatrices-1,hm[1+h]/h!);
  forstep(k=hmatrices,terms,hmatrices,
           ff2=ff2*ff1;
           ff3=ff2*sum(h=0,hmatrices-1,hm[1+h]/(k+h)!);
           s=s+ff3;
           if(prec<>0,if(MASum(ff3)<prec,print(k," terms used");break()));
          );
  return(s); }
Gottfried Helms, Kassel
Reply




Users browsing this thread: 1 Guest(s)