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