4 #include "base/cl_sysdep.h"
7 #include "complex/cl_C.h"
12 #include "cln/ffloat.h"
13 #include "float/ffloat/cl_FF.h"
17 const cl_C_FF cl_C_recip (const cl_FF& a, const cl_FF& b)
19 // a=0.0 -> liefere die Komponenten a=0.0 und -1/b.
20 // b=0.0 -> liefere die Komponenten 1/a und b=0.0.
21 // e:=max(exponent(a),exponent(b)).
22 // a':=a/2^e bzw. 0.0 bei Underflowmöglichkeit (beim Skalieren a':=a/2^e
23 // oder beim Quadrieren a'*a': 2*(e-exponent(a))>exp_mid-exp_low-1
24 // d.h. exponent(b)-exponent(a)>floor((exp_mid-exp_low-1)/2) ).
25 // b':=b/2^e bzw. 0.0 bei Underflowmöglichkeit (beim Skalieren b':=b/2^e
26 // oder beim Quadrieren b'*b': 2*(e-exponent(b))>exp_mid-exp_low-1
27 // d.h. exponent(a)-exponent(b)>floor((exp_mid-exp_low-1)/2) ).
29 // liefere die beiden Komponenten 2^(-e)*a'/c' und -2^(-e)*b'/c'.
33 // Exponenten von a holen:
34 var uintL uexp = FF_uexp(cl_ffloat_value(a));
36 // a=0.0 -> liefere (complex a (- (/ b))) :
37 return cl_C_FF(a,-recip(b));
38 a_exp = (sintL)(uexp - FF_exp_mid);
41 // Exponenten von b holen:
42 var uintL uexp = FF_uexp(cl_ffloat_value(b));
44 // b=0.0 -> liefere (complex (/ a) b) :
45 return cl_C_FF(recip(a),b);
46 b_exp = (sintL)(uexp - FF_exp_mid);
48 // Nun a_exp = float_exponent(a), b_exp = float_exponent(b).
49 var sintL e = (a_exp > b_exp ? a_exp : b_exp); // Maximum der Exponenten
50 // a und b durch 2^e dividieren:
51 var cl_FF na = (b_exp-a_exp > floor(FF_exp_mid-FF_exp_low-1,2) ? cl_FF_0 : scale_float(a,-e));
52 var cl_FF nb = (a_exp-b_exp > floor(FF_exp_mid-FF_exp_low-1,2) ? cl_FF_0 : scale_float(b,-e));
53 // c' := a'*a'+b'*b' berechnen:
54 var cl_FF nc = square(na) + square(nb);
55 // 2^(-e)*a'/c' + i * -2^(-e)*b'/c'
56 return cl_C_FF(scale_float(na/nc,-e), scale_float(-(nb/nc),-e));