[Rd] Different results for cos,sin,tan and cospi,sinpi,tanpi
Ei-ji Nakama
nakama at ki.rim.or.jp
Mon Dec 5 14:04:02 CET 2016
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>
More information about the R-devel
mailing list