The R Project SVN R

Rev

Rev 79080 | Show entire file | Ignore whitespace | Details | Blame | Last modification | View Log | RSS feed

Rev 79080 Rev 89335
Line 81... Line 81...
81 
 
 81 
 
82 
 
 82 
 
83 
## Real Trig.:
 83 
## Real Trig.:
84 
cos(0) == 1
 84 
cos(0) == 1
85 
sin(3*pi/2) == cos(pi)
 85 
sin(3*pi/2) == cos(pi)
86 
x <- rnorm(99)
 - 
 
87 
all( sin(-x) == - sin(x))
 - 
 
88 
all( cos(-x) == cos(x))
 - 
 
89 
 
 86 
 
90 
x <- 1:99/100
 87 
x <- 1:99/100
- 
 
 88 
## IGNORE_RDIFF_BEGIN
91 
all(abs(1 - x / asin(sin(x))) <= 2*Meps)# "== 2*" for HP-UX
 89 
table(Isin <- abs(1 - x / asin(sin(x)))/Meps)# mostly 0, often <= 1
92 
all(abs(1 - x / atan(tan(x))) <  2*Meps)
 90 
table(Itan <- abs(1 - x / atan(tan(x)))/Meps)#  (ditto)
- 
 
 91 
## IGNORE_RDIFF_END
- 
 
 92 
all(Isin <= 4) # typically <= 1, 3 needed for [mingw-w64 v12 on Win]
- 
 
 93 
all(Itan <= 4)
93 
 
 94 
 
94 
## Sun has asin(.) = acos(.) = 0 for these:
 95 
## Sun has asin(.) = acos(.) = 0 for these:
95 
##	is.nan(acos(1.1)) && is.nan(asin(-2)) [!]
 96 
##	is.nan(acos(1.1)) && is.nan(asin(-2)) [!]
96 
 
 97 
 
97 
## gamma()
 98 
## gamma()
Line 132... Line 133...
132 
          pi^(2*(1:3)) * c(1/6, 1/15, 8/63), tolerance = 32*Meps)
 133 
          pi^(2*(1:3)) * c(1/6, 1/15, 8/63), tolerance = 32*Meps)
133 
x <- c(-100,-3:2, -99.9, -7.7, seq(-3,3, length=61), 5.1, 77)
 134 
x <- c(-100,-3:2, -99.9, -7.7, seq(-3,3, length=61), 5.1, 77)
134 
## Intel icc showed a < 1ulp difference in the second.
 135 
## Intel icc showed a < 1ulp difference in the second.
135 
stopifnot(all.equal( digamma(x), psigamma(x,0), tolerance = 2*Meps),
 136 
stopifnot(all.equal( digamma(x), psigamma(x,0), tolerance = 2*Meps),
136 
          all.equal(trigamma(x), psigamma(x,1), tolerance = 2*Meps))# TRUE (+ NaN warnings)
 137 
          all.equal(trigamma(x), psigamma(x,1), tolerance = 2*Meps))# TRUE (+ NaN warnings)
137 
## very large x:
 138 
## very large x,  compare with asymptotic approximation:
138 
x <- 1e30 ^ (1:10)
 139 
x <- 1e30 ^ (1:10)
139 
a.relE <- function(appr, true) abs(1 - appr/true)
 140 
arelEps <- function(appr, true) abs(1 - appr/true) / Meps
- 
 
 141 
## IGNORE_RDIFF_BEGIN
140 
stopifnot(a.relE(digamma(x),   log(x)) < 1e-13,
 142 
sort( digamE <- arelEps( digamma(x), log(x))) # all 0 for glibc
141 
          a.relE(trigamma(x),     1/x) < 1e-13)
 143 
sort(trigamE <- arelEps(trigamma(x),  1/x  )) # 11 22 ... 108 159 for glibc
142 
x <- sqrt(x[2:6]); stopifnot(a.relE(psigamma(x,2), - 1/x^2) < 1e-13)
 144 
x <- sqrt(x[2:6]); sort(psiga2E <- arelEps(psigamma(x,2), - 1/x^2)) # 22 33 43 55 65
143 
x <- 10^(10*(2:6));stopifnot(a.relE(psigamma(x,5), +24/x^5) < 1e-13)
 145 
x <- 10^(10*(2:6));sort(psiga5E <- arelEps(psigamma(x,5), +24/x^5)) #  5 50.5 .. 106
- 
 
 146 
## IGNORE_RDIFF_END
- 
 
 147 
stopifnot(digamE  <   64, # seen '0'
- 
 
 148 
          trigamE < 1300, # mfuzz: upto 671
- 
 
 149 
          psiga2E <  430, # mfuzz: upto 212
- 
 
 150 
          psiga5E < 1010) # mfuzz: upto 505.5
144 
 
 151 
 
145 
## fft():
 152 
## fft():
146 
ok <- TRUE
 153 
ok <- TRUE
147 
##test EXTENSIVELY:	for(N in 1:100) {
 154 
##test EXTENSIVELY:	for(N in 1:100) {
148 
    cat(".")
 155 
    cat(".")
Line 159... Line 166...
159 
ok
 166 
ok
160 
 
 167 
 
161 
## var():
 168 
## var():
162 
for(n in 2:10)
 169 
for(n in 2:10)
163 
    print(all.equal(n*(n-1)*var(diag(n)),
 170 
    print(all.equal(n*(n-1)*var(diag(n)),
164 
		    matrix(c(rep(c(n-1,rep(-1,n)),n-1), n-1), nr=n, nc=n),
 171 
		    matrix(c(rep(c(n-1,rep(-1,n)),n-1), n-1), nrow=n, ncol=n),
- 
 
 172 
                                ##----- use 0*Meps to see non-0 rel.error ==> all TRUE for glibc
165 
		    tolerance = 20*Meps)) # use tolerance = 0 to see rel.error
 173 
		    tolerance = 20*Meps))
166 
 
 174 
 
167 
## pmin() & pmax() -- "attributes" !
 175 
## pmin() & pmax() -- "attributes" !
168 
v1 <- c(a=2)
 176 
v1 <- c(a=2)
169 
m1 <- cbind(  2:4,3)
 177 
m1 <- cbind(  2:4,3)
170 
m2 <- cbind(a=2:4,2)
 178 
m2 <- cbind(a=2:4,2)