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 "cln/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 namespace cln {
  21
  22 ALL_cl_LF_OPERATIONS_SAME_PRECISION()
  23
  24 const cl_LF cl_hypot (const cl_LF& a, const cl_LF& b)
  25 {
  26 //  Zuerst a und b auf gleiche Länge bringen: den längeren runden.
  27 //  a=0.0 -> liefere abs(b).
  28 //  b=0.0 -> liefere abs(a).
  29 //  e:=max(exponent(a),exponent(b)).
  30 //  a':=a/2^e bzw. 0.0 bei Underflowmöglichkeit (beim Skalieren a':=a/2^e
  31 //      oder beim Quadrieren a'*a':  2*(e-exponent(a))>exp_mid-exp_low-1
  32 //      d.h. exponent(b)-exponent(a)>floor((exp_mid-exp_low-1)/2) ).
  33 //  b':=b/2^e bzw. 0.0 bei Underflowmöglichkeit (beim Skalieren b':=b/2^e
  34 //      oder beim Quadrieren b'*b':  2*(e-exponent(b))>exp_mid-exp_low-1
  35 //      d.h. exponent(a)-exponent(b)>floor((exp_mid-exp_low-1)/2) ).
  36 //  c':=a'*a'+b'*b', c':=sqrt(c'), liefere 2^e*c'.
  37  {      Mutable(cl_LF,a);
  38         Mutable(cl_LF,b);
  39         {
  40                 var uintC a_len = TheLfloat(a)->len;
  41                 var uintC b_len = TheLfloat(b)->len;
  42                 if (!(a_len == b_len)) {
  43                         if (a_len < b_len)
  44                                 b = shorten(b,a_len);
  45                         else
  46                                 a = shorten(a,b_len);
  47                 }
  48         }
  49         var sintL a_exp;
  50         var sintL b_exp;
  51         {
  52                 // Exponenten von a holen:
  53                 var uintL uexp = TheLfloat(a)->expo;
  54                 if (uexp == 0)
  55                         // a=0.0 -> liefere (abs b) :
  56                         return (minusp(b) ? -b : b);
  57                 a_exp = (sintL)(uexp - LF_exp_mid);
  58         }
  59         {
  60                 // Exponenten von b holen:
  61                 var uintL uexp = TheLfloat(b)->expo;
  62                 if (uexp == 0)
  63                         // b=0.0 -> liefere (abs a) :
  64                         return (minusp(a) ? -a : a);
  65                 b_exp = (sintL)(uexp - LF_exp_mid);
  66         }
  67         // Nun a_exp = float_exponent(a), b_exp = float_exponent(b).
  68         var sintL e = (a_exp > b_exp ? a_exp : b_exp); // Maximum der Exponenten
  69         // a und b durch 2^e dividieren:
  70         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));
  71         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));
  72         // c' := a'*a'+b'*b' berechnen:
  73         var cl_LF nc = square(na) + square(nb);
  74         return scale_float(sqrt(nc),e); // c' := sqrt(c'), 2^e*c'
  75 }}
  76
  77 }  // namespace cln