12 #include "cl_ffloat.h"
16 #define MAYBE_INLINE inline
17 #include "cl_FF_minusp.cc"
19 const cl_FF cl_hypot (const cl_FF& a, const cl_FF& b)
21 // a=0.0 -> liefere abs(b).
22 // b=0.0 -> liefere abs(a).
23 // e:=max(exponent(a),exponent(b)).
24 // a':=a/2^e bzw. 0.0 bei Underflowmöglichkeit (beim Skalieren a':=a/2^e
25 // oder beim Quadrieren a'*a': 2*(e-exponent(a))>exp_mid-exp_low-1
26 // d.h. exponent(b)-exponent(a)>floor((exp_mid-exp_low-1)/2) ).
27 // b':=b/2^e bzw. 0.0 bei Underflowmöglichkeit (beim Skalieren b':=b/2^e
28 // oder beim Quadrieren b'*b': 2*(e-exponent(b))>exp_mid-exp_low-1
29 // d.h. exponent(a)-exponent(b)>floor((exp_mid-exp_low-1)/2) ).
30 // c':=a'*a'+b'*b', c':=sqrt(c'), liefere 2^e*c'.
34 // Exponenten von a holen:
35 var uintL uexp = FF_uexp(cl_ffloat_value(a));
37 // a=0.0 -> liefere (abs b) :
38 return (minusp(b) ? -b : b);
39 a_exp = (sintL)(uexp - FF_exp_mid);
42 // Exponenten von b holen:
43 var uintL uexp = FF_uexp(cl_ffloat_value(b));
45 // b=0.0 -> liefere (abs a) :
46 return (minusp(a) ? -a : a);
47 b_exp = (sintL)(uexp - FF_exp_mid);
49 // Nun a_exp = float_exponent(a), b_exp = float_exponent(b).
50 var sintL e = (a_exp > b_exp ? a_exp : b_exp); // Maximum der Exponenten
51 // a und b durch 2^e dividieren:
52 var cl_FF na = (b_exp-a_exp > floor(FF_exp_mid-FF_exp_low-1,2) ? cl_FF_0 : scale_float(a,-e));
53 var cl_FF nb = (a_exp-b_exp > floor(FF_exp_mid-FF_exp_low-1,2) ? cl_FF_0 : scale_float(b,-e));
54 // c' := a'*a'+b'*b' berechnen:
55 var cl_FF nc = square(na) + square(nb);
56 return scale_float(sqrt(nc),e); // c' := sqrt(c'), 2^e*c'