From: Bill Gosper <billgosper@gmail.com <http://gosper.org/webmail/src/compose.php?send_to=billgosper%40gmail.com>>> Wow, this Remezes <http://www.tweedledum.com/rwg/recipoddremez.png> to> mu=2*10^-36 with 14 parameters!
Oop, more like 3*10^-36
--rwg
But noting |a[15]/a[1]| > |b[12]/1|, switch a[17] for b[12] <http://www.tweedledum.com/rwg/a17b10oddrecip.png> (which is a fairly notorious placebo anyway). μ = .6 * 10^-36 ! wds> --there is little need to Remez a fast-converging Chebyshev series since it already is Remezed to decent accuracy.
Oh, but I see, Gosper is not making a polynomial, but rather a rational approximation. The way he did it, though, he's Remezing something different -- really you want the approximation error of the rational function to be equiripple, while Gosper is making something else be equiripple. [That probably cost Gosper 50% more error, or something?] There are known convergent algorithms that do this right.
I think we just needed a different weight function. (Although this produced a lovely plot, I'm completely confused by NSolve producing several solutions, with the only desirable one claiming μ = 0. By the way, I also disunderstand your half integer offset.) Despite the anomalies, I'll stick my neck out and claim this will be hard to beat. With two surplus digits of precision, you might hope we could punt b[10], but its > a[17]. --rwg wds> Thinking about it, the difference-based method I proposed (odd sym) seems superior to the sum-based (even sym) one because later when solving a*b=p, a+b=s for a,b... with the sum you will need to do a bad subtraction which loses a lot of significant digits in cancellation sqrt(s*s/4 - p) is subtracting two nearly equal quantities while with difference a-b=d that issue does not arise: sqrt(d*d/4+p) Also to correct what I said before the solution of a*b = p, a-b = d is: 1/a = 1/ [ -d/2 +- sqrt(d*d/4+p) ] = [ d/2 +- sqrt(d*d/4+p) ] / p where the last expression had been erroneously reciprocated. Also, Gosper earlier had the idea, which is still a good one with this approach, that he could get symmetry for the Gamma function around x=1/2, OR get it around x=0, either way, and thus make a PAIR of symmetric approximants now each with half the interval length. That is, the interval is now 1/2 wide. The two functions you then need to approximate are d1(x) and d2(x) where |x|<=1/4... or perhaps you prefer to approximate d1(x)/x and d2(x)/x to keep the relative error small near x=0...: F(x) = 1/GAMMA(x+1/2). F(x) * F(-x) = cos(Pi*x)/Pi d1(x) = F(x) - F(-x). G(x) = 1/(x!) G(x) * G(-x) = sin(Pi*x) / (Pi*x) d2(x) = G(x) - G(-x) . d1(0) = d2(0) = 0. d1(1/4) = -d1(-1/4) = (1/2)*(2*Pi-sqrt(2)*GAMMA(3/4)^2)/(GAMMA(3/4)*Pi) = 0.54023327626805366671714040736252461351838415204125 d2(1/4) = -d2(-1/4) = (sqrt(8)*GAMMA(3/4)^2-Pi)/(GAMMA(3/4)*Pi) = 0.28721371222257427636269621260241438575441350426278 d1(1/2) = 1 d2(1/2) = Pi^(-1/2) d1'(0) = (2*gamma+4*ln(2))/Pi^(1/2) = 2.2155838077457420463420706011273861532011964034195 d2'(0) = 2*gamma = 1.1544313298030657212130241801648048620843186718798