R devel - Dec 2016 - Different results for cos,sin,tan and cospi,sinpi,tanpi

Hi,

i try sin, cos, and tan.
> sapply(c(cos,sin,tan),function(x,y)x(y),1.23e45*pi)
[1] 0.5444181 0.8388140 1.5407532

However, *pi results the following
> sapply(c(cospi,sinpi,tanpi),function(x,y)x(y),1.23e45)
[1] 1 0 0

Please try whether the following becomes all right.

diff -ruN R-3.3.2.orig/src/nmath/cospi.c R-3.3.2/src/nmath/cospi.c
--- R-3.3.2.orig/src/nmath/cospi.c      2016-09-15 07:15:31.000000000 +0900
+++ R-3.3.2/src/nmath/cospi.c   2016-12-01 13:54:38.863357149 +0900
@@ -35,7 +35,11 @@
 #endif
     if(!R_FINITE(x)) ML_ERR_return_NAN;

-    x = fmod(fabs(x), 2.);// cos() symmetric; cos(pi(x + 2k)) =cos(pi x) for
all integer k
+    x = fabs(x);
+    if ( x > 9007199254740991 ) /* 2^53-1 */
+        return cos(M_PI * x);
+
+    x = fmod(x, 2.);// cos() symmetric; cos(pi(x + 2k)) == cos(pi x)
for all integer k
     if(fmod(x, 1.) == 0.5) return 0.;
     if( x == 1.)       return -1.;
     if( x == 0.)       return  1.;
@@ -57,6 +61,9 @@
 #endif
     if(!R_FINITE(x)) ML_ERR_return_NAN;

+    if (( x >  9007199254740991 )||  /*  2^53-1 */
+        ( x < -9007199254740991 )  ) /* -2^53-1 */
+        return sin(M_PI * x);
     x = fmod(x, 2.); // sin(pi(x + 2k)) == sin(pi x)  for all integer k
     // map (-2,2) --> (-1,1] :
     if(x <= -1) x += 2.; else if (x > 1.) x -= 2.;
@@ -81,6 +88,10 @@
 #endif
     if(!R_FINITE(x)) ML_ERR_return_NAN;

+    if (( x >  9007199254740991 )||  /*  2^53-1 */
+        ( x < -9007199254740991 )  ) /* -2^53-1 */
+        return tan(M_PI * x);
+
     x = fmod(x, 1.); // tan(pi(x + k)) == tan(pi x)  for all integer k
     // map (-1,1) --> (-1/2, 1/2] :
     if(x <= -0.5) x++; else if(x > 0.5) x--;

-- 
Best Regards,
--
Eiji NAKAMA <nakama (a) ki.rim.or.jp>
"\u4e2d\u9593\u6804\u6cbb"  <nakama (a) ki.rim.or.jp>

>>>>> Ei-ji Nakama <nakama at ki.rim.or.jp>
>>>>>     on Thu, 1 Dec 2016 14:39:55 +0900 writes:
    > Hi,
    > i try sin, cos, and tan.

    >> sapply(c(cos,sin,tan),function(x,y)x(y),1.23e45*pi)
    > [1] 0.5444181 0.8388140 1.5407532

    > However, *pi results the following

    >> sapply(c(cospi,sinpi,tanpi),function(x,y)x(y),1.23e45)
    > [1] 1 0 0

    > Please try whether the following becomes all right.

    [..............................]

Yes, it does  -- the fix will be in all future versions of R.

Thank you very much Ei-ji Nakama, for this valuable contribution
to make R better!

Martin Maechler,
ETH Zurich


    > -- 
    > Best Regards,
    > --
    > Eiji NAKAMA <nakama (a) ki.rim.or.jp>
    > "\u4e2d\u9593\u6804\u6cbb"  <nakama (a) ki.rim.or.jp>

>>>>> Martin Maechler <maechler at stat.math.ethz.ch>
>>>>>     on Thu, 1 Dec 2016 09:36:10 +0100 writes:
>>>>> Ei-ji Nakama <nakama at ki.rim.or.jp>
>>>>>     on Thu, 1 Dec 2016 14:39:55 +0900 writes:
    >> Hi,
    >> i try sin, cos, and tan.

    >>> sapply(c(cos,sin,tan),function(x,y)x(y),1.23e45*pi)
    >> [1] 0.5444181 0.8388140 1.5407532

    >> However, *pi results the following

    >>> sapply(c(cospi,sinpi,tanpi),function(x,y)x(y),1.23e45)
    >> [1] 1 0 0

    >> Please try whether the following becomes all right.

    > [..............................]

    > Yes, it does  -- the fix will be in all future versions of R.

oops.... not so quickly, Martin!

Of course, the results then coincide,  by sheer implementation.

*BUT* it is not at all clear which of the two results is better;
e.g., if you replace '1.23' by '1' in the above examples, the
result of the unchnaged  *pi() functions is 100% accurate,
whereas

 R> sapply(c(cos,sin,tan), function(Fn) Fn(1e45*pi))
 [1] -0.8847035 -0.4661541  0.5269043

is "garbage".  After all,  1e45 is an even integer and so, the
(2pi)-periodic functions should give the same as for 0  which
*is*  (1, 0, 0).

For such very large arguments, the results of all of sin() ,
cos() and tan()  are in some sense "random garbage" by
necessity:
Such large numbers have zero information about the resolution modulo
[0, 2pi)  or (-pi, pi]  and hence any (non-trivial) periodic
function with such a "small" period can only return "random
noise".


    > Thank you very much Ei-ji Nakama, for this valuable contribution
    > to make R better!

That is still true!  It raises the issue to all of us and will
improve the documentation at least!

At the moment, I'm not sure where we should go.
Of course, I could start experiments using my own 'Rmpfr'
package where I can (with increasing computational effort!) get
correct values (for increasingly larger arguments) but at the
moment, I don't see how this would help.

Martin

    > Martin Maechler,
    > ETH Zurich


    >> -- 
    >> Best Regards,
    >> --
    >> Eiji NAKAMA <nakama (a) ki.rim.or.jp>
    >> "\u4e2d\u9593\u6804\u6cbb"  <nakama (a)
ki.rim.or.jp>

    > ______________________________________________
    > R-devel at r-project.org mailing list
    > https://stat.ethz.ch/mailman/listinfo/r-devel