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

______________________________________________
R-devel@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel

Reply via email to