• 1 Vote(s) - 5 Average
• 1
• 2
• 3
• 4
• 5
 fast accurate Kneser sexp algorithm sheldonison Long Time Fellow Posts: 683 Threads: 24 Joined: Oct 2008 08/20/2010, 01:05 AM (08/19/2010, 02:35 AM)sheldonison Wrote: Later, I will probably post results showing how the Taylor and Fourier series terms converge, and how many terms it takes to get appropriate precision in the theta(z) fourier series so as to generate sufficiently accurate samples for the sexp(z) series......Here's an example, settings which I like as the next step up from the term term double precision results, which take less than a minute. Type "init(exp(1));loop(14);" and this next example runs for 7-8 minutes on my laptop, and reports that resulting sexp(z) Taylor series is accurate to 108 binary bits, or just a bit over 32 decimal digits. That precision would be within a radius of abs(z)<=0.5. \p 67; ctrm = 75; cvec = 150; rtrm = 1000; rvec = 2000; sfuncprec = 1E-32; Here, gp's precision is set to 67 digits. The goal is 32 decimal digits of precision. The sfuncprec (superfunction precision) is set to 1E-32 digits, to give superf/isuperf results accurate to 32 digits. That is, the isuperf/superf iteration count was set to 235 iterations. If you take the $\log^{[235]}(0.5)$, the result will be within 10^-32 of the fixed point of L. That is roughly what it takes to get 32 digits of accuracy from the superfunction. We compare exp(L+1E-32) with the fixed point approximation, L+L*1E-32, and we notice that the exponent function has an x^2 term around 1E-64. Thus, the maximum precision we're going to be left with is just a bit over 32 decimal digits. 67-32 leaves 35 digits of precision for gp to work with, or 3 spare bits for rounding errors, which surprisingly seems sufficient! We're asking for a Taylor series with 75 terms. By the Nyquist sampling criteria, we need at least 150 sample points. So cvec (the number of sample terms used to generate the sexp(z) Taylor series) =150. There will be a small sampling error term, due to higher frequencies in the superfunction being sampled. But, for base e sexp(z), the higher frequency terms roll off pretty fast, with each term being roughly half the magnitude of the term before. We want results accurate to 32 decimal digits, so we need 150 sample points, accurate to 32 decimal digits. The sample points will be in the upper half of the complex plane, centered around z=0. The points closest to the real axis are $\exp(\pi i /300)$ and $\exp(299\pi i /300)$. This is at $\Im(z)=0.0105$ rtrm = 1000 is driven by the requirement that the Riemann approximation, which is $\text{superf}(z+\theta(z))$ be accurate to 32 decimal digits, at $\Im(z)=0.0105$. superf is accurate to 32 decimal digits, so the error term is determined by theta(z), which has a radius of convergence at Im(z)=0. The 1000th term in theta(z) has a magnitude of 2.26E-05. The theta(z)=a_n*exp(z*2pi*i). With the known value for Im(z) closest to the real axis, we get magnitude of the a_n term is 6.02*10E-34, which is smaller than the goal of 10E-32. So, using rtrm~=1000 is ok. Sampling for the $\theta(z)$ function requires 2000 samples by the Nyquist criterion, so we need rvec (sample points for the RiemannCircle/theta(z) function) >= 2000. The high frequency roll off is much more gradual, so we need a little bit of extra precision here. We also need to show that our 75 term sexp(z) Taylor series supports 32 decimal digits of precision at a radius with abs(z)<=0.5. This is for the unit length we use to generate the RiemannCircle/Theta function from. The unit length will be from -0.5+idelta to +0.5+idelta. We want all terms over the unit length accurate to 32 decimal digits. The 75th term is -7.15 E-25. 7.15E-25*(0.5^75)=1.9E-47, which is much much less then 1E-32. The sexp(z) Taylor series has 32 decimal digits of precision out to a radius of about 0.78, which is greater then 1/2. So, all of the parameters are sufficient to get around 32 bits of precision. Observationally, when iterating the "loop(14)", we get about 108 binary bits equal in the sexp(z) algorithm versus the riemaprx(z) algorithm. This is about 33 decimal digits. Graphing the error gainst a golden copy with higher precision shows similar results. - Sheldon « Next Oldest | Next Newest »

 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 RE: fast accurate Knesser sexp algorithm - by bo198214 - 08/08/2010, 03:54 PM RE: fast accurate Knesser sexp algorithm - by sheldonison - 08/08/2010, 06:46 PM RE: fast accurate Kneser sexp algorithm - by nuninho1980 - 08/15/2010, 12:09 AM 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 RE: fast accurate Kneser sexp algorithm - by sheldonison - 10/14/2010, 10:00 PM RE: fast accurate Kneser sexp algorithm - by nuninho1980 - 10/15/2010, 04:03 PM RE: fast accurate Kneser sexp algorithm - by sheldonison - 10/15/2010, 09:20 PM kneser.gp modified for large bases - by sheldonison - 10/19/2010, 03:33 AM small update to allow more flexibility, faster - by sheldonison - 10/30/2010, 09:47 PM update to support B

 Possibly Related Threads... Thread Author Replies Views Last Post Kneser-iteration on n-periodic-points (base say \sqrt(2)) Gottfried 11 4,092 05/05/2021, 04:53 AM Last Post: Gottfried "Kneser"/Riemann mapping method code for *complex* bases mike3 2 10,035 08/15/2011, 03:14 PM Last Post: Gottfried Attempt to make own implementation of "Kneser" algorithm: trouble mike3 9 24,209 06/16/2011, 11:48 AM Last Post: mike3 Numerical algorithm for Fourier continuum sum tetration theory mike3 12 30,435 09/18/2010, 04:12 AM Last Post: mike3 regular sexp: curve near h=-2 (h=-2 + eps*I) Gottfried 2 9,082 03/10/2010, 07:52 AM Last Post: Gottfried Attempting to compute the kslog numerically (i.e., Kneser's construction) jaydfox 11 28,758 10/26/2009, 05:56 PM Last Post: bo198214 regular sexp:different fixpoints Gottfried 6 18,311 08/11/2009, 06:47 PM Last Post: jaydfox sexp(strip) is winding around the fixed points Kouznetsov 8 20,608 06/29/2009, 10:05 AM Last Post: bo198214 sexp and slog at a microcalculator Kouznetsov 0 4,812 01/08/2009, 08:51 AM Last Post: Kouznetsov

Users browsing this thread: 1 Guest(s)