14 #include "cl_F_tran.h"
18 #define MAYBE_INLINE inline
19 #include "cl_F_from_R_def.cc"
23 // Hilfsfunktion für atanh und atan: u+iv := artanh(x+iy). Liefert cl_C_R(u,v).
25 const cl_C_R atanh (const cl_R& x, const cl_R& y)
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).
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,
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).
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,
51 // Ergebnis ist reell nur, wenn z reell.
52 // Real- und Imaginärteil des Ergebnisses sind Floats, außer wenn z reell oder
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));
59 var cl_F xf = cl_float(x); // (float x)
63 // x=0.0 -> x als Ergebnis
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)
74 { throw division_by_0_exception(); }
77 // (1+x)/(1-x) < 0 -> Betrag nehmen, Imaginärteil berechnen:
79 v = scale_float(pi(),-1); // (scale-float pi -1) = pi/2
81 // 1-x<0 -> dann -pi/2
85 return cl_C_R(scale_float(ln(xx),-1), v);
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.)
100 yf = cl_somefloat(y,xf);
102 var cl_F yf_2 = square(yf);
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);
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
115 { throw division_by_0_exception(); }
116 u = scale_float(ln(num/den),-2);
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))
126 v = scale_float(v,-1); // 1/2 * atan(...) * +-1