src/complex/transcendental/cl_C_asinh_aux.cc

   1 // asinh().
   2
   3 // General includes.
   4 #include "cl_sysdep.h"
   5
   6 // Specification.
   7 #include "cl_C.h"
   8
   9
  10 // Implementation.
  11
  12 #include "cln/real.h"
  13 #include "cl_F_tran.h"
  14 #include "cl_R.h"
  15 #include "cln/rational.h"
  16 #include "cl_RA.h"
  17 #include "cln/float.h"
  18
  19 #undef MAYBE_INLINE
  20 #define MAYBE_INLINE inline
  21 #include "cl_F_from_R_def.cc"
  22
  23 namespace cln {
  24
  25 // Hilfsfunktion für asinh und asin: u+iv := arsinh(x+iy). Liefert cl_C_R(u,v).
  26
  27 const cl_C_R asinh (const cl_R& x, const cl_R& y)
  28 {
  29 // Methode:
  30 // Wert und Branch Cuts nach der Formel CLTL2, S. 313:
  31 //   arsinh(z) = log(z+sqrt(1+z^2))
  32 // z=x+iy, Ergebnis u+iv.
  33 // Falls x=0 und y=0: u=0, v=0.
  34 // Falls x=0: arsinh(iy) = i arcsin(y).
  35 //   y rational ->
  36 //     Bei y=1: u = 0, v = pi/2.
  37 //     Bei y=1/2: u = 0, v = pi/6.
  38 //     Bei y=0: u = 0, v = 0.
  39 //     Bei y=-1/2: u = 0, v = -pi/6.
  40 //     Bei y=-1: u = 0, v = -pi/2.
  41 //     Sonst y in Float umwandeln.
  42 //   e := Exponent aus (decode-float y), d := (float-digits y)
  43 //   Bei y=0.0 oder e<=-d/2 liefere u = 0, v = y
  44 //     (denn bei e<=-d/2 ist y^2/3 < y^2/2 < 2^(-d)/2 = 2^(-d-1), also
  45 //     1 <= asin(y)/y < 1+y^2/3 < 1+2^(-d-1) < 1+2^(-d),
  46 //     also ist asin(y)/y, auf d Bits gerundet, gleich 1.0).
  47 //   Berechne 1-y^2.
  48 //   Bei y>1 liefere  u = ln(y+sqrt(y^2-1)), v = pi/2.
  49 //   Bei y<-1 liefere  u = -ln(|y|+sqrt(|y|^2-1)), v = -pi/2.
  50 //   Bei |y|<=1 liefere  u = 0, v = atan(X=sqrt(1-y^2),Y=y).
  51 // Falls y=0:
  52 //   x rational -> x in Float umwandeln.
  53 //   |x|<1/2: u = atanh(x/sqrt(1+x^2)),
  54 //   x>=1/2: u = ln(x+sqrt(1+x^2)),
  55 //   x<=-1/2: u = -ln(-x+sqrt(1+x^2)).
  56 //   v = 0.
  57 // Sonst:
  58 //   z in Bild(sqrt) -> log(sqrt(1+z^2)+z) = (!) = 2 artanh(z/(1+sqrt(1+z^2))).
  59 //   z nicht in Bild(sqrt) ->
  60 //     arsinh(z) = -arsinh(-z).
  61 //     (Denn arsinh(z)+arsinh(-z) == log((z+sqrt(1+z^2))(-z+sqrt(1+z^2)))
  62 //           = log((1+z^2)-z^2) = log(1) = 0 mod 2 pi i, und links ist
  63 //      der Imaginärteil betragsmäßig <=pi.)
  64 //     Also arsinh(z) = -arsinh(-z) = - 2 artanh(-z/(1+sqrt(1+z^2)))
  65 //          = (wegen -artanh(-w) = artanh(w)) = 2 artanh(z/(1+sqrt(1+z^2))).
  66 // Real- und Imaginärteil des Ergebnisses sind Floats, außer wenn z reell oder
  67 // rein imaginär ist.
  68
  69 // Um für zwei Zahlen u,v mit u^2-v^2=1 und u,v beide in Bild(sqrt)
  70 // (d.h. Realteil>0.0 oder Realteil=0.0 und Imaginärteil>=0.0)
  71 // log(u+v) zu berechnen:
  72 //               log(u+v) = 2 artanh(v/(u+1))                            (!)
  73 // (Beweis: 2 artanh(v/(u+1)) = log(1+(v/(u+1))) - log(1-(v/(u+1)))
  74 //  = log((1+u+v)/(u+1)) - log((1+u-v)/(u+1)) == log((1+u+v)/(1+u-v))
  75 //  = log(u+v) mod 2 pi i, und beider Imaginärteil ist > -pi und <= pi.)
  76
  77         if (eq(x,0)) {
  78                 // x=0
  79                 var cl_F yf;
  80                 if (rationalp(y)) {
  81                         DeclareType(cl_RA,y);
  82                         // y rational
  83                         if (eq(y,0)) // x=0, y=0 -> u=0, v=0
  84                                 return cl_C_R(0,0);
  85                         if (integerp(y)) {
  86                                 DeclareType(cl_I,y);
  87                                 // y Integer
  88                                 if (eq(y,1)) // x=0, y=1 -> v = pi/2
  89                                         return cl_C_R(0,scale_float(pi(),-1));
  90                                 if (eq(y,-1)) // x=0, y=-1 -> v = -pi/2
  91                                         return cl_C_R(0,-scale_float(pi(),-1));
  92                                 yf = cl_float(y); // y in Float umwandeln
  93                         } else {
  94                                 DeclareType(cl_RT,y);
  95                                 // y Ratio
  96                                 if (eq(denominator(y),2)) { // Nenner = 2 ?
  97                                         if (eq(numerator(y),1)) // x=0, y=1/2 -> v = pi/6
  98                                                 return cl_C_R(0,pi()/6);
  99                                         if (eq(numerator(y),-1)) // x=0, y=-1/2 -> v = -pi/6
 100                                                 return cl_C_R(0,-(pi()/6));
 101                                 }
 102                                 yf = cl_float(y); // y in Float umwandeln
 103                         }
 104                 } else {
 105                         DeclareType(cl_F,y);
 106                         yf = y;
 107                 }
 108                 // y Float
 109                 var cl_F& y = yf;
 110                 if (zerop(y)) // y=0.0 -> arcsin(y) = y als Ergebnis
 111                         return cl_C_R(0,y);
 112                 if (float_exponent(y) <= (-(sintC)float_digits(y))>>1)
 113                         // e <= -d/2 <==> e <= -ceiling(d/2)
 114                         return cl_C_R(0,y);
 115                 var cl_F temp = 1-square(y);
 116                 if (!minusp(temp))
 117                         // 1-y*y>=0, also |y|<=1
 118                         // v = atan(X=sqrt(1-y*y),Y=y)
 119                         return cl_C_R(0,atan(sqrt(temp),y));
 120                 else {
 121                         // 1-y*y<0, also |y|>1
 122                         temp = sqrt(-temp); // sqrt(y*y-1)
 123                         if (minusp(y))
 124                                 temp = temp - y;
 125                         else
 126                                 temp = temp + y;
 127                         // temp = sqrt(y^2-1)+|y|, ein Float >1
 128                         var cl_F u = ln(temp); // ln(|y|+sqrt(y^2-1)), ein Float >0
 129                         var cl_F v = scale_float(pi(),-1); // (scale-float pi -1) = pi/2
 130                         if (!minusp(y))
 131                                 return cl_C_R(u,v); // y>1 -> v = pi/2
 132                         else
 133                                 return cl_C_R(-u,-v); // y<-1 -> v = -pi/2, u = -ln(...)
 134                 }
 135         }
 136         if (eq(y,0)) {
 137                 // y=0
 138                 var cl_F xf = cl_float(x); // x in Float umwandeln
 139                 var cl_F& x = xf;
 140                 // x Float
 141                 if (zerop(x))
 142                         return cl_C_R(x,0); // x=0.0 -> u=x, v=0.
 143                 var cl_F temp = sqrt(1+square(x)); // sqrt(1+x^2)
 144                 if (float_exponent(x) < 0) // Exponent e (von x/=0) <0 ?
 145                         // |x|<1/2
 146                         return cl_C_R(atanhx(x/temp),0);
 147                 else
 148                         // |x|>=1/2
 149                         if (!minusp(x))
 150                                 // x>=1
 151                                 return cl_C_R(ln(temp+x),0); // u = ln(x+sqrt(1+x^2))
 152                         else
 153                                 // x<=-1
 154                                 return cl_C_R(-ln(temp-x),0); // u = -ln(-x+sqrt(1+x^2))
 155         }
 156         var cl_N z = complex_C(x,y); // z=x+iy
 157         var cl_N w = z/(1+sqrt(1+square(z))); // z/(1+sqrt(1+z^2))
 158         // Da z=x+iy weder reell noch rein imaginär ist, ist auch
 159         // w := z/(1+sqrt(1+z^2)) weder reell noch rein imaginär.
 160         // (Beweis: Sollte sqrt(1+z^2) rationalen Real- und Imaginärteil haben,
 161         // so auch z, also auch w, und die Formel z = 2w/(1-w^2) zeigt, daß dann
 162         // z reell oder rein imaginär sein müßte. Also hat sqrt(1+z^2) ein
 163         // Float als Real- oder Imaginärteil, das Betragsquadrat des Nenners
 164         // ist also ein Float, und da Real- und Imaginärteil von z /=0 sind,
 165         // sind Real- und Imaginärteil von w Floats.)
 166         // Daher hat dann atanh(...) Floats als Realteil u und Imaginärteil v.
 167  {      DeclareType(cl_C,w);
 168         cl_C_R u_v = atanh(realpart(w),imagpart(w));
 169         var cl_R& u = u_v.realpart;
 170         var cl_R& v = u_v.imagpart;
 171   {     DeclareType(cl_F,u);
 172         DeclareType(cl_F,v);
 173         return cl_C_R(scale_float(u,1),scale_float(v,1)); // u:=2*u, v:=2*v
 174 }}}
 175
 176 }  // namespace cln