# wch/r-source

### Subversion checkout URL

You can clone with
or
.
Fetching contributors…

Cannot retrieve contributors at this time

147 lines (121 sloc) 4.45 kB
 ####=== 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) } abs(1- .Machine\$double.xmin * 10^(-.Machine\$double.min.exp*log10(2)))/Meps < 1e3 ##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) is.infinite( i1) & is.infinite( i2) & i1 > 12 & i2 > 12 is.infinite(-i1) & is.infinite(-i2) & (-i1) < -12 & (-i2) < -12 is.nan(n1 <- 0 / 0) is.nan( - n1) i1 == i1 + i1 i1 == i1 * i1 is.nan(i1 - i1) is.nan(i1 / i1) 1/0 == Inf & 0 ^ -1 == Inf 1/Inf == 0 & Inf ^ -1 == 0 !is.na(Inf) & !is.nan(Inf) & is.infinite(Inf) & !is.finite(Inf) !is.na(-Inf)& !is.nan(-Inf)& is.infinite(-Inf)& !is.finite(-Inf) is.na(NA) & !is.nan(NA) & !is.infinite(NA) & !is.finite(NA) is.na(NaN) & is.nan(NaN) & !is.infinite(NaN) & !is.finite(NaN) all(!is.nan(c(1,NA))) all(c(F,T,F) == is.nan(c (1,NaN,NA))) all(c(F,T,F) == is.nan(list(1,NaN,NA))) ## log() and "pow()" -- POSIX is not specific enough.. log(0) == -Inf is.nan(log(-1))# TRUE and warning rp <- c(1:2,Inf); rn <- rev(- rp) r <- c(rn, 0, rp) all(r^0 == 1) all((rn ^ -3) == -((-rn) ^ -3)) # all(c(1.1,2,Inf) ^ Inf == Inf) all(c(1.1,2,Inf) ^ -Inf == 0) .9 ^ Inf == 0 .9 ^ -Inf == Inf ## Wasn't ok in 0.64: all(1^c(-Inf,Inf) == 1) all(is.nan(rn ^ .5))# in some C's : (-Inf) ^ .5 gives Inf, instead of NaN ## Real Trig.: cos(0) == 1 sin(3*pi/2) == cos(pi) x <- rnorm(99) all( sin(-x) == - sin(x)) all( cos(-x) == cos(x)) x <- 1:99/100 all(abs(1 - x / asin(sin(x))) <= Meps) all(abs(1 - x / atan(tan(x))) <= Meps) ## Complex Trig.: abs(Im(cos(acos(1i))) - 1) < 2*Meps abs(Im(sin(asin(1i))) - 1) < 2*Meps ##P (1 - Im(sin(asin(Ii))))/Meps ##P (1 - Im(cos(acos(Ii))))/Meps abs(Im(asin(sin(1i))) - 1) < 2*Meps cos(1i) == cos(-1i)# i.e. Im(acos(*)) gives + or - 1i: abs(abs(Im(acos(cos(1i)))) - 1) < 4*Meps .Random.seed <- c(0, 629, 6137, 22167) # want reproducible output Isi <- Im(sin(asin(1i + rnorm(100)))) all(abs(Isi-1) < 100* Meps) ##P table(2*abs(Isi-1) / Meps) Isi <- Im(cos(acos(1i + rnorm(100)))) all(abs(Isi-1) < 100* Meps) ##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) ##P table(2*abs(Isi-1) / Meps) ## gamma(): abs(gamma(1/2)^2 - pi) < 4* Meps r <- rlnorm(5000) all(abs(rErr(gamma(r+1), r*gamma(r))) < 500 * Meps) n <- 10; all( gamma(1:n) == cumprod(c(1,1:(n-1)))) n <- 20; all(abs(rErr( gamma(1:n), cumprod(c(1,1:(n-1))))) < 100*Meps) n <- 120; all(abs(rErr( gamma(1:n), cumprod(c(1,1:(n-1))))) < 1000*Meps) n <- 10000;all(abs(rErr(lgamma(1:n),cumsum(log(c(1,1:(n-1)))))) < 100*Meps) ## 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 ## 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 ## 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)) all( pmin(v1, 1:3) == pmin(1:3, v1) & pmin(1:3, v1) == c(1,2,2)) 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)) all( pmin(m1, 1:7) == pmin(1:7, m1) & pmin(1:7, m1) == c(1:3,3,3,3,2)) all( pmax(m2, 1:7) == pmax(1:7, m2) & pmax(1:7, m2) == pmax(1:7, m1)) all( pmin(m2, 1:7) == pmin(1:7, m2) & pmin(1:7, m2) == c(1:3,2,2,2,2)) options(oo)
Something went wrong with that request. Please try again.