src/complex/transcendental/cl_C_asinh.cc

   1 // asinh().
   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
  15 // Methode:
  16 // Wert und Branch Cuts nach der Formel CLTL2, S. 313:
  17 //   arsinh(z) = log(z+sqrt(1+z^2))
  18 // z=x+iy, Ergebnis u+iv.
  19 // Falls x=0 und y=0: u=0, v=0.
  20 // Falls x=0: arsinh(iy) = i arcsin(y).
  21 //   y rational ->
  22 //     Bei y=1: u = 0, v = pi/2.
  23 //     Bei y=1/2: u = 0, v = pi/6.
  24 //     Bei y=0: u = 0, v = 0.
  25 //     Bei y=-1/2: u = 0, v = -pi/6.
  26 //     Bei y=-1: u = 0, v = -pi/2.
  27 //     Sonst y in Float umwandeln.
  28 //   e := Exponent aus (decode-float y), d := (float-digits y)
  29 //   Bei y=0.0 oder e<=-d/2 liefere u = 0, v = y
  30 //     (denn bei e<=-d/2 ist y^2/3 < y^2/2 < 2^(-d)/2 = 2^(-d-1), also
  31 //     1 <= asin(y)/y < 1+y^2/3 < 1+2^(-d-1) < 1+2^(-d),
  32 //     also ist asin(y)/y, auf d Bits gerundet, gleich 1.0).
  33 //   Berechne 1-y^2.
  34 //   Bei y>1 liefere  u = ln(y+sqrt(y^2-1)), v = pi/2.
  35 //   Bei y<-1 liefere  u = -ln(|y|+sqrt(|y|^2-1)), v = -pi/2.
  36 //   Bei |y|<=1 liefere  u = 0, v = atan(X=sqrt(1-y^2),Y=y).
  37 // Falls y=0:
  38 //   x rational -> x in Float umwandeln.
  39 //   |x|<1/2: u = atanh(x/sqrt(1+x^2)),
  40 //   x>=1/2: u = ln(x+sqrt(1+x^2)),
  41 //   x<=-1/2: u = -ln(-x+sqrt(1+x^2)).
  42 //   v = 0.
  43 // Sonst:
  44 //   z in Bild(sqrt) -> log(sqrt(1+z^2)+z) = (!) = 2 artanh(z/(1+sqrt(1+z^2))).
  45 //   z nicht in Bild(sqrt) ->
  46 //     arsinh(z) = -arsinh(-z).
  47 //     (Denn arsinh(z)+arsinh(-z) == log((z+sqrt(1+z^2))(-z+sqrt(1+z^2)))
  48 //           = log((1+z^2)-z^2) = log(1) = 0 mod 2 pi i, und links ist
  49 //      der Imaginärteil betragsmäßig <=pi.)
  50 //     Also arsinh(z) = -arsinh(-z) = - 2 artanh(-z/(1+sqrt(1+z^2)))
  51 //          = (wegen -artanh(-w) = artanh(w)) = 2 artanh(z/(1+sqrt(1+z^2))).
  52 // Real- und Imaginärteil des Ergebnisses sind Floats, außer wenn z reell oder
  53 // rein imaginär ist.
  54
  55 // Um für zwei Zahlen u,v mit u^2-v^2=1 und u,v beide in Bild(sqrt)
  56 // (d.h. Realteil>0.0 oder Realteil=0.0 und Imaginärteil>=0.0)
  57 // log(u+v) zu berechnen:
  58 //               log(u+v) = 2 artanh(v/(u+1))                            (!)
  59 // (Beweis: 2 artanh(v/(u+1)) = log(1+(v/(u+1))) - log(1-(v/(u+1)))
  60 //  = log((1+u+v)/(u+1)) - log((1+u-v)/(u+1)) == log((1+u+v)/(1+u-v))
  61 //  = log(u+v) mod 2 pi i, und beider Imaginärteil ist > -pi und <= pi.)
  62
  63 inline const cl_C_R _asinh (const cl_N& z)
  64 {
  65         if (realp(z)) {
  66                 DeclareType(cl_R,z);
  67                 return asinh(z,0);
  68         } else {
  69                 DeclareType(cl_C,z);
  70                 return asinh(realpart(z),imagpart(z));
  71         }
  72 }
  73
  74 const cl_N asinh (const cl_N& z)
  75 {
  76         var cl_C_R u_v = _asinh(z);
  77         var cl_R& u = u_v.realpart;
  78         var cl_R& v = u_v.imagpart;
  79         return complex(u,v);
  80 }