src/complex/transcendental/cl_C_acosh.cc

   1 // acosh().
   2
   3 // General includes.
   4 #include "cl_sysdep.h"
   5
   6 // Specification.
   7 #include "cl_complex.h"
   8
   9
  10 // Implementation.
  11
  12 #include "cl_C.h"
  13 #include "cl_real.h"
  14 #include "cl_R.h"
  15 #include "cl_rational.h"
  16 #include "cl_RA.h"
  17 #include "cl_float.h"
  18
  19 #undef MAYBE_INLINE
  20 #define MAYBE_INLINE inline
  21 #include "cl_F_from_R_def.cc"
  22
  23 const cl_N acosh (const cl_N& z)
  24 {
  25 // Methode:
  26 // Wert und Branch Cuts nach der Formel CLTL2, S. 314:
  27 //   arcosh(z) = 2 log(sqrt((z+1)/2)+sqrt((z-1)/2))
  28 // Sei z=x+iy.
  29 // Falls y=0:
  30 //   Falls x rational:
  31 //     Bei x=1: Ergebnis 0.
  32 //     Bei x=1/2: Ergebnis pi/3 i.
  33 //     Bei x=0: Ergebnis pi/2 i.
  34 //     Bei x=-1/2: Ergebnis 2pi/3 i.
  35 //     Bei x=-1: Ergebnis pi i.
  36 //   Falls x<-1:
  37 //     x in Float umwandeln, Ergebnis log(sqrt(x^2-1)-x) + i pi.
  38 // Sonst nach (!) mit u = sqrt((z+1)/2) und v = sqrt((z-1)/2) :
  39 // arcosh(z) = 4 artanh(v/(u+1)) = 4 artanh(sqrt((z-1)/2)/(1+sqrt((z+1)/2)))
  40
  41 // Um für zwei Zahlen u,v mit u^2-v^2=1 und u,v beide in Bild(sqrt)
  42 // (d.h. Realteil>0.0 oder Realteil=0.0 und Imaginärteil>=0.0)
  43 // log(u+v) zu berechnen:
  44 //               log(u+v) = 2 artanh(v/(u+1))                            (!)
  45 // (Beweis: 2 artanh(v/(u+1)) = log(1+(v/(u+1))) - log(1-(v/(u+1)))
  46 //  = log((1+u+v)/(u+1)) - log((1+u-v)/(u+1)) == log((1+u+v)/(1+u-v))
  47 //  = log(u+v) mod 2 pi i, und beider Imaginärteil ist > -pi und <= pi.)
  48
  49         cl_C_R u_v;
  50         if (realp(z)) {
  51                 DeclareType(cl_R,z);
  52                 // y=0
  53                 var const cl_R& x = z;
  54                 if (rationalp(x)) {
  55                         DeclareType(cl_RA,x);
  56                         // x rational
  57                         if (integerp(x)) {
  58                                 DeclareType(cl_I,x);
  59                                 // x Integer
  60                                 if (eq(x,0)) // x=0 -> Ergebnis pi/2 i
  61                                         return complex_C(0,scale_float(cl_pi(),-1));
  62                                 if (eq(x,1)) // x=1 -> Ergebnis 0
  63                                         return 0;
  64                                 if (eq(x,-1)) // x=-1 -> Ergebnis pi i
  65                                         return complex_C(0,cl_pi());
  66                         } else {
  67                                 DeclareType(cl_RT,x);
  68                                 // x Ratio
  69                                 if (eq(denominator(x),2)) { // Nenner = 2 ?
  70                                         if (eq(numerator(x),1)) // x=1/2 -> Ergebnis pi/3 i
  71                                                 return complex_C(0,cl_pi()/3);
  72                                         if (eq(numerator(x),-1)) // x=-1/2 -> Ergebnis 2pi/3 i
  73                                                 return complex_C(0,scale_float(cl_pi(),1)/3);
  74                                 }
  75                         }
  76                 }
  77                 if (x < cl_I(-1)) {
  78                         // x < -1
  79                         var cl_F xf = cl_float(x);
  80                         var cl_F& x = xf;
  81                         // x Float <= -1
  82                         // log(sqrt(x^2-1)-x), ein Float >=0, Imaginärteil pi
  83                         return complex_C(ln(sqrt(square(x)-1)-x),cl_pi());
  84                 }
  85         }
  86         return 4 * atanh( sqrt(minus1(z)/2) / plus1(sqrt(plus1(z)/2)) );
  87 }