The R Project SVN R

Rev

Rev 54435 | Rev 69410 | Go to most recent revision | Details | Compare with Previous | Last modification | View Log | RSS feed

Rev Author Line No. Line
35226 ripley 1 
### Tests of complex arithemetic.
2 
 
49068 maechler 3 
Meps <- .Machine$double.eps
35226 ripley 4 
## complex
5 
z <- 0i ^ (-3:3)
6 
stopifnot(Re(z) == 0 ^ (-3:3))
7 
 
8 
 
9 
## powers, including complex ones
10 
a <- -4:12
11 
m <- outer(a +0i, b <- seq(-.5,2, by=.5), "^")
12 
dimnames(m) <- list(paste(a), "^" = sapply(b,format))
13 
round(m,3)
49068 maechler 14 
stopifnot(m[,as.character(0:2)] == cbind(1,a,a*a),
15 
                                        # latter were only approximate
16 
          all.equal(unname(m[,"0.5"]),
17 
                    sqrt(abs(a))*ifelse(a < 0, 1i, 1),
18 
                    tol= 20*Meps))
54249 ripley 19 
 
20 
## 2.10.0-2.12.1 got z^n wrong in the !HAVE_C99_COMPLEX case
21 
z <- 0.2853725+0.3927816i
22 
z2 <- z^(1:20)
23 
z3 <- z^-(1:20)
24 
z0 <- cumprod(rep(z, 20))
25 
stopifnot(all.equal(z2, z0), all.equal(z3, 1/z0))
26 
## was z^3 had value z^2 ....
27 
 
35253 ripley 28 
## fft():
29 
for(n in 1:30) cat("\nn=",n,":", round(fft(1:n), 8),"\n")
35226 ripley 30 
 
35253 ripley 31 
 
54249 ripley 32 
## polyroot():
33 
stopifnot(abs(1 + polyroot(choose(8, 0:8))) < 1e-10)# maybe smaller..
34 
 
35 
## precision of complex numbers
36 
signif(1.678932e80+0i, 5)
37 
signif(1.678932e-300+0i, 5)
38 
signif(1.678932e-302+0i, 5)
39 
signif(1.678932e-303+0i, 5)
40 
signif(1.678932e-304+0i, 5)
41 
signif(1.678932e-305+0i, 5)
42 
signif(1.678932e-306+0i, 5)
43 
signif(1.678932e-307+0i, 5)
44 
signif(1.678932e-308+0i, 5)
45 
signif(1.678932-1.238276i, 5)
46 
signif(1.678932-1.238276e-1i, 5)
47 
signif(1.678932-1.238276e-2i, 5)
48 
signif(1.678932-1.238276e-3i, 5)
49 
signif(1.678932-1.238276e-4i, 5)
50 
signif(1.678932-1.238276e-5i, 5)
51 
signif(8.678932-9.238276i, 5)
52 
## prior to 2.2.0 rounded real and imaginary parts separately.
53 
 
54 
 
35226 ripley 55 
## Complex Trig.:
56 
abs(Im(cos(acos(1i))) -	 1) < 2*Meps
57 
abs(Im(sin(asin(1i))) -	 1) < 2*Meps
58 
##P (1 - Im(sin(asin(Ii))))/Meps
59 
##P (1 - Im(cos(acos(Ii))))/Meps
60 
abs(Im(asin(sin(1i))) -	 1) < 2*Meps
57524 ripley 61 
all.equal(cos(1i), cos(-1i)) # i.e. Im(acos(*)) gives + or - 1i:
35226 ripley 62 
abs(abs(Im(acos(cos(1i)))) - 1) < 4*Meps
63 
 
64 
 
42034 ripley 65 
set.seed(123) # want reproducible output
35226 ripley 66 
Isi <- Im(sin(asin(1i + rnorm(100))))
67 
all(abs(Isi-1) < 100* Meps)
68 
##P table(2*abs(Isi-1)	/ Meps)
69 
Isi <- Im(cos(acos(1i + rnorm(100))))
70 
all(abs(Isi-1) < 100* Meps)
71 
##P table(2*abs(Isi-1)	/ Meps)
72 
Isi <- Im(atan(tan(1i + rnorm(100)))) #-- tan(atan(..)) does NOT work (Math!)
73 
all(abs(Isi-1) < 100* Meps)
74 
##P table(2*abs(Isi-1)	/ Meps)
75 
 
54249 ripley 76 
set.seed(123)
77 
z <- complex(real = rnorm(100), imag = rnorm(100))
78 
stopifnot(Mod ( 1 -  sin(z) / ( (exp(1i*z)-exp(-1i*z))/(2*1i) )) < 20 * Meps)
79 
## end of moved from complex.Rd
35226 ripley 80 
 
81 
 
82 
## PR#7781
83 
## This is not as given by e.g. glibc on AMD64
84 
(z <- tan(1+1000i)) # 0+1i from R's own code.
85 
stopifnot(is.finite(z))
86 
##
87 
 
88 
 
89 
## Branch cuts in complex inverse trig functions
90 
atan(2)
91 
atan(2+0i)
92 
tan(atan(2+0i))
35227 ripley 93 
## should not expect exactly 0i in result
94 
round(atan(1.0001+0i), 7)
95 
round(atan(0.9999+0i), 7)
35226 ripley 96 
## previously not as in Abramowitz & Stegun.
97 
 
98 
 
99 
## typo in z_atan2.
100 
(z <- atan2(0+1i, 0+0i))
101 
stopifnot(all.equal(z, pi/2+0i))
102 
## was NA in 2.1.1
35250 ripley 103 
 
104 
 
54249 ripley 105 
## Hyperbolic
106 
x <- seq(-3, 3, len=200)
107 
Meps <- .Machine$double.eps
108 
stopifnot(
109 
 Mod(cosh(x) - cos(1i*x))	< 20*Meps,
54435 ripley 110 
 Mod(sinh(x) - sin(1i*x)/1i)	< 20*Meps
54249 ripley 111 
)
112 
## end of moved from Hyperbolic.Rd
54158 ripley 113 
 
54262 ripley 114 
## values near and on branch cuts
115 
options(digits=5)
116 
z <- c(2+0i, 2-0.0001i, -2+0i, -2+0.0001i)
117 
asin(z)
118 
acos(z)
119 
atanh(z)
120 
z <- c(0+2i, 0.0001+2i, 0-2i, -0.0001i-2i)
121 
asinh(z)
122 
acosh(z)
123 
atan(z)
54435 ripley 124 
## According to C99, should have continuity from the side given if there
125 
## are not signed zeros
54262 ripley 126 
## Both glibc 2.12 and Mac OS X 10.6 use continuity from above in the first set
54435 ripley 127 
## but they seem to assume signed zeros.
54262 ripley 128 
## Windows gave incorrect (NaN) values on the cuts.
129