12 #include "cl_lfloat.h"
14 #include "cl_LF_impl.h"
17 #define MAYBE_INLINE inline
18 #include "cl_LF_minusp.cc"
20 ALL_cl_LF_OPERATIONS_SAME_PRECISION()
22 const cl_LF cl_hypot (const cl_LF& a, const cl_LF& b)
24 // Zuerst a und b auf gleiche Länge bringen: den längeren runden.
25 // a=0.0 -> liefere abs(b).
26 // b=0.0 -> liefere abs(a).
27 // e:=max(exponent(a),exponent(b)).
28 // a':=a/2^e bzw. 0.0 bei Underflowmöglichkeit (beim Skalieren a':=a/2^e
29 // oder beim Quadrieren a'*a': 2*(e-exponent(a))>exp_mid-exp_low-1
30 // d.h. exponent(b)-exponent(a)>floor((exp_mid-exp_low-1)/2) ).
31 // b':=b/2^e bzw. 0.0 bei Underflowmöglichkeit (beim Skalieren b':=b/2^e
32 // oder beim Quadrieren b'*b': 2*(e-exponent(b))>exp_mid-exp_low-1
33 // d.h. exponent(a)-exponent(b)>floor((exp_mid-exp_low-1)/2) ).
34 // c':=a'*a'+b'*b', c':=sqrt(c'), liefere 2^e*c'.
38 var uintC a_len = TheLfloat(a)->len;
39 var uintC b_len = TheLfloat(b)->len;
40 if (!(a_len == b_len)) {
50 // Exponenten von a holen:
51 var uintL uexp = TheLfloat(a)->expo;
53 // a=0.0 -> liefere (abs b) :
54 return (minusp(b) ? -b : b);
55 a_exp = (sintL)(uexp - LF_exp_mid);
58 // Exponenten von b holen:
59 var uintL uexp = TheLfloat(b)->expo;
61 // b=0.0 -> liefere (abs a) :
62 return (minusp(a) ? -a : a);
63 b_exp = (sintL)(uexp - LF_exp_mid);
65 // Nun a_exp = float_exponent(a), b_exp = float_exponent(b).
66 var sintL e = (a_exp > b_exp ? a_exp : b_exp); // Maximum der Exponenten
67 // a und b durch 2^e dividieren:
68 var cl_LF na = ((b_exp > a_exp) && ((uintL)(b_exp-a_exp) > (uintL)floor(LF_exp_mid-LF_exp_low-1,2)) ? encode_LF0(TheLfloat(a)->len) : scale_float(a,-e));
69 var cl_LF nb = ((a_exp > b_exp) && ((uintL)(a_exp-b_exp) > (uintL)floor(LF_exp_mid-LF_exp_low-1,2)) ? encode_LF0(TheLfloat(b)->len) : scale_float(b,-e));
70 // c' := a'*a'+b'*b' berechnen:
71 var cl_LF nc = square(na) + square(nb);
72 return scale_float(sqrt(nc),e); // c' := sqrt(c'), 2^e*c'