R version 3.0.3 RC (2014-02-27 r65086) -- "Warm Puppy"
Copyright (C) 2014 The R Foundation for Statistical Computing
Platform: x86_64-unknown-linux-gnu (64-bit)

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.

> ####=== Numerical / Arithmetic Tests
> ####--- ALL tests here should return  TRUE !
> ###
> ### '##P': These lines don't give TRUE but relevant ``Print output''
> 
> ### --> d-p-q-r-tests.R  for distribution things
> 
> .proctime00 <- proc.time()
> opt.conformance <- 0
> Meps <- .Machine $ double.eps
> 
> ## this uses random inputs, so set the seed
> set.seed(1)
> 
> options(rErr.eps = 1e-30)
> rErr <- function(approx, true, eps = .Options$rErr.eps)
+ {
+     if(is.null(eps)) { eps <- 1e-30; options(rErr.eps = eps) }
+     ifelse(Mod(true) >= eps,
+ 	   1 - approx / true, # relative error
+ 	   true - approx)     # absolute error (e.g. when true=0)
+ }
> 
> abs(1- .Machine$double.xmin * 10^(-.Machine$double.min.exp*log10(2)))/Meps < 1e3
[1] TRUE
> ##P (1- .Machine$double.xmin * 10^(-.Machine$double.min.exp*log10(2)))/Meps
> if(opt.conformance)#fails at least on SGI/IRIX 6.5
+ abs(1- .Machine$double.xmax * 10^(-.Machine$double.max.exp*log10(2)))/Meps < 1e3
> 
> ## More IEEE  Infinity/NaN checks
> i1 <- pi / 0
> i1 == (i2 <- 1:1 / 0:0)
[1] TRUE
> is.infinite( i1) & is.infinite( i2) &   i1 > 12   &   i2 > 12
[1] TRUE
> is.infinite(-i1) & is.infinite(-i2) & (-i1) < -12 & (-i2) < -12
[1] TRUE
> 
> is.nan(n1 <- 0 / 0)
[1] TRUE
> is.nan( - n1)
[1] TRUE
> 
> i1 ==  i1 + i1
[1] TRUE
> i1 ==  i1 * i1
[1] TRUE
> is.nan(i1 - i1)
[1] TRUE
> is.nan(i1 / i1)
[1] TRUE
> 
> 1/0 == Inf & 0 ^ -1 == Inf
[1] TRUE
> 1/Inf == 0 & Inf ^ -1 == 0
[1] TRUE
> 
> iNA <- as.integer(NA)
> !is.na(Inf) & !is.nan(Inf) &   is.infinite(Inf) & !is.finite(Inf)
[1] TRUE
> !is.na(-Inf)& !is.nan(-Inf)&   is.infinite(-Inf)& !is.finite(-Inf)
[1] TRUE
>  is.na(NA)  & !is.nan(NA)  &  !is.infinite(NA)  & !is.finite(NA)
[1] TRUE
>  is.na(NaN) &  is.nan(NaN) &  !is.infinite(NaN) & !is.finite(NaN)
[1] TRUE
>  is.na(iNA) & !is.nan(iNA) &  !is.infinite(iNA) & !is.finite(iNA)
[1] TRUE
> 
> ## These are "double"s:
> all(!is.nan(c(1.,NA)))
[1] TRUE
> all(c(FALSE,TRUE,FALSE) == is.nan(c   (1.,NaN,NA)))
[1] TRUE
> ## lists are no longer allowed
> ## all(c(FALSE,TRUE,FALSE) == is.nan(list(1.,NaN,NA)))
> 
> 
> ##  log() and "pow()" -- POSIX is not specific enough..
> log(0) == -Inf
[1] TRUE
> is.nan(log(-1))# TRUE and warning
[1] TRUE
Warning message:
In log(-1) : NaNs produced
> 
> rp <- c(1:2,Inf); rn <- rev(- rp)
> r <- c(rn, 0, rp, NA, NaN)
> all(r^0 == 1)
[1] TRUE
> ir <- suppressWarnings(as.integer(r))
> all(ir^0  == 1)
[1] TRUE
> all(ir^0L == 1)# not in R <= 2.15.0
[1] TRUE
> all( 1^r  == 1)# not in R 0.64
[1] TRUE
> all(1L^r  == 1)
[1] TRUE
> all(1L^ir == 1)# not in R <= 2.15.0
[1] TRUE
> all((rn ^ -3) == -((-rn) ^ -3))
[1] TRUE
> #
> all(c(1.1,2,Inf) ^ Inf == Inf)
[1] TRUE
> all(c(1.1,2,Inf) ^ -Inf == 0)
[1] TRUE
> .9 ^ Inf == 0
[1] TRUE
> .9 ^ -Inf == Inf
[1] TRUE
> ## Wasn't ok in 0.64:
> all(is.nan(rn ^ .5))# in some C's : (-Inf) ^ .5 gives Inf, instead of NaN
[1] TRUE
> 
> 
> ## Real Trig.:
> cos(0) == 1
[1] TRUE
> sin(3*pi/2) == cos(pi)
[1] TRUE
> x <- rnorm(99)
> all( sin(-x) == - sin(x))
[1] TRUE
> all( cos(-x) == cos(x))
[1] TRUE
> 
> x <- 1:99/100
> all(abs(1 - x / asin(sin(x))) <= 2*Meps)# "== 2*" for HP-UX
[1] TRUE
> all(abs(1 - x / atan(tan(x))) <  2*Meps)
[1] TRUE
> 
> ## Sun has asin(.) = acos(.) = 0 for these:
> ##	is.nan(acos(1.1)) && is.nan(asin(-2)) [!]
> 
> ## gamma()
> abs(gamma(1/2)^2 - pi) < 4* Meps
[1] TRUE
> r <- rlnorm(5000) # NB random, and next has failed for some seed
> all(abs(rErr(gamma(r+1), r*gamma(r))) < 500 * Meps)
[1] TRUE
> ## more accurate for integers n <= 50 since R 1.8.0	Sol8: perfect
> n <-   20; all(		 gamma(1:n) == cumprod(c(1,1:(n-1))))# Lnx: up too n=28
[1] TRUE
> n <-   50; all(abs(rErr( gamma(1:n), cumprod(c(1,1:(n-1))))) < 20*Meps)#Lnx: f=2
[1] TRUE
> n <-  120; all(abs(rErr( gamma(1:n), cumprod(c(1,1:(n-1))))) < 1000*Meps)
[1] TRUE
> n <- 10000;all(abs(rErr(lgamma(1:n),cumsum(log(c(1,1:(n-1)))))) < 100*Meps)
[1] TRUE
> 
> n <-   10; all(		 gamma(1:n) == cumprod(c(1,1:(n-1))))
[1] TRUE
> n <-   20; all(abs(rErr( gamma(1:n), cumprod(c(1,1:(n-1))))) < 100*Meps)
[1] TRUE
> n <-  120; all(abs(rErr( gamma(1:n), cumprod(c(1,1:(n-1))))) < 1000*Meps)
[1] TRUE
> n <- 10000;all(abs(rErr(lgamma(1:n),cumsum(log(c(1,1:(n-1)))))) < 100*Meps)
[1] TRUE
> 
> all(is.nan(gamma(0:-47))) # + warn.
[1] TRUE
Warning message:
In gamma(0:-47) : NaNs produced
> 
> ## choose() {and lchoose}:
> n51 <- c(196793068630200, 229591913401900, 247959266474052)
> abs(c(n51, rev(n51))- choose(51, 23:28)) <= 2
[1] TRUE TRUE TRUE TRUE TRUE TRUE
> all(choose(0:4,2) == c(0,0,1,3,6))
[1] TRUE
> ## 3 to 8 units off and two NaN's in 1.8.1
> 
> ## psi[gamma](x) and derivatives:
> ## psi == digamma:
> gEuler <- 0.577215664901532860606512# = Euler's gamma
> abs(digamma(1) + gEuler) <   32*Meps # i386 Lx: = 2.5*Meps
[1] TRUE
> all.equal(digamma(1) - digamma(1/2), log(4), tolerance = 32*Meps)# Linux: < 1*Meps!
[1] TRUE
> n <- 1:12
> all.equal(digamma(n),
+          - gEuler + c(0, cumsum(1/n)[-length(n)]),tolerance = 32*Meps)#i386 Lx: 1.3 Meps
[1] TRUE
> all.equal(digamma(n + 1/2),
+           - gEuler - log(4) + 2*cumsum(1/(2*n-1)),tolerance = 32*Meps)#i386 Lx: 1.8 Meps
[1] TRUE
> ## higher psigamma:
> all.equal(psigamma(1, deriv=c(1,3,5)),
+           pi^(2*(1:3)) * c(1/6, 1/15, 8/63), tolerance = 32*Meps)
[1] TRUE
> x <- c(-100,-3:2, -99.9, -7.7, seq(-3,3, length=61), 5.1, 77)
> ## Intel icc showed a < 1ulp difference in the second.
> stopifnot(all.equal( digamma(x), psigamma(x,0), tolerance = 2*Meps),
+           all.equal(trigamma(x), psigamma(x,1), tolerance = 2*Meps))# TRUE (+ NaN warnings)
Warning messages:
1: In digamma(x) : NaNs produced
2: In psigamma(x, 0) : NaNs produced
> ## very large x:
> x <- 1e30 ^ (1:10)
> a.relE <- function(appr, true) abs(1 - appr/true)
> stopifnot(a.relE(digamma(x),   log(x)) < 1e-13,
+           a.relE(trigamma(x),     1/x) < 1e-13)
> x <- sqrt(x[2:6]); stopifnot(a.relE(psigamma(x,2), - 1/x^2) < 1e-13)
> x <- 10^(10*(2:6));stopifnot(a.relE(psigamma(x,5), +24/x^5) < 1e-13)
> 
> ## fft():
> ok <- TRUE
> ##test EXTENSIVELY:	for(N in 1:100) {
>     cat(".")
.>     for(n in c(1:30, 1000:1050)) {
+ 	x <- rnorm(n)
+ 	er <- Mod(rErr(fft(fft(x), inverse = TRUE)/n, x*(1+0i)))
+ 	n.ok <- all(er < 1e-8) & quantile(er, 0.95, names=FALSE) < 10000*Meps
+ 	if(!n.ok) cat("\nn=",n,": quantile(rErr, c(.95,1)) =",
+ 		      formatC(quantile(er, prob= c(.95,1))),"\n")
+ 	ok <- ok & n.ok
+     }
>     cat("\n")

