src/complex/transcendental/cl_C_atanh.cc

   1 // atanh().
   2
   3 // General includes.
   4 #include "base/cl_sysdep.h"
   5
   6 // Specification.
   7 #include "cln/complex.h"
   8
   9
  10 // Implementation.
  11
  12 #include "complex/cl_C.h"
  13 #include "cln/real.h"
  14
  15 namespace cln {
  16
  17 // Methode:
  18 // Wert und Branch Cuts nach der Formel CLTL2, S. 315:
  19 //   artanh(z) = (log(1+z)-log(1-z)) / 2
  20 // Sei z=x+iy, Ergebnis u+iv.
  21 // Falls x=0 und y=0: u=0, v=0.
  22 // Falls x=0: u = 0, v = atan(X=1,Y=y).
  23 // Falls y=0:
  24 //   x rational -> x in Float umwandeln.
  25 //   |x|<1/2: u = atanh(x), v = 0.
  26 //   |x|>=1/2: (1+x)/(1-x) errechnen,
  27 //             =0 -> Error,
  28 //             >0 (also |x|<1) -> u = 1/2 log((1+x)/(1-x)), v = 0.
  29 //             <0 (also |x|>1) -> u = 1/2 log(-(1+x)/(1-x)),
  30 //                                v = (-pi/2 für x>1, pi/2 für x<-1).
  31 // Sonst:
  32 //   1+x und 1-x errechnen.
  33 //   x und y in Floats umwandeln.
  34 //   |4x| und 1+x^2+y^2 errechnen,
  35 //   |4x| < 1+x^2+y^2 -> u = 1/2 atanh(2x/(1+x^2+y^2)),
  36 //   |4x| >= 1+x^2+y^2 -> u = 1/4 ln ((1+x^2+y^2)+2x)/((1+x^2+y^2)-2x)
  37 //                        oder besser (an der Singularität: |x|-1,|y| klein):
  38 //                        u = 1/4 ln ((1+x)^2+y^2)/((1-x)^2+y^2).
  39 //   v = 1/2 atan(X=(1-x)(1+x)-y^2,Y=2y) * (-1 falls Y=0.0 und X<0.0 und x>=0.0,
  40 //                                          1 sonst)
  41 // Ergebnis ist reell nur, wenn z reell.
  42 // Real- und Imaginärteil des Ergebnisses sind Floats, außer wenn z reell oder
  43 // rein imaginär ist.
  44
  45 inline const cl_C_R _atanh (const cl_N& z)
  46 {
  47         if (realp(z)) {
  48                 DeclareType(cl_R,z);
  49                 return atanh(z,0);
  50         } else {
  51                 DeclareType(cl_C,z);
  52                 return atanh(realpart(z),imagpart(z));
  53         }
  54 }
  55
  56 const cl_N atanh (const cl_N& z)
  57 {
  58         var cl_C_R u_v = _atanh(z);
  59         var cl_R& u = u_v.realpart;
  60         var cl_R& v = u_v.imagpart;
  61         return complex(u,v);
  62 }
  63
  64 }  // namespace cln