Post Reply 
Thread Rating:
  • 0 Votes - 0 Average
  • 1
  • 2
  • 3
  • 4
  • 5
fast accurate Kneser sexp algorithm
10/30/2010, 09:47 PM (This post was last modified: 11/21/2011 09:38 PM by sheldonison.)
Post: #13
small update to allow more flexibility, faster
This is a small update, which improves performance, and allows more flexibility. Its a little faster (25 seconds for 105-110 binary bit precision, as opposed to 1 minute). I added an slog function, as well as a gentaylorseries function that will generate the sexp taylor series about any point in the complex plane. I cleaned up the code, and made it a little simpler, in that the "n" parameter for loop(n) is now optional, and the default base for init; is base e. You can still "loop(1,2...)" at a time if desired.

.gp (Size: 24.09 KB / Downloads: 187)
For the most recent code version: go to the Nov 21st, 2011 thread.
   base          2.71828182845904523536029
   fixed point   0.318131505204764135312654 +
.... 25 seconds later, after 13 loop iterations
gp > sexp(0.5)
gp > slog(I)
-1.3937809364969141717386716853218 + 1.0171347771806901474061430995874*I
gp > gentaylor(3*I,1)
complex sexp Taylor series centered at 3*I
gp > for (s=1,4,print(tseries[s])); /* prints first 4 series terms */
0.37090658903228120608281774428121 + 1.3368216707889122043183467452060*I
0.018300482687991735872774054492565 + 0.069611076949751462108351480615434*I
-0.042221079601603961951042313227903 + 0.024296334049072988793474267457353*I
-0.015851643810851135752181652414535 - 0.014789535958785875826049470725178*I
There is a much larger performance increase for "\p 134", which now takes about 5 1/2 minutes, as opposed to the previous version which took 32 minutes. "\p 134 gives sexp results accurate to 64 decimal digits. The new algorithm cleverly converges the sexp(z) with the riemaprx(z) at imag(z)=0.12*i, thereby avoiding all of the Riemann mapping singularity problems requiring larger and larger numbers of Riemann approximation Fourier series terms. That's where the speed up comes from. I had to fix a couple of pari-gp precision loss problems to make this new algorithm reliable for all bases, init(1.47) - init(100000). The earlier algorithm is still available via "slowmode=1". One reason for the new code is to allow a secondary program, to generate pentation, which I will be posting separately.

Another code example, showing \p 134, and the taylor series function.
\p 134
.... wait 5 1/2 minutes, after 26 iterations ....
gp > sexp(0.5)

type "morestuff" on the command line to see new functions
new functions:
slog(z,est)      inverse sexp(z), using optional est as initial estimate
gentaylor(w,r)   generates tseries sexp taylor series, centered at w, radius r
taprx(z)         evaulates tseries taylor series from gentaylor at z
This implementation shows that excellent convergence can be had at 0.12*I. Here the Cauchy unit circle which is sampled to generate the sexp(z) taylor series is divided into a region where imag(z) at the unit circle is >= 0.12*I. For that region, the Cauchy unit circle is sampled from the Riemann approximation generated from the 1-cyclic transformation of the super function. But, for points where imag(z)<0.12*I, the sample is taken from the exp or logarithm of previous version of sexp(z) itself. Then these samples are stitched together, to generate the next sexp(z) series. Then the next Riemann approximation Fourier series is generated at 0.12*I. The new algorithm required identifying and fixing the pari-gp precision problems, which results in much more predictable stable behavior, which convinces me that the convergence will continue with more and more iterations.
- Sheldon
Find all posts by this user
Quote this message in a reply
Post Reply 

Messages In This Thread
The pari-GP code - sheldonison - 08/07/2010, 09:17 PM
updated code - sheldonison - 08/19/2010, 02:35 AM
RE: updated code - nuninho1980 - 08/19/2010, 12:08 PM
RE: updated code - sheldonison - 08/20/2010, 01:05 AM
small update to allow more flexibility, faster - sheldonison - 10/30/2010 09:47 PM
update to support B<eta - sheldonison - 11/15/2010, 02:53 PM
RE: update to support B<eta - nuninho1980 - 11/15/2010, 03:26 PM
another new version - sheldonison - 11/17/2010, 06:52 PM

Possibly Related Threads...
Thread: Author Replies: Views: Last Post
  "Kneser"/Riemann mapping method code for *complex* bases mike3 2 3,373 08/15/2011 03:14 PM
Last Post: Gottfried
  Attempt to make own implementation of "Kneser" algorithm: trouble mike3 9 7,431 06/16/2011 11:48 AM
Last Post: mike3
  Numerical algorithm for Fourier continuum sum tetration theory mike3 12 9,818 09/18/2010 04:12 AM
Last Post: mike3
  regular sexp: curve near h=-2 (h=-2 + eps*I) Gottfried 2 3,182 03/10/2010 07:52 AM
Last Post: Gottfried
  Attempting to compute the kslog numerically (i.e., Kneser's construction) jaydfox 11 9,002 10/26/2009 05:56 PM
Last Post: bo198214
  regular sexp:different fixpoints Gottfried 6 5,779 08/11/2009 06:47 PM
Last Post: jaydfox
  sexp(strip) is winding around the fixed points Kouznetsov 8 6,485 06/29/2009 10:05 AM
Last Post: bo198214
  sexp and slog at a microcalculator Kouznetsov 0 1,831 01/08/2009 08:51 AM
Last Post: Kouznetsov

User(s) browsing this thread: 1 Guest(s)