Thread Rating:
  • 0 Vote(s) - 0 Average
  • 1
  • 2
  • 3
  • 4
  • 5
fast accurate Kneser sexp algorithm
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: 191)
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

Messages In This Thread
The pari-GP code - by sheldonison - 08/07/2010, 09:17 PM
updated code - by sheldonison - 08/19/2010, 02:35 AM
RE: updated code - by nuninho1980 - 08/19/2010, 12:08 PM
RE: updated code - by sheldonison - 08/20/2010, 01:05 AM
small update to allow more flexibility, faster - by sheldonison - 10/30/2010, 09:47 PM
update to support B<eta - by sheldonison - 11/15/2010, 02:53 PM
RE: update to support B<eta - by nuninho1980 - 11/15/2010, 03:26 PM
another new version - by 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,438 08/15/2011, 03:14 PM
Last Post: Gottfried
  Attempt to make own implementation of "Kneser" algorithm: trouble mike3 9 7,585 06/16/2011, 11:48 AM
Last Post: mike3
  Numerical algorithm for Fourier continuum sum tetration theory mike3 12 10,061 09/18/2010, 04:12 AM
Last Post: mike3
  regular sexp: curve near h=-2 (h=-2 + eps*I) Gottfried 2 3,268 03/10/2010, 07:52 AM
Last Post: Gottfried
  Attempting to compute the kslog numerically (i.e., Kneser's construction) jaydfox 11 9,236 10/26/2009, 05:56 PM
Last Post: bo198214
  regular sexp:different fixpoints Gottfried 6 5,913 08/11/2009, 06:47 PM
Last Post: jaydfox
  sexp(strip) is winding around the fixed points Kouznetsov 8 6,620 06/29/2009, 10:05 AM
Last Post: bo198214
  sexp and slog at a microcalculator Kouznetsov 0 1,866 01/08/2009, 08:51 AM
Last Post: Kouznetsov

Users browsing this thread: 1 Guest(s)