src/complex/transcendental/cl_C_acosh.cc

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