[Rd] 1/tan(-0) != 1/tan(0)
Martin Maechler
maechler at stat.math.ethz.ch
Wed Jun 1 10:21:06 CEST 2005
Testing the code that Morten Welinder suggested for improving
extreme tail behavior of qcauchy(),
I found what you can read in the subject.
namely that the tan() + floating-point implementation on all
four different versions of Redhat linux, I have access to on
i686 and amd64 architectures,
> 1/tan(c(-0,0))
gives
-Inf Inf
and of course, that can well be considered a feature, since
after all, the tan() function does jump from -Inf to +Inf at 0.
I was still surprised that this even happens on the R level,
and I wonder if this distinction of "-0" and "0" shouldn't be
mentioned in some place(s) of the R documentation.
For the real problem, the R source (in C), It's simple
to work around the fact that
qcauchy(0, log=TRUE)
for Morten's code proposal gives -Inf instead of +Inf.
Martin
>>>>> "MM" == Martin Maechler <maechler at stat.math.ethz.ch>
>>>>> on Wed, 1 Jun 2005 08:57:18 +0200 (CEST) writes:
>>>>> "Morten" == Morten Welinder <mwelinder at gmail.com>
>>>>> on Fri, 27 May 2005 20:24:36 +0200 (CEST) writes:
.............
Morten> Now that pcauchy has been fixed, it is becoming
Morten> clear that qcauchy suffers from the same problems.
Morten>
Morten> qcauchy(pcauchy(1e100,0,1,FALSE,TRUE),0,1,FALSE,TRUE)
Morten> should yield 1e100 back, but I get 1.633178e+16.
Morten> The code below does much better. Notes:
Morten> 1. p need not be finite. -Inf is ok in the log_p
Morten> case and R_Q_P01_check already checks things.
MM> yes
Morten> 2. No need to disallow scale=0 and infinite
Morten> location.
MM> yes
Morten> 3. The code below uses isnan and finite directly.
Morten> It needs to be adapted to the R way of doing that.
MM> I've done this, and started testing the new code; a version will
MM> be put into the next version of R.
MM> Thank you for the suggestions.
>>> double
>>> qcauchy (double p, double location, double scale, int lower_tail, int log_p)
>>> {
>>> if (isnan(p) || isnan(location) || isnan(scale))
>>> return p + location + scale;
>>> R_Q_P01_check(p);
>>> if (scale < 0 || !finite(scale)) ML_ERR_return_NAN;
>>> if (log_p) {
>>> if (p > -1)
>>> lower_tail = !lower_tail, p = -expm1 (p);
>>> else
>>> p = exp (p);
>>> }
>>> if (lower_tail) scale = -scale;
>>> return location + scale / tan(M_PI * p);
>>> }
More information about the R-devel
mailing list