From: Bill Gosper <billgosper@gmail.com> Wow, this Remezes <http://www.tweedledum.com/rwg/recipoddremez.png> to mu=2*10^-36 with 14 parameters! --rwg

--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. 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