[Rd] branch cuts in atan(), asin(), acos() (PR#7871)
r.hankin at noc.soton.ac.uk
r.hankin at noc.soton.ac.uk
Mon May 16 15:48:48 CEST 2005
The complex versions of atain(), asin(), acos() are not documented
except in the source
code.
The branch cuts that are implemented do not tally with those of
Abramowitz and Stegun,
figure 4.4, p79. The current branch cuts appear to be inherited from
sqrt() and log(),
in a non-obvious manner.
Also, the current behaviour of atan() leads to the following:
R> tan(atan(2))
[1] 2
R> tan(atan(2+0i))
[1] -0.5+0i
(I would expect a value close to 2 in both these)
R> atan(2)
[1] 1.107149
R> atan(2+0i)
[1] -0.4636476+0i
>
(I would expect these two expressions to evaluate to numbers with small
absolute
difference)
> atan(1.0001+0i)
[1] -0.7853482+0i
> atan(0.9999+0i)
[1] 0.7853482+0i
>
(I would expect these two expressions to evaluate to numbers with small
absolute
difference)
The following patch to src/main/complex.c
appears to implement the desired behaviour:<
octopus:~/downloads/R-2.1.0/src/main% diff -c complex.c
/Users/rksh/unix/rksh/complex.c
*** complex.c Mon Apr 18 22:30:00 2005
--- /Users/rksh/unix/rksh/complex.c Mon May 16 14:33:15 2005
***************
*** 367,372 ****
--- 367,376 ----
bet = t1 - t2;
r->r = asin(bet);
r->i = log(alpha + sqrt(alpha*alpha - 1));
+
+ if(y<0){
+ r->i *= -1;
+ }
}
static void z_acos(Rcomplex *r, Rcomplex *z)
***************
*** 388,393 ****
--- 392,404 ----
r->r = 0.5 * atan(2 * x / ( 1 - x * x - y * y));
r->i = 0.25 * log((x * x + (y + 1) * (y + 1)) /
(x * x + (y - 1) * (y - 1)));
+
+ if((x*x+y*y > 1) && x>0){
+ r->r += M_PI_2;
+ }
+ if((x*x+y*y > 1) && x<0){
+ r->r -= M_PI_2;
+ }
}
static void z_atan2(Rcomplex *r, Rcomplex *csn, Rcomplex *ccs)
octopus:~/downloads/R-2.1.0/src/main%
--
Robin Hankin
Uncertainty Analyst
National Oceanography Centre, Southampton
European Way, Southampton SO14 3ZH, UK
tel 023-8059-7743
More information about the R-devel
mailing list