Thread Rating:
  • 0 Vote(s) - 0 Average
  • 1
  • 2
  • 3
  • 4
  • 5
fast accurate Kneser sexp algorithm
#1
Here is the most recent code, Nov 12th, 2013:
.gp   kneser.gp (Size: 35.84 KB / Downloads: 434)
updates Nov 12th 2013
o updated to be compatible with latest pari-gp
o removed eval, and kill
o enhanced prtpoly routine
o added lambertL routine
end updates Nov 12th, 2013

This is the first of several posts. update: changed the names of the routines to match current code.

I found a fast accurate way to implement the Knesser sexp, and implemented it in Pari-GP, which is a free mathematics program available on the web. The algorithm generates the Knesser sexp accurate to 15 decimal digits in about a minute. This includes a 400 term Riemann mapping approximation, and a 50 term Cauchy unit circle approximation to generate the Taylor series. I've tried the algorithm for base 2, base e, and base 10, and they are all well behaved. I will post the pari-GP code, as well as the Cauchy series it generated. The algorithm is for bases > eta.

The algorithm is to start with the linear approximation for the sexp, that has a smooth first derivative. For base e, this is the approximation that the sexp(z) from [-1 to 0]=z+1. This region is extended with the exponent/logarithm to extend the initial unit length estimate.

This is the initial trivial Taylor series. In my program, this Taylor series approximation for the sexp is called sexp.

Then iterate as follows:

1) generate the Riemann mapping of the unit length segment
This involves the Laurent series of the inverse superfunction of the unit length segment, mapped to the unit circle. The complex valued superfunction and its inverse are developed from the complex fixed point. Specifically, over a unit length, generate inverse_superfunction(sexp(z))-z, which is then mapped to the unit circle. Use the Cauchy Integral formula to generate the Taylor series. Strictly speaking, this function has a Laurent series, not a Taylor series, but we ignore the x^-n terms in the Laurent series, and only keep the x^+n terms, where n>=0. This is then trivially mapped back into a complex valued 1-cyclic fourier series, theta(z).

2) riemaprx(z). This 1-cyclic fourier series is used to generate a sexp approximation defined only in the imag(z)>0 portion of the plane. riemaprx(z)=superfunction(z+theta(z)). riemaprx(z) is poorly behaved at the real axis, with high frequency oscillations due to the slowly converging theta function. But the high frequency terms decay very rapidly as imag(z) increases in the complex plane, and if imag(z)>=1*I, the function behaves close to the complex superfunction.

3) sexp(z). Using riemaprx(z), generate the upper half of a unit half circle, conventionally centered on z=0, in the upper half of the complex plane. Treat the lower half of the unit circle as the complex conjugate of the upper half of the unit circle. Use the Cauchy integral to generate a Taylor series. This is a real valued Taylor series. The Taylor series is renormalized and recentered, so that sexp(0)=1, and so that sexp(+0.5)=exponent(sexp(-0.5).

Iterate steps 1) 2) and 3). The claim is that this infinite iteration converges to Knessder's sexp(z), which is real valued at the real axis, with sexp(0)=1, and sexp(z+1)=exp(sexp(z)). Seven iterations were required to get 15 digits of accuracy, which each iteration adding approximately 8 binary bits of accuracy. In the Pari-GP program, each iteration takes less than 10 seconds.

Implementing this algorithm is highly problematic and requires very large numbers of terms in the 1-cyclic theta function, due to the very slow convergence of the Riemann mapping. I found a trick to dramatically improve the convergence of the Riemann mapping, to handle the singularity, and to allow fast convergence.

Here a few more details, before I call it a night. In my implementation, I wanted a 50 term Taylor series for the sexp approximation. This requires 200 points around a unit circle. Since we only generate the upper half of the unit circle, that's 100 points. These points are evenly spaced, with all points having imag(z)>0. The first point has z=e^(I*pi/200), which has imag(z)=0.015. This means that in riemprx(z), the imag(z) will always be >=0.015. That in turns means that the high frequency oscillations decay. For Riemann mapping theta terms>~300 is not computationally significant for any of these 100 points on the unit circle.

This is the first of several posts. The next post will deal with the ideas required to improve the Riemann mapping convergence problems due to the singularity. I need to get some sleep, and it may take a few days to post all of my results. I want to clean up the code a little bit, since there's a lot of junk in the code that I needed to debug it and get it all working that I want to remove before posting, to avoid unnecessary confusion.

Teaser: This is the output of my program, to generate the Taylor series terms I generated for base e, base 2, and base 10.
- Sheldon

