13 #include "cl_F_tran.h"
15 #include "cl_rational.h"
20 #define MAYBE_INLINE inline
21 #include "cl_F_from_R_def.cc"
23 // Hilfsfunktion für asinh und asin: u+iv := arsinh(x+iy). Liefert cl_C_R(u,v).
25 const cl_C_R asinh (const cl_R& x, const cl_R& y)
28 // Wert und Branch Cuts nach der Formel CLTL2, S. 313:
29 // arsinh(z) = log(z+sqrt(1+z^2))
30 // z=x+iy, Ergebnis u+iv.
31 // Falls x=0 und y=0: u=0, v=0.
32 // Falls x=0: arsinh(iy) = i arcsin(y).
34 // Bei y=1: u = 0, v = pi/2.
35 // Bei y=1/2: u = 0, v = pi/6.
36 // Bei y=0: u = 0, v = 0.
37 // Bei y=-1/2: u = 0, v = -pi/6.
38 // Bei y=-1: u = 0, v = -pi/2.
39 // Sonst y in Float umwandeln.
40 // e := Exponent aus (decode-float y), d := (float-digits y)
41 // Bei y=0.0 oder e<=-d/2 liefere u = 0, v = y
42 // (denn bei e<=-d/2 ist y^2/3 < y^2/2 < 2^(-d)/2 = 2^(-d-1), also
43 // 1 <= asin(y)/y < 1+y^2/3 < 1+2^(-d-1) < 1+2^(-d),
44 // also ist asin(y)/y, auf d Bits gerundet, gleich 1.0).
46 // Bei y>1 liefere u = ln(y+sqrt(y^2-1)), v = pi/2.
47 // Bei y<-1 liefere u = -ln(|y|+sqrt(|y|^2-1)), v = -pi/2.
48 // Bei |y|<=1 liefere u = 0, v = atan(X=sqrt(1-y^2),Y=y).
50 // x rational -> x in Float umwandeln.
51 // |x|<1/2: u = atanh(x/sqrt(1+x^2)),
52 // x>=1/2: u = ln(x+sqrt(1+x^2)),
53 // x<=-1/2: u = -ln(-x+sqrt(1+x^2)).
56 // z in Bild(sqrt) -> log(sqrt(1+z^2)+z) = (!) = 2 artanh(z/(1+sqrt(1+z^2))).
57 // z nicht in Bild(sqrt) ->
58 // arsinh(z) = -arsinh(-z).
59 // (Denn arsinh(z)+arsinh(-z) == log((z+sqrt(1+z^2))(-z+sqrt(1+z^2)))
60 // = log((1+z^2)-z^2) = log(1) = 0 mod 2 pi i, und links ist
61 // der Imaginärteil betragsmäßig <=pi.)
62 // Also arsinh(z) = -arsinh(-z) = - 2 artanh(-z/(1+sqrt(1+z^2)))
63 // = (wegen -artanh(-w) = artanh(w)) = 2 artanh(z/(1+sqrt(1+z^2))).
64 // Real- und Imaginärteil des Ergebnisses sind Floats, außer wenn z reell oder
67 // Um für zwei Zahlen u,v mit u^2-v^2=1 und u,v beide in Bild(sqrt)
68 // (d.h. Realteil>0.0 oder Realteil=0.0 und Imaginärteil>=0.0)
69 // log(u+v) zu berechnen:
70 // log(u+v) = 2 artanh(v/(u+1)) (!)
71 // (Beweis: 2 artanh(v/(u+1)) = log(1+(v/(u+1))) - log(1-(v/(u+1)))
72 // = log((1+u+v)/(u+1)) - log((1+u-v)/(u+1)) == log((1+u+v)/(1+u-v))
73 // = log(u+v) mod 2 pi i, und beider Imaginärteil ist > -pi und <= pi.)
81 if (eq(y,0)) // x=0, y=0 -> u=0, v=0
86 if (eq(y,1)) // x=0, y=1 -> v = pi/2
87 return cl_C_R(0,scale_float(cl_pi(),-1));
88 if (eq(y,-1)) // x=0, y=-1 -> v = -pi/2
89 return cl_C_R(0,-scale_float(cl_pi(),-1));
90 yf = cl_float(y); // y in Float umwandeln
94 if (eq(denominator(y),2)) { // Nenner = 2 ?
95 if (eq(numerator(y),1)) // x=0, y=1/2 -> v = pi/6
96 return cl_C_R(0,cl_pi()/6);
97 if (eq(numerator(y),-1)) // x=0, y=-1/2 -> v = -pi/6
98 return cl_C_R(0,-(cl_pi()/6));
100 yf = cl_float(y); // y in Float umwandeln
108 if (zerop(y)) // y=0.0 -> arcsin(y) = y als Ergebnis
110 if (float_exponent(y) <= (-(sintL)float_digits(y))>>1)
111 // e <= -d/2 <==> e <= -ceiling(d/2)
113 var cl_F temp = 1-square(y);
115 // 1-y*y>=0, also |y|<=1
116 // v = atan(X=sqrt(1-y*y),Y=y)
117 return cl_C_R(0,atan(sqrt(temp),y));
119 // 1-y*y<0, also |y|>1
120 temp = sqrt(-temp); // sqrt(y*y-1)
125 // temp = sqrt(y^2-1)+|y|, ein Float >1
126 var cl_F u = ln(temp); // ln(|y|+sqrt(y^2-1)), ein Float >0
127 var cl_F v = scale_float(cl_pi(),-1); // (scale-float pi -1) = pi/2
129 return cl_C_R(u,v); // y>1 -> v = pi/2
131 return cl_C_R(-u,-v); // y<-1 -> v = -pi/2, u = -ln(...)
136 var cl_F xf = cl_float(x); // x in Float umwandeln
140 return cl_C_R(x,0); // x=0.0 -> u=x, v=0.
141 var cl_F temp = sqrt(1+square(x)); // sqrt(1+x^2)
142 if (float_exponent(x) < 0) // Exponent e (von x/=0) <0 ?
144 return cl_C_R(atanhx(x/temp),0);
149 return cl_C_R(ln(temp+x),0); // u = ln(x+sqrt(1+x^2))
152 return cl_C_R(-ln(temp-x),0); // u = -ln(-x+sqrt(1+x^2))
154 var cl_N z = complex_C(x,y); // z=x+iy
155 var cl_N w = z/(1+sqrt(1+square(z))); // z/(1+sqrt(1+z^2))
156 // Da z=x+iy weder reell noch rein imaginär ist, ist auch
157 // w := z/(1+sqrt(1+z^2)) weder reell noch rein imaginär.
158 // (Beweis: Sollte sqrt(1+z^2) rationalen Real- und Imaginärteil haben,
159 // so auch z, also auch w, und die Formel z = 2w/(1-w^2) zeigt, daß dann
160 // z reell oder rein imaginär sein müßte. Also hat sqrt(1+z^2) ein
161 // Float als Real- oder Imaginärteil, das Betragsquadrat des Nenners
162 // ist also ein Float, und da Real- und Imaginärteil von z /=0 sind,
163 // sind Real- und Imaginärteil von w Floats.)
164 // Daher hat dann atanh(...) Floats als Realteil u und Imaginärteil v.
165 { DeclareType(cl_C,w);
166 cl_C_R u_v = atanh(realpart(w),imagpart(w));
167 var cl_R& u = u_v.realpart;
168 var cl_R& v = u_v.imagpart;
169 { DeclareType(cl_F,u);
171 return cl_C_R(scale_float(u,1),scale_float(v,1)); // u:=2*u, v:=2*v