[R] Numerical stability of: 1/(1 - cos(x)) - 2/x^2

[Martin Maechler]

>>>>> Iris Simmons 
>>>>>     on Wed, 16 Aug 2023 02:57:48 -0400 writes:

    > You might also be able to rewrite
    > log(1 - cos(x))

    > as

    > log1p(-cos(x))

    > On Wed, Aug 16, 2023, 02:51 Iris Simmons <ikwsimmo using gmail.com> wrote:

    >> You could rewrite
    >> 
    >> 1 - cos(x)
    >> 
    >> as
    >> 
    >> 2 * sin(x/2)^2
    >> 
    >> and that might give you more precision?

Thank you, Iris,
yes, both suggestions are excellent *and* necessary.

Note that this is part  "numerical analysis 101"
and is called "cancellation":

Direct evaluation of  1 - cos(x)  for small 'x'  *cannot* ever
be numerically accurate and suffers from cancellation.

log(1+x) for small x  is slightly more subtle than pure
cancellation, but exactly the same reason we introduced  log1p()
{and expm1(); then even cospi(), sinpi() etc}
into R  even before it was widely adopted in C math libraries.

Martin


    >> On Wed, Aug 16, 2023, 01:50 Leonard Mada via R-help <r-help using r-project.org>
    >> wrote:
    >> 
    >>> Dear R-Users,
    >>> 
    >>> I tried to compute the following limit:
    >>> x = 1E-3;
    >>> (-log(1 - cos(x)) - 1/(cos(x)-1)) / 2 - 1/(x^2) + log(x)
    >>> # 0.4299226
    >>> log(2)/2 + 1/12
    >>> # 0.4299069
    >>> 
    >>> However, the result diverges as x decreases:
    >>> x = 1E-4
    >>> (-log(1 - cos(x)) - 1/(cos(x)-1)) / 2 - 1/(x^2) + log(x)
    >>> # 0.9543207
    >>> # correct: 0.4299069
    >>> 
    >>> I expected the precision to remain good with x = 1E-4 or x = 1E-5.
    >>> 
    >>> This part blows up - probably some significant loss of precision of
    >>> cos(x) when x -> 0:
    >>> 1/(cos(x) - 1) - 2/x^2
    >>> 
    >>> Maybe there is some room for improvement.
    >>> 
    >>> Sincerely,
    >>> 
    >>> Leonard
    >>> ==========
    >>> The limit was part of the integral:
    >>> up = pi/5;
    >>> integrate(function(x) 1 / sin(x)^3 - 1/x^3 - 1/(2*x), 0, up)
    >>> (log( (1 - cos(up)) / (1 + cos(up)) ) +
    >>> + 1/(cos(up) - 1) + 1/(cos(up) + 1) + 2*log(2) - 1/3) / 4 +
    >>> + (1/(2*up^2) - log(up)/2);
    >>> 
    >>> # see:
    >>> 
    >>> https://github.com/discoleo/R/blob/master/Math/Integrals.Trig.Fractions.Poly.R
    >>> 
    >>> ______________________________________________
    >>> R-help using r-project.org mailing list -- To UNSUBSCRIBE and more, see
    >>> https://stat.ethz.ch/mailman/listinfo/r-help
    >>> PLEASE do read the posting guide
    >>> http://www.R-project.org/posting-guide.html
    >>> and provide commented, minimal, self-contained, reproducible code.
    >>> 
    >> 

    > [[alternative HTML version deleted]]

    > ______________________________________________
    > R-help using r-project.org mailing list -- To UNSUBSCRIBE and more, see
    > https://stat.ethz.ch/mailman/listinfo/r-help
    > PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
    > and provide commented, minimal, self-contained, reproducible code.