The R Project SVN R

R : Copyright 1999, The R Development Core Team
Version 0.90.0 Under development (unstable) (October 22, 1999)

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 a list.

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

> ####=== Numerical / Arithmetic Tests
> ####--- ALL tests here should return  TRUE !
> ###
> ### '##P': This lines give not 'T' but relevant ``Print output''
> 
> ### --> d-p-q-r-tests.R  for distribution things
> 
> opt.conformance <- 0
> Meps <- .Machine $ double.eps
> 
> 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)
+ }
> 
> is.infinite(.Machine$double.base ^ .Machine$double.max.exp)# overflow
[1] TRUE
> 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
> 
> log10(.Machine$double.xmax) / log10(2) == .Machine$double.max.exp
[1] TRUE
> log10(.Machine$double.xmin) / log10(2) == .Machine$double.min.exp
[1] TRUE
> 
> ##  log() and "pow()" -- POSIX is not specific enough..
> log(0) == -Inf
[1] TRUE
> is.nan(log(-1))# TRUE and warning
[1] TRUE
Warning message: 
NaNs produced in: log(x) 
> 
> rp <- c(1:2,Inf); rn <- rev(- rp)
> r <- c(rn, 0, rp)
> all(r^0 == 1)
[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(1^c(-Inf,Inf) == 1)
[1] TRUE
> 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))) <= Meps)
[1] TRUE
> all(abs(1 - x / atan(tan(x))) <= Meps)
[1] TRUE
> 
> ## Complex Trig.:
> abs(Im(cos(acos(1i))) -    1) < 2*Meps
[1] TRUE
> abs(Im(sin(asin(1i))) -    1) < 2*Meps
[1] TRUE
> ##P (1 - Im(sin(asin(Ii))))/Meps
> ##P (1 - Im(cos(acos(Ii))))/Meps
> abs(Im(asin(sin(1i))) -    1) < 2*Meps
[1] TRUE
> cos(1i) == cos(-1i)# i.e. Im(acos(*)) gives + or - 1i:
[1] TRUE
> abs(abs(Im(acos(cos(1i)))) - 1) < 4*Meps
[1] TRUE
> 
> .Random.seed <- c(0, 629, 6137, 22167) # want reproducible output
> Isi <- Im(sin(asin(1i + rnorm(100))))
> all(abs(Isi-1) < 100* Meps)
[1] TRUE
> ##P table(2*abs(Isi-1)    / Meps)
> Isi <- Im(cos(acos(1i + rnorm(100))))
> all(abs(Isi-1) < 100* Meps)
[1] TRUE
> ##P table(2*abs(Isi-1)    / Meps)
> Isi <- Im(atan(tan(1i + rnorm(100)))) #-- tan(atan(..)) does NOT work (Math!)
> all(abs(Isi-1) < 100* Meps)
[1] TRUE
> ##P table(2*abs(Isi-1)    / Meps)
> 
> ## gamma():
> abs(gamma(1/2)^2 - pi) < 4* Meps
[1] TRUE
> r <- rlnorm(5000)
> all(abs(rErr(gamma(r+1), r*gamma(r))) < 500 * 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
> 
> ## 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),
+           tol = 20*Meps))# use tol=0  to see rel.error
[1] TRUE
[1] TRUE
[1] TRUE
[1] TRUE
[1] TRUE
[1] TRUE
[1] TRUE
[1] TRUE
[1] TRUE
>