Re: Minimax approximation for cosine

From: Richard Mathar (mathar_at_amer.strw.leidenuniv.nl)
Date: 11/18/04 

Date: Thu, 18 Nov 2004 21:01:19 +0000 (UTC)

In article <9215d7ac.0411141747.23763c56@posting.google.com>,
 glenlow@pixelglow.com (Glen Low) writes:
>> It might be better to use identities to reduce x to the interval [0,
>> Pi/2], and then divide that interval into two subintervals, say [0,
>> beta*Pi] and [beta*Pi, Pi/2].
>
>The Maclaurin expansion of cos(x) around x = 0 is a polynomial with
>even powers. Therefore it's highly likely that a minimax polynomial
>with even powers only will be a better fit than one with all powers.
>Such a polynomial will have the additional benefit of P(x) = P(-x)
>because of the even powers, so that I can just use the interval [0,
>Pi/2] to get a result that works correctly with [-Pi/2, Pi/2].
>
>After some thought on the bus to the day job, I came up with this.
>
>Usually in minimax approximations, you divide out any zeros of the
>function because otherwise the algorithm chokes e.g. if f(x) = 0 for x
>= k, then you can't do a minimax approximation on f(x), but you have
>to do it on f(x)/(x-k).
>
>The obvious function for cos(x) in [0, Pi/2] is cos(x)/(x-Pi/2) since
>cos(Pi/2)=0. But the resulting polynomial has all powers, and then
>multiply by (x-Pi/2) still gives me all powers back.
>
>So I thought I'd try: cos(x)/{(x-Pi/2)(x+Pi/2)} = cos(x)/(x^2 -
>Pi^2/4)... at x = Pi/2, this tends to 2/Pi. Substituing x -> x^.5, we
>have cos(x^.5)/(x - Pi^2/4). Assuming that P'(x) is the minimax
>polynomial for this, then P(x) = P'(x^2)(x^2 - Pi^2/4), so nicely P(x)
>has all even powers too.
>
>> Use a minimax polynomial for cos(x) for
>> x in [0, beta*Pi] and a minimax polynomial for sin(Pi/2 - x) for x in
>> [beta*Pi, Pi/2], where beta is chosen to balance the accuracy and the
>> execution time of the two approximations.
>
>I did the analysis for sin(x) and found a polynomial of degree 9 with
>odd powers only was sufficient to get enough accuracy for IEEE 754
>float (24 bits). I tried a couple of splits of the range [0, Pi/2] but
>none actually reduced the polynomial degree while maintaining
>accuracy. Of course if there is such a split, the hardware pipeline
>for fma would allow me to do both intervals almost simultaneously, so
>there'd be no loss in speed.
>
>I'm assuming that a similar situation would apply with cosine, we'll
>see...
>
>> How much accuracy do you need? I'm sure that there are newer books on
>> the topic, but I might be able to suggest values of beta and the
>> coefficients from an ancient (1967) book on approximation of
>> functions.
>
>The relative error in sin(x) or cos(x) must be < 2^-24 or about
>5.96e-8.
>
>Cheers,
>Glen Low, Pixelglow Software
>www.pixelglow.com

One can extend the example 5.6 in my script in
http://arXiv.org/abs/math.CA/0403344
for the Chebyshev polynomial for cos(pi/2*x)/(1-x^2) by two higher orders
than actually quoted in the example in the script.
This yields the expansion coefficients for the Chebyshev polynomials

b[0] = .890365196792210693146129746935
b[2] = -.107274469488517592983676091401
b[4] = .233370052015861433579008625011e-2
b[6] = -.264479401843810274578975201202e-4
b[8] = .184330753236829605876190693440e-6

or
.999999999071824518322659698297-.233700504645726537695303066267*x^2+.\
199685982106370967412399614777e-1*x^4-.893522758728821257757025461367e-3*x^6+.\
235943364143141895521524087603e-4*x^8

which should be good for a relative error of about 1.05e-9 in this
combined function. If one wants to go for the Chebyshev approximation
(some sort of minimum *absolute* error) for cos(pi*x/2)/(1-x^2):
these are also given in terms of Bessel functions there.

Relevant Pages

  • Re: Programmers unpaid overtime.
    ... >> involving a sum of powers in one or more variables multiplied ... The phrase "your accuracy and attention to detail" is lifted whole ... >> that polynomials are restricted to a highest power of two. ... That's why I felt free to reverse them from the ...
    (comp.programming)
  • Re: Programmers unpaid overtime.
    ... > involving a sum of powers in one or more variables multiplied ... > that polynomials are restricted to a highest power of two. ... Is your life so empty and lonely that negative attention ... I sure don't think I'd hang out with people who had no regard ...
    (comp.programming)
  • Re: ambiguities
    ... >2) polynomial can only have terms with powers that are positive ... Nobody thinks polynomials can only have positive powers: ... The constant term is both a term and a coefficient. ...
    (sci.math)
  • Re: Fractional part
    ... >with integer a, b, c, p, q and r have a same root. ... polynomials from each other to reduce the degree. ... This would mean that the original trinomials have ... made by someone else that the even powers of sqrtall have the ...
    (sci.math)
  • Re: Polynomial Approximation L-inf norm
    ... approximation for a similar problem. ... What is the numerical rank of A? ... interval with high order polynomials is a killer. ...
    (comp.soft-sys.matlab)