base 2.71828182845904523536028747135266249775724709369
Cauchy Taylor series approximation terms
0 0.999999999999999920606916968922030701567136645654
1 1.09176735125832094016082238370083804714
2 0.271483212901694628277009772276018058438
3 0.212453248176256341682340133816642561735
4 0.069540376139987455601815633290052837540
5 0.0442919520904733787858836715649205635588
6 0.0147367420963894641781675687873762988951
7 0.0086687818172253171998963830218805136384
8 0.00279647939838551204530926929318637318059
9 0.00161063129058431027448748226836099353371
10 0.00048992723148441335050668409792945192050
11 0.000288181071154071176125207353574308091920
12 0.000080094612538568642532817921403719623613
13 0.0000502911417938227014378026179786713377326
14 0.0000121837903449185625957879474875634578856
15 0.0000086655336673937799848419995967770476324
16 0.00000168778231933152083906008897805134285921
17 0.00000149325324858236409519824063327981138611
18 0.000000198760764215837789507617671297185393027
19 0.000000260867356010921328462837051102110834560
20 0.0000000147099541512649833075949006297264385980
21 0.0000000468344973323925542131416564470746086702
22 -0.00000000154924165843564586518006754343312538792
23 0.0000000087415107851444402981222259166563437222
24 -0.00000000112578730421439340057613203197171894664
25 0.00000000170795927016981912380839845219103887952
26 -0.000000000377858310553104789819392071420497066106
27 0.000000000349577878719151948745500164259927239404
28 -1.05377008157789121410377229095750742668 E-10
29 7.4590973139432053510241775086048122696 E-11
30 -2.71759784858227219934013489755749519135 E-11
31 1.64607673327765883270714032926236921945 E-11
32 -6.7418700690910586159590314585805404498 E-12
33 3.7253295931126574070795526274392527304 E-12
34 -1.63908465536567938669004255576025126318 E-12
35 8.5836440640244217926800645809445203856 E-13
36 -3.94371548069373997097518693401613112984 E-13
37 2.00252641630788814440567087736426238381 E-13
38 -9.4417589560533107254019594050715409934 E-14
39 4.71206677912917615338890731001808112440 E-14
40 -2.25611366289088681618159333858518060425 E-14
41 1.11546306823408543027263377496376191447 E-14
42 -5.3891796465832969287083377311762631343 E-15
43 2.65194707551149745587030410085037067028 E-15
44 -1.28753286810996291447909130936705727867 E-15
45 6.3232296660338767222943776801526973230 E-16
46 -3.07332373088892725649620921357557234127 E-16
47 1.50855996920297588144091655941600050227 E-16
48 -7.2896911298653544932621055870671863118 E-17
49 3.57053526604856938153249432482269325976 E-17

base 2
Cauchy Taylor series approximation terms
0 1.00000000000000007750805374695916476506324315386
1 0.88936495462097636471560015559202608188
2 0.0086765489653699778435866428210914580477
3 0.095238800075181821015052159076811795525
4 -0.0057523485401261058429044155126967586843
5 0.0129665820200371929900639963746185426267
6 -0.00219604962303099108685649788446966214674
7 0.00199674684791144956471972891306784395245
8 -0.00056335481487852293422791227955443116661
9 0.000348242328188164833713322066658348490588
10 -0.000128532441264722069102630177525101178612
11 0.000067081924420528746629841393628806788277
12 -0.0000282987528228002605436935419049427621654
13 0.0000138001319906297777277157382455652644180
14 -0.0000062019093983767935657716936885382308667
15 0.00000295556146480613946456248019931704026334
16 -0.00000136867922453685047021475603425081919190
17 0.00000064905707564831820958280933643834369294
18 -0.000000305166939330940495799879426965798767859
19 0.000000144948206147741006787309611484943294178
20 -0.000000068746643115668186072914364608901276865
21 0.0000000327674451745520233209052672738147657033
22 -0.0000000156310874697440742179774639803902393806
23 0.0000000074781283860060906963656804741326980492
24 -0.00000000358267681486784107606081902426567798118
25 0.00000000171986774278787862816161264617020453732
26 -0.00000000082681596398826102575453499011874838983
27 0.000000000398110465534940913092926440510368851600
28 -1.91942994009799320930467536291965520002 E-10
29 9.2663184086753526923394526399596161631 E-11
30 -4.47869896241156433452311745776051738114 E-11
31 2.16712006481128547956338987778919331607 E-11
32 -1.04969755427424430603006383413954325851 E-11
33 5.0894457415080947575828909582705030799 E-12
34 -2.46987949612974961256130795951078775734 E-12
35 1.19965332790083051126762500373890740311 E-12
36 -5.8316699693946087927290168849662635343 E-13
37 2.83700134396762172106948455139976794164 E-13
38 -1.38119296763538037439839695675052294198 E-13
39 6.7286246573391428225154736432311797155 E-14
40 -3.28040714229361317312282338553709457619 E-14
41 1.59994527980264052587204787980778074081 E-14
42 -7.8111883641343962658543820790964098310 E-15
43 3.8123345996113862398910808416038811110 E-15
44 -1.86468916410329958272266244696629139897 E-15
45 9.0928731462620461748731226835758948880 E-16
46 -4.4652219382963819619261489249575291192 E-16
47 2.16258309701983021090373383869474626354 E-16
48 -1.07562115928439172115561558722739104721 E-16
49 5.05103639405910509662889857211651054178 E-17

