12 #include "cl_ffloat.h"
15 const cl_C_FF cl_C_recip (const cl_FF& a, const cl_FF& b)
17 // a=0.0 -> liefere die Komponenten a=0.0 und -1/b.
18 // b=0.0 -> liefere die Komponenten 1/a und b=0.0.
19 // e:=max(exponent(a),exponent(b)).
20 // a':=a/2^e bzw. 0.0 bei Underflowmöglichkeit (beim Skalieren a':=a/2^e
21 // oder beim Quadrieren a'*a': 2*(e-exponent(a))>exp_mid-exp_low-1
22 // d.h. exponent(b)-exponent(a)>floor((exp_mid-exp_low-1)/2) ).
23 // b':=b/2^e bzw. 0.0 bei Underflowmöglichkeit (beim Skalieren b':=b/2^e
24 // oder beim Quadrieren b'*b': 2*(e-exponent(b))>exp_mid-exp_low-1
25 // d.h. exponent(a)-exponent(b)>floor((exp_mid-exp_low-1)/2) ).
27 // liefere die beiden Komponenten 2^(-e)*a'/c' und -2^(-e)*b'/c'.
31 // Exponenten von a holen:
32 var uintL uexp = FF_uexp(cl_ffloat_value(a));
34 // a=0.0 -> liefere (complex a (- (/ b))) :
35 return cl_C_FF(a,-recip(b));
36 a_exp = (sintL)(uexp - FF_exp_mid);
39 // Exponenten von b holen:
40 var uintL uexp = FF_uexp(cl_ffloat_value(b));
42 // b=0.0 -> liefere (complex (/ a) b) :
43 return cl_C_FF(recip(a),b);
44 b_exp = (sintL)(uexp - FF_exp_mid);
46 // Nun a_exp = float_exponent(a), b_exp = float_exponent(b).
47 var sintL e = (a_exp > b_exp ? a_exp : b_exp); // Maximum der Exponenten
48 // a und b durch 2^e dividieren:
49 var cl_FF na = (b_exp-a_exp > floor(FF_exp_mid-FF_exp_low-1,2) ? cl_FF_0 : scale_float(a,-e));
50 var cl_FF nb = (a_exp-b_exp > floor(FF_exp_mid-FF_exp_low-1,2) ? cl_FF_0 : scale_float(b,-e));
51 // c' := a'*a'+b'*b' berechnen:
52 var cl_FF nc = square(na) + square(nb);
53 // 2^(-e)*a'/c' + i * -2^(-e)*b'/c'
54 return cl_C_FF(scale_float(na/nc,-e), scale_float(-(nb/nc),-e));