> ##test EXTENSIVELY:	}
> ok
[1] TRUE
> 
> ## var():
> for(n in 2:10)
+     print(all.equal(n*(n-1)*var(diag(n)),
+ 		    matrix(c(rep(c(n-1,rep(-1,n)),n-1), n-1), nr=n, nc=n),
+ 		    tolerance = 20*Meps)) # use tolerance = 0 to see rel.error
[1] TRUE
[1] TRUE
[1] TRUE
[1] TRUE
[1] TRUE
[1] TRUE
[1] TRUE
[1] TRUE
[1] TRUE
> 
> ## pmin() & pmax() -- "attributes" !
> v1 <- c(a=2)
> m1 <- cbind(  2:4,3)
> m2 <- cbind(a=2:4,2)
> 
> all( pmax(v1, 1:3) == pmax(1:3, v1) & pmax(1:3, v1) == c(2,2,3))
[1] TRUE
> all( pmin(v1, 1:3) == pmin(1:3, v1) & pmin(1:3, v1) == c(1,2,2))
[1] TRUE
> 
> oo <- options(warn = -1)# These four lines each would give 3-4 warnings :
>  all( pmax(m1, 1:7) == pmax(1:7, m1) & pmax(1:7, m1) == c(2:4,4:7))
[1] TRUE
>  all( pmin(m1, 1:7) == pmin(1:7, m1) & pmin(1:7, m1) == c(1:3,3,3,3,2))
[1] TRUE
>  all( pmax(m2, 1:7) == pmax(1:7, m2) & pmax(1:7, m2) == pmax(1:7, m1))
[1] TRUE
>  all( pmin(m2, 1:7) == pmin(1:7, m2) & pmin(1:7, m2) == c(1:3,2,2,2,2))
[1] TRUE
> options(oo)
> 
> ## pretty()
> stopifnot(pretty(1:15)	    == seq(0,16, by=2),
+ 	  pretty(1:15, h=2) == seq(0,15, by=5),
+ 	  pretty(1)	    == 0:1,
+ 	  pretty(pi)	    == c(2,4),
+ 	  pretty(pi, n=6)   == 2:4,
+ 	  pretty(pi, n=10)  == 2:5,
+ 	  pretty(pi, shr=.1)== c(3, 3.5))
> 
> ## gave infinite loop [R 0.64; Solaris], seealso PR#390 :
> all(pretty((1-1e-5)*c(1,1+3*Meps), 7) == seq(0,1,len=3))
[1] TRUE
> 
> n <- 1000
> x12 <- matrix(NA, 2,n); x12[,1] <- c(2.8,3) # Bug PR#673
> for(j in 1:2) x12[j, -1] <- round(rnorm(n-1), dig = rpois(n-1, lam=3.5) - 2)
> for(i in 1:n) {
+     lp <- length(p <- pretty(x <- sort(x12[,i])))
+     stopifnot(p[1] <= x[1] & x[2] <= p[lp],
+               all(x==0) || all.equal(p, rev(-pretty(-x)), tolerance = 10*Meps))
+ }
> 
> ## PR#741:
> pi != (pi0 <- pi + 2*.Machine$double.eps)
[1] TRUE
> is.na(match(c(1,pi,pi0), pi)[3])
[1] TRUE
> 
> ## PR#749:
> all(is.na(c(NA && TRUE, TRUE  && NA, NA && NA,
+             NA || FALSE,FALSE || NA, NA || NA)))
[1] TRUE
> 
> all((c(NA || TRUE,  TRUE || NA,
+     !c(NA && FALSE,FALSE && NA))))
[1] TRUE
> 
> 
> ## not sure what the point of this is: it gives mean(numeric(0)), that is NaN
> (z <- mean(rep(NA_real_, 2), trim = .1, na.rm = TRUE))
[1] NaN
> is.na(z)
[1] TRUE
> 
> ## Last Line:
> cat('Time elapsed: ', proc.time() - .proctime00,'\n')
Time elapsed:  0.237 0.003 0.24 0 0 
>