r.hankin@noc.soton.ac.uk
2005-May-16 15:48 UTC
[Rd] branch cuts in atan(), asin(), acos() (PR#7871)
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
ripley@stats.ox.ac.uk
2005-May-17 08:02 UTC
[Rd] branch cuts in atan(), asin(), acos() (PR#7871)
On Mon, 16 May 2005 r.hankin@noc.soton.ac.uk wrote:> The complex versions of atain(), asin(), acos() are not documented > except in the source > code.Thank you, but your patch does not address that so they would become undocumented (they no longer follow the comment in the source code, which BTW is regarded as documentation). Could you please supply documentation to match this change, including the exact definition you have implemented. Also: have you looked at atan2?> 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 > >-- Brian D. Ripley, ripley@stats.ox.ac.uk Professor of Applied Statistics, http://www.stats.ox.ac.uk/~ripley/ University of Oxford, Tel: +44 1865 272861 (self) 1 South Parks Road, +44 1865 272866 (PA) Oxford OX1 3TG, UK Fax: +44 1865 272595
r.hankin@noc.soton.ac.uk
2005-May-17 15:16 UTC
[Rd] branch cuts in atan(), asin(), acos() (PR#7871)
On May 17, 2005, at 07:01 am, Prof Brian Ripley wrote:> On Mon, 16 May 2005 r.hankin@noc.soton.ac.uk wrote: > >> The complex versions of atain(), asin(), acos() are not documented >> except in the source >> code. > > Thank you, but your patch does not address that so they would become > undocumented (they no longer follow the comment in the source code, > which BTW is regarded as documentation). > > Could you please supply documentation to match this change, including > the exact definition you have implemented. > > Also: have you looked at atan2? >OK, this follows PR#7871, discussing branch cuts of inverse trig functions. In this email, I provide a context diff between complex.c (R-2.1.0) and a version that moves the branch cuts to the standard places. The changes are documented in the code. I also give updated Trig.Rd and Hyperbolic.Rd files which document the behaviour of the new functions. R-2.1.0 passes "make check" with the new files (after editing Rdutils.Rd, which caused a problem). I have not attempted to alter z_atan2(): there appear to be unrelated issues here. For example, atan2(0,0) gives 0 as documented, but atan2(0i,0i) gives NA. Also, atan2(0,1i) returns an error (I would expect a zero here too). I feel it is better to consider z_atan2() as a separate issue. I will file a separate bugreport for this. diff -c <old complex.c> <new complex.c> *** /Users/rksh/downloads/R-2.1.0/src/main/complex.c Mon Apr 18 22:30:00 2005 --- complex.c Tue May 17 13:35:27 2005 *************** *** 355,360 **** --- 355,365 ---- /* Complex Arcsin and Arccos Functions */ /* Equation (4.4.37) Abramowitz and Stegun */ + /* with additional terms to force the branch + /* to agree with figure 4.4, p79. Continuity */ + /* on the branch cuts (real axis; y==0, |x| > 1) is */ + /* standard: z_asin() is continuous from below if x >= 1 */ + /* and continuous from above if x <= -1. */ static void z_asin(Rcomplex *r, Rcomplex *z) { *************** *** 367,372 **** --- 372,382 ---- bet = t1 - t2; r->r = asin(bet); r->i = log(alpha + sqrt(alpha*alpha - 1)); + + if(y<0 || (y==0 && x > 1)){ + r->i *= -1; + } + } static void z_acos(Rcomplex *r, Rcomplex *z) *************** *** 379,384 **** --- 389,399 ---- /* Complex Arctangent Function */ /* Equation (4.4.39) Abramowitz and Stegun */ + /* with additional terms to force the branch cuts */ + /* to agree with figure 4.4, p79. Continuity */ + /* on the branch cuts (pure imaginary axis; x==0, |y|>1) */ + /* is standard: z_asin() is continuous from the right */ + /* if y >= 1, and continuous from the left if y <= -1. */ static void z_atan(Rcomplex *r, Rcomplex *z) { *************** *** 388,393 **** --- 403,416 ---- 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){ + r->r += M_PI_2; + if(x < 0 || (x==0 && y<0)){ + r->r -= M_PI; + } + } + } static void z_atan2(Rcomplex *r, Rcomplex *csn, Rcomplex *ccs) ====FILE Hyperbolic.Rd STARTS===\name{Hyperbolic} \title{Hyperbolic Functions} \usage{ cosh(x) sinh(x) tanh(x) acosh(x) asinh(x) atanh(x) } \alias{cosh} \alias{sinh} \alias{tanh} \alias{acosh} \alias{asinh} \alias{atanh} \description{ These functions give the obvious hyperbolic functions. They respectively compute the hyperbolic cosine, sine, tangent, and their inverses, arc-cosine, arc-sine, arc-tangent (or \dQuote{\emph{area cosine}}, etc). } \arguments{ \item{x}{a numeric or complex vector} } \details{ These are generic functions: methods can be defined for them individually or via the \code{\link{Math}} group generic. Branch cuts are consistent with the inverse trigonometric functions \code{asin()} et seq, and agree with those defined in Abramowitz and Stegun, figure 4.7, page 86. } \seealso{ The trigonometric functions, \code{\link{cos}}, \code{\link{sin}}, \code{\link{tan}}, and their inverses \code{\link{acos}}, \code{\link{asin}}, \code{\link{atan}}. The logistic distribution function \code{\link{plogis}} is a shifted version of \code{tanh()}. } \keyword{math} ====FILE Hyperbolic.Rd ENDS=== ====File Trig.Rd STARTS===\name{Trig} \alias{Trig} \alias{cos} \alias{sin} \alias{tan} \alias{acos} \alias{asin} \alias{atan} \alias{atan2} \title{Trigonometric Functions} \description{ These functions give the obvious trigonometric functions. They respectively compute the cosine, sine, tangent, arc-cosine, arc-sine, arc-tangent, and the two-argument arc-tangent. } \usage{ cos(x) sin(x) tan(x) acos(x) asin(x) atan(x) atan2(y, x) } \arguments{ \item{x, y}{numeric or complex vector} } \details{ The arc-tangent of two arguments \code{atan2(y,x)} returns the angle between the x-axis and the vector from the origin to \eqn{(x,y)}, i.e., for positive arguments \code{atan2(y,x) == atan(y/x)}. Angles are in radians, not degrees (i.e., a right angle is \eqn{\pi/2}). All except \code{atan2} are generic functions: methods can be defined for them individually or via the \code{\link{Math}} group generic. For the inverse trigonometric functions, branch cuts are defined as in Abramowitz and Stegun, figure 4.4, page 79. Continuity on the branch cuts is standard. \item For \code{asin()} and \code{acos()}, there are two cuts, both along the real axis: \eqn{\left(-\infty,-1\right]}{\(-infinity,1\]} and \eqn{\left[1,\infty\right)}{\[1,infinity\)}. Functions \code{asin()} and \code{acos()} are continuous from above on the interval \eqn{\left(-\infty,-1\right]}{\(-infinity,1\]} and continuous from below on \eqn{\left[1,\infty\right)}{\[1,infinity\)}. \item For \code{atan()} there are two cuts, both along the pure imaginary axis: \eqn{\left(-i\infty,-i\right]}{\(-i*infinity,-i\]} and \eqn{\left[i,i\infty\right)}{\[i,i*infinity\)}. Function \code{atan()} is continuous from the left on the interval \eqn{\left(-i\infty,-i\right]}{\(-i*infinity,-i\]} and from the right on the interval \eqn{\left[i,i\infty\right)}{\[i,i*infinity\)}. } \references{ Becker, R. A., Chambers, J. M. and Wilks, A. R. (1988) \emph{The New S Language}. Wadsworth \& Brooks/Cole. Abramowitz, M. and Stegun, I. (1964). \emph{Handbook of Mathematical Functions}, Dover. } \examples{ #Visualize branch cuts: x <- seq(from = -4, to = 4, len = 100) z <- outer(x,1i*x,"+") persp(x,x,Im(asin(z)),theta= 40,phi=40,r=1e9, expand=0.6, ticktype="detailed") persp(x,x,Re(atan(z)),theta=-40,phi=40,r=1e9, expand=0.6, ticktype="detailed") #check continuity on branch cuts: asin( 3+ 1e-9i*(-1:1)) asin(-3+ 1e-9i*(-1:1)) atan( 3i + 1e-9*(-1:1)) atan(-3i + 1e-9*(-1:1)) } \keyword{math} ====FILE Trig.Rd ENDS=== -- Robin Hankin Uncertainty Analyst National Oceanography Centre, Southampton European Way, Southampton SO14 3ZH, UK tel 023-8059-7743