On a pi-related note: one algorithm that would be useful in Base is a remainder modulo-pi. At the moment we call out to a routine in openspecfun, but this requires allocating a 2-element array on each call (and also has some issues on certain compilers (see, for example https://github.com/JuliaLang/julia/issues/8799).
Having this would also make it possible to implement a faster "sincos" function (https://github.com/JuliaLang/julia/issues/10442). I've made a start on translating it here: https://github.com/JuliaLang/julia/compare/master...simonbyrne:sincos However it still requires an implementation of the Payne-Hanek range reduction for large (>1e6) arguments. For reference, the openspecfun code is available here: https://github.com/JuliaLang/openspecfun/blob/master/rem_pio2/k_rem_pio2.c Anyone interested? -Simon On Monday, 9 March 2015 17:47:23 UTC, Hans W Borchers wrote: > > Steven mentioned he would be interested in computing pi through a slowly > converging series. Well, I think the Leibniz series might be a good example: > > pi/4 = 1 - 1/3 + 1/5 - 1/7 + 1/9 -+ ... > > Let us define the first hundred terms of this series and see what their > sum is: > > julia> k = [1:100]; > julia> leibniz = (-1).^(k-1) ./ (2*k-1); > > julia> 4*sum(leibniz) > 3.1315929035585537 > > This is not even correct in the second digit after the decimal point. > > What can we do to accelerate the convergence of an alternating sum? A very > powerful method has been proposed by Henri Cohen, F. Rodriguez Villegas, > and Don Zagier (2000). The following function is my implementation of > 'Algorithm 1' from their paper: > > function sumalt(s_alt, n) > b = 2^(2*n-1) > c = b > s = 0.0 > for k in (n-1):-1:0 > t = s_alt[k+1] > s = s + c*t > b = b * (2*k+1) * (k+1) / (2 * (n-k) * (n+k)) > c = c + b > end > s = s / c > return(s) > end > > Let us apply this to the first 20(!) terms of the Leibniz series: > > julia> 4*sumalt(leibniz, 20) > 3.1415926535897936 > > The result is pi with an accuracy of 2*eps() = 4.44e-16 or 15 digits. > > Using Big Floats with 256 bit, we could easily reach 75 correct digits. Of > course, `sumalt` is general and can be applied in other contexts as well, > for example to compute the `[z]eta` function for complex values very > accurately. > >