base 10
Cauchy Taylor series approximation terms
0 1.00000000003413666821707397610280580973
1 1.6842890280105318114702665361
2 1.4130640093559017784867641066
3 1.3295945895120850984491868131
4 1.1104423522130456306141162585
5 0.90791475960484107525431757972
6 0.70590812180827017012794979492
7 0.53596957102484588479532491192
8 0.39514882149633253685483274509
9 0.28568024972253825344818666302
10 0.20237131278470927324534438107
11 0.14110168269002467755684498200
12 0.096870196895171185635561049276
13 0.065642683799764598157183797364
14 0.043932559903871435514499480872
15 0.029083045516276459146860315440
16 0.019054947684195538710084470697
17 0.0123686846345218426147071658882
18 0.0079582306458762837173385697464
19 0.0050792074578440317507369752633
20 0.0032170388379653970425789785613
21 0.0020231559289094629314623265337
22 0.0012638006922446651516872615141
23 0.00078448455296479211241357653392
24 0.00048404347686816718054617660014
25 0.00029697714889618699433841203816
26 0.00018122446592630282867291500326
27 0.000110023183879937876785387661614
28 0.000066469809662926642710985803857
29 0.000039970210188298586919385209485
30 0.000023927948827218540851314911165
31 0.0000142631537294666506633793631158
32 0.0000084672146369634588559435171282
33 0.0000050067182976642942211046734889
34 0.0000029492933862829666416493809347
35 0.0000017310105129704227353019917697
36 0.0000010124091578563321764991987405
37 0.00000059013055820603257123183162307
38 0.00000034287053679750024168609816182
39 0.00000019859459586612750224033089360
40 0.000000114688324996137110089024467356
41 0.000000066049601906226671780697958378
42 0.000000037941008758273721620067345461
43 0.000000021745359672035340985615756584
44 0.0000000124403916641799613324962201452
45 0.0000000071074328972264520946865200883
46 0.0000000040598591783922772751609136103
47 0.0000000023194687970727546265839419100
48 0.0000000013297332794947389880115713152
49 0.00000000076366914316723959202509556617
Reply


Messages In This Thread
fast accurate Kneser sexp algorithm - by sheldonison - 08/07/2010, 06:53 AM
The pari-GP code - by sheldonison - 08/07/2010, 09:17 PM
updated kneser.gp code - by sheldonison - 08/19/2010, 02:35 AM
RE: updated kneser.gp code - by nuninho1980 - 08/19/2010, 12:08 PM
RE: updated kneser.gp code - by sheldonison - 08/20/2010, 01:05 AM
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 6,632 08/15/2011, 03:14 PM
Last Post: Gottfried
  Attempt to make own implementation of "Kneser" algorithm: trouble mike3 9 15,568 06/16/2011, 11:48 AM
Last Post: mike3
  Numerical algorithm for Fourier continuum sum tetration theory mike3 12 19,414 09/18/2010, 04:12 AM
Last Post: mike3
  regular sexp: curve near h=-2 (h=-2 + eps*I) Gottfried 2 5,889 03/10/2010, 07:52 AM
Last Post: Gottfried
  Attempting to compute the kslog numerically (i.e., Kneser's construction) jaydfox 11 18,122 10/26/2009, 05:56 PM
Last Post: bo198214
  regular sexp:different fixpoints Gottfried 6 11,627 08/11/2009, 06:47 PM
Last Post: jaydfox
  sexp(strip) is winding around the fixed points Kouznetsov 8 12,869 06/29/2009, 10:05 AM
Last Post: bo198214
  sexp and slog at a microcalculator Kouznetsov 0 3,136 01/08/2009, 08:51 AM
Last Post: Kouznetsov



Users browsing this thread: 2 Guest(s)