src/complex/transcendental/cl_C_atanh_aux.cc

   1 // atanh().
   2
   3 // General includes.
   4 #include "cl_sysdep.h"
   5
   6 // Specification.
   7 #include "cl_C.h"
   8
   9
  10 // Implementation.
  11
  12 #include "cl_N.h"
  13 #include "cln/real.h"
  14 #include "cl_F_tran.h"
  15 #include "cl_R.h"
  16
  17 #undef MAYBE_INLINE
  18 #define MAYBE_INLINE inline
  19 #include "cl_F_from_R_def.cc"
  20
  21 namespace cln {
  22
  23 // Hilfsfunktion für atanh und atan: u+iv := artanh(x+iy). Liefert cl_C_R(u,v).
  24
  25 const cl_C_R atanh (const cl_R& x, const cl_R& y)
  26 {
  27 // Methode:
  28 // Wert und Branch Cuts nach der Formel CLTL2, S. 315:
  29 //   artanh(z) = (log(1+z)-log(1-z)) / 2
  30 // Sei z=x+iy, Ergebnis u+iv.
  31 // Falls x=0 und y=0: u=0, v=0.
  32 // Falls x=0: u = 0, v = atan(X=1,Y=y).
  33 // Falls y=0:
  34 //   x rational -> x in Float umwandeln.
  35 //   |x|<1/2: u = atanh(x), v = 0.
  36 //   |x|>=1/2: (1+x)/(1-x) errechnen,
  37 //             =0 -> Error,
  38 //             >0 (also |x|<1) -> u = 1/2 log((1+x)/(1-x)), v = 0.
  39 //             <0 (also |x|>1) -> u = 1/2 log(-(1+x)/(1-x)),
  40 //                                v = (-pi/2 für x>1, pi/2 für x<-1).
  41 // Sonst:
  42 //   1+x und 1-x errechnen.
  43 //   x und y in Floats umwandeln.
  44 //   |4x| und 1+x^2+y^2 errechnen,
  45 //   |4x| < 1+x^2+y^2 -> u = 1/2 atanh(2x/(1+x^2+y^2)),
  46 //   |4x| >= 1+x^2+y^2 -> u = 1/4 ln ((1+x^2+y^2)+2x)/((1+x^2+y^2)-2x)
  47 //                        oder besser (an der Singularität: |x|-1,|y| klein):
  48 //                        u = 1/4 ln ((1+x)^2+y^2)/((1-x)^2+y^2).
  49 //   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,
  50 //                                          1 sonst)
  51 // Ergebnis ist reell nur, wenn z reell.
  52 // Real- und Imaginärteil des Ergebnisses sind Floats, außer wenn z reell oder
  53 // rein imaginär ist.
  54
  55         if (eq(x,0))
  56                 // x=0 -> u=0, v=atan(X=1,Y=y) (Fall y=0 ist inbegriffen)
  57                 return cl_C_R(0, atan(1,y));
  58         if (eq(y,0)) {
  59                 var cl_F xf = cl_float(x); // (float x)
  60                 var cl_F& x = xf;
  61                 // x Float
  62                 if (zerop(x))
  63                         // x=0.0 -> x als Ergebnis
  64                         return cl_C_R(x, 0);
  65                 if (float_exponent(x) < 0)
  66                         // Exponent e<0, also |x|<1/2
  67                         return cl_C_R(atanhx(x), 0);
  68                 // e>=0, also |x|>=1/2
  69                 var cl_F xx_den = 1 - x;
  70                 var cl_F xx = (1 + x) / xx_den; // (1+x)/(1-x)
  71                 var cl_R v;
  72                 if (!minusp(xx)) {
  73                         if (zerop(xx))
  74                                 { cl_error_division_by_0(); }
  75                         v = 0;
  76                 } else {
  77                         // (1+x)/(1-x) < 0 -> Betrag nehmen, Imaginärteil berechnen:
  78                         xx = - xx;
  79                         v = scale_float(pi(),-1); // (scale-float pi -1) = pi/2
  80                         if (minusp(xx_den))
  81                                 // 1-x<0 -> dann -pi/2
  82                                 v = -v;
  83                 }
  84                 // ln bilden, durch 2
  85                 return cl_C_R(scale_float(ln(xx),-1), v);
  86         }
  87         var cl_R _1_plus_x = 1+x;
  88         var cl_R _1_minus_x = 1-x;
  89         // x und y in Floats umwandeln: (Diese Fallunterscheidung ist
  90         // symmetrisch in x und y, auch wenn's nicht so aussieht.)
  91         var cl_F xf;
  92         var cl_F yf;
  93         if (rationalp(x)) {
  94                 DeclareType(cl_RA,x);
  95                 yf = cl_float(y);
  96                 xf = cl_float(x,yf);
  97         } else {
  98                 DeclareType(cl_F,x);
  99                 xf = x;
 100                 yf = cl_somefloat(y,xf);
 101         }
 102         var cl_F yf_2 = square(yf);
 103         var cl_F u;
 104         {
 105                 var cl_F temp1 = abs(scale_float(xf,2)); // |4x|
 106                 var cl_F temp2 = 1 + (square(xf) + yf_2); // 1+x^2+y^2
 107                 if (temp1 < temp2) // |4x| < 1+x^2+y^2 ?
 108                         // u = 1/2 atanh(2x/(1+x^2+y^2))
 109                         u = scale_float(atanhx(scale_float(xf,1)/temp2),-1);
 110                 else {
 111                         // u = 1/4 ln ((1+x)^2+y^2)/((1-x)^2+y^2)
 112                         var cl_F num = _1_plus_x*_1_plus_x + yf_2; // (1+x)^2+y^2, ein Float >=0
 113                         var cl_F den = _1_minus_x*_1_minus_x + yf_2; // (1-x)^2+y^2, ein Float >=0
 114                         if (zerop(den))
 115                                 { cl_error_division_by_0(); }
 116                         u = scale_float(ln(num/den),-2);
 117                 }
 118         }
 119         var cl_F v;
 120         {
 121                 var cl_F X = _1_plus_x*_1_minus_x-yf_2;
 122                 var cl_F Y = scale_float(yf,1);
 123                 v = atan(X,Y); // atan(X=(1-x)(1+x)-y^2,Y=2y), ein Float
 124                 if (minusp(X) && !minusp(x) && zerop(Y))
 125                         v = -v;
 126                 v = scale_float(v,-1); // 1/2 * atan(...) * +-1
 127         }
 128         return cl_C_R(u,v);
 129 }
 130
 131 }  // namespace cln