src/complex/algebraic/cl_LF_hypot.cc

   1 // cl_hypot().
   2
   3 // General includes.
   4 #include "cl_sysdep.h"
   5
   6 // Specification.
   7 #include "cl_C.h"
   8
   9
  10 // Implementation.
  11
  12 #include "cl_lfloat.h"
  13 #include "cl_LF.h"
  14 #include "cl_LF_impl.h"
  15
  16 #undef MAYBE_INLINE
  17 #define MAYBE_INLINE inline
  18 #include "cl_LF_minusp.cc"
  19
  20 ALL_cl_LF_OPERATIONS_SAME_PRECISION()
  21
  22 const cl_LF cl_hypot (const cl_LF& a, const cl_LF& b)
  23 {
  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'.
  35  {      Mutable(cl_LF,a);
  36         Mutable(cl_LF,b);
  37         {
  38                 var uintC a_len = TheLfloat(a)->len;
  39                 var uintC b_len = TheLfloat(b)->len;
  40                 if (!(a_len == b_len)) {
  41                         if (a_len < b_len)
  42                                 b = shorten(b,a_len);
  43                         else
  44                                 a = shorten(a,b_len);
  45                 }
  46         }
  47         var sintL a_exp;
  48         var sintL b_exp;
  49         {
  50                 // Exponenten von a holen:
  51                 var uintL uexp = TheLfloat(a)->expo;
  52                 if (uexp == 0)
  53                         // a=0.0 -> liefere (abs b) :
  54                         return (minusp(b) ? -b : b);
  55                 a_exp = (sintL)(uexp - LF_exp_mid);
  56         }
  57         {
  58                 // Exponenten von b holen:
  59                 var uintL uexp = TheLfloat(b)->expo;
  60                 if (uexp == 0)
  61                         // b=0.0 -> liefere (abs a) :
  62                         return (minusp(a) ? -a : a);
  63                 b_exp = (sintL)(uexp - LF_exp_mid);
  64         }
  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'
  73 }}