From nakama at ki.rim.or.jp  Mon Dec  5 14:04:02 2016
From: nakama at ki.rim.or.jp (Ei-ji Nakama)
Date: Mon, 5 Dec 2016 22:04:02 +0900
Subject: [Rd] Different results for cos,sin,tan and cospi,sinpi,tanpi
In-Reply-To: <CAJqeyYY+boNQGDOrH2uSyAg9zFGFW5WFsGEgSDWJkptED5LXdw@mail.gmail.com>
References: <CAJqeyYZ2jozkWyvpBk9m=zeunh25ciFCfJdBx10ge3hbRUV_cg@mail.gmail.com>
	<22591.57594.122801.713377@stat.math.ethz.ch>
	<22591.59795.613756.145797@stat.math.ethz.ch>
	<af030acd-b3fb-be11-3449-52301a46b053@stats.ox.ac.uk>
	<CAJqeyYY+boNQGDOrH2uSyAg9zFGFW5WFsGEgSDWJkptED5LXdw@mail.gmail.com>
Message-ID: <CAJqeyYbbwABt-Ro1eAgsY6u4d4dPcJiFPwuocCPwGn9wyNMfkg@mail.gmail.com>

hello,

i read this pdf (http://www.open-std.org/jtc1/sc22/wg14/www/docs/n1950.pdf)
i think ... x to be integer only(more accurate value).

following test code and patch.
### F.10.1.12 The cospi functions
## cospi(+-0)                      : 1, 1
cospi(c(+0,-0))
## cospi(n + 1/2)                  : all 0
cospi((-4:4)+.5)
## cospi(+-Inf)                    : NaN
cospi(c(+Inf,-Inf))
## cospi(n)
cospi((-4:4))

### F.10.1.13 The sinpi functions
## sinpi(+-0)                      :+0,-0
sprintf("%a",sinpi(c(+0,-0)))
## sinpi(n + 1/2)                  : all 0
sinpi((-4:4))
## sinpi(+-Inf)                    : NaN
sinpi(c(+Inf,-Inf))
## sinpi(n)                        :1 -1  1 -1  1 -1  1 -1  1
sinpi((-4:4+.5))

### F.10.1.14 The tanpi functions
## tanpi(+-0)                       :+0,-0
sprintf("%a",tanpi(c(+0,-0)))
## tanpi(pos even and neg odd)      :+0
sprintf("%a",tanpi(c(-3,-1,2,4)))
## tanpi(pos odd and ned even)      :-0
sprintf("%a",tanpi(c(-4,-2,1,3)))
## tanpi(n+1/2) n = even            :+Inf
tanpi(c(1:3*2)+.5)
tanpi(c(1:3*2)*-1+.5)
## tanpi(n+1/2) n = odd             :-Inf
tanpi(c(1:3*2+1)+.5)
tanpi(c(1:3*2+1)*-1+.5)

## tanpi(+-Inf)                     :NaN NaN
tanpi(c(+Inf,-Inf))

## no integer
sinpi(1.23e23) # 0.4652223
cospi(1.23e23) # 0.8851939
tanpi(1.23e23) # 0.5255597



--- 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-05 21:29:20.764593514 +0900
@@ -21,13 +21,10 @@
    The __cospi etc variants are from macOS (and perhaps other
BSD-based systems).
 */

-#ifdef HAVE_COSPI
-#elif defined HAVE___COSPI
-double cospi(double x) {
-    return __cospi(x);
-}
+
+#if defined(__STDC_WANT_IEC_60559_FUNCS_EXT__) &&
__STDC_WANT_IEC_60559_FUNCS_EXT__ >= 201506L
+/* use standard cospi */
 #else
-// cos(pi * x)  -- exact when x = k/2  for all integer k
 double cospi(double x) {
 #ifdef IEEE_754
     /* NaNs propagated correctly */
@@ -35,7 +32,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.;
@@ -44,11 +45,8 @@
 }
 #endif

-#ifdef HAVE_SINPI
-#elif defined HAVE___SINPI
-double sinpi(double x) {
-    return __sinpi(x);
-}
+#if defined(__STDC_WANT_IEC_60559_FUNCS_EXT__) &&
__STDC_WANT_IEC_60559_FUNCS_EXT__ >= 201506L
+/* use standard cospi */
 #else
 // sin(pi * x)  -- exact when x = k/2  for all integer k
 double sinpi(double x) {
@@ -57,6 +55,12 @@
 #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);
+
+    if( x == 0 || x == -0 )
+        return(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.;
@@ -69,26 +73,50 @@
 #endif

 // tan(pi * x)  -- exact when x = k/2  for all integer k
-#if defined(HAVE_TANPI) || defined(HAVE___TANPI)
+#if defined(__STDC_WANT_IEC_60559_FUNCS_EXT__) &&
__STDC_WANT_IEC_60559_FUNCS_EXT__ >= 201506L
+/* use standard cospi */
 // for use in arithmetic.c, half-values documented to give NaN
-double Rtanpi(double x)
 #else
 double tanpi(double x)
-#endif
 {
+  int _sig=0;
+  int _even=0;
+  int _odd=0;
+  int _int=0;
 #ifdef IEEE_754
     if (ISNAN(x)) return x;
 #endif
     if(!R_FINITE(x)) ML_ERR_return_NAN;

-    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--;
-    return (x == 0.) ? 0. : ((x == 0.5) ? ML_NAN : tan(M_PI * x));
+    if (( x >  9007199254740991 )||  /*  2^53-1 */
+        ( x < -9007199254740991 )  ) /* -2^53-1 */
+        return tan(M_PI * x);
+
+    if( x == 0. || x == -0. )
+        return(x);
+    if(x>0) _sig = 1;
+    if(x<0) _sig =-1;
+
+    x = fmod(x, 2.);
+    if(( x == 0.0 )||( x == -0.0)){ _even = 1; _int=1;}
+    if(( x == 0.5 )||( x == -0.5))  _even = 1;
+    if(( x == 1.0 )||( x == -1.0)){ _odd = 1;  _int=1;}
+    if(( x == 1.5 )||( x == -1.5))  _odd = 1;
+    if(_int){
+      if( _sig ==  1 && _even ) return(  0.);
+      if( _sig == -1 && _odd  ) return(  0.);
+      if( _sig ==  1 && _odd  ) return( -0.);
+      if( _sig == -1 && _even ) return( -0.);
+    }
+    if(_even){
+        if ( x ==  0.5 ) return(R_PosInf);
+    if ( x == -0.5 ) return(R_NegInf);
+    }else if (_odd){
+    if ( x ==  1.5 ) return(R_NegInf);
+    if ( x == -1.5 ) return(R_PosInf);
+    }
+    // otherwise
+    return tan(M_PI * x);
 }

-#if !defined(HAVE_TANPI) && defined(HAVE___TANPI)
-double tanpi(double x) {
-    return __tanpi(x);
-}
 #endif




2016-12-01 18:45 GMT+09:00 Ei-ji Nakama <nakama at ki.rim.or.jp>:
> hi,
>
> my environment...
>> sessionInfo()
> R version 3.3.2 (2016-10-31)
> Platform: x86_64-pc-linux-gnu (64-bit)
> Running under: Debian GNU/Linux 8 (jessie)
>
> locale:
>  [1] LC_CTYPE=ja_JP.UTF-8       LC_NUMERIC=C
>  [3] LC_TIME=ja_JP.UTF-8        LC_COLLATE=ja_JP.UTF-8
>  [5] LC_MONETARY=ja_JP.UTF-8    LC_MESSAGES=ja_JP.UTF-8
>  [7] LC_PAPER=ja_JP.UTF-8       LC_NAME=C
>  [9] LC_ADDRESS=C               LC_TELEPHONE=C
> [11] LC_MEASUREMENT=ja_JP.UTF-8 LC_IDENTIFICATION=C
>
> attached base packages:
> [1] stats     graphics  grDevices utils     datasets  methods   base
>
> It's not a very good example...
>
> f0<-function(x,y)exp(complex(real=x,imag=y))
> f1<-function(x,y)complex(real=exp(1)^x*cos(y),imag=exp(1)^x*sin(y))
> f2<-function(x,y)complex(real=exp(1)^x*cospi(y/pi),imag=exp(1)^x*sinpi(y/pi))
>
> f0(700,1.23)
> f1(700,1.23)
> f2(700,1.23)
>
> f0(700,1.23e23)
> f1(700,1.23e23)
> f2(700,1.23e23)
>
> Garbage number is required.
>
> Thank you!
>
> 2016-12-01 18:31 GMT+09:00 Prof Brian Ripley <ripley at stats.ox.ac.uk>:
>> Please note that you need to report your platforms (as per the posting
>> guide), as the C function starts
>>
>> #ifdef HAVE_COSPI
>> #elif defined HAVE___COSPI
>> double cospi(double x) {
>>     return __cospi(x);
>> }
>>
>> And AFAICS the system versions on Solaris and OS X behave the same way as
>> R's substitute.
>>
>>
>>
>>
>> On 01/12/2016 09:12, Martin Maechler wrote:
>>>>>>>>
>>>>>>>> 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
>>
>>
>>
>> --
>> Brian D. Ripley,                  ripley at stats.ox.ac.uk
>> Emeritus Professor of Applied Statistics, University of Oxford
>
>
>
> --
> Best Regards,
> --
> Eiji NAKAMA <nakama (a) ki.rim.or.jp>
> "\u4e2d\u9593\u6804\u6cbb"  <nakama (a) ki.rim.or.jp>



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



[Rd] Different results for cos,sin,tan and cospi,sinpi,tanpi
Ei-ji Nakama