src/complex/algebraic/cl_SF_hypot.cc

   1 // cl_hypot().
   2
   3 // General includes.
   4 #include "base/cl_sysdep.h"
   5
   6 // Specification.
   7 #include "complex/cl_C.h"
   8
   9
  10 // Implementation.
  11
  12 #include "cln/sfloat.h"
  13 #include "float/sfloat/cl_SF.h"
  14
  15 /* For inline copy of minusp */
  16 #include "base/cl_inline.h"
  17 #include "float/sfloat/elem/cl_SF_minusp.cc"
  18
  19 namespace cln {
  20
  21 const cl_SF cl_hypot (const cl_SF& a, const cl_SF& b)
  22 {
  23 //  a=0.0 -> liefere abs(b).
  24 //  b=0.0 -> liefere abs(a).
  25 //  e:=max(exponent(a),exponent(b)).
  26 //  a':=a/2^e bzw. 0.0 bei Underflowmöglichkeit (beim Skalieren a':=a/2^e
  27 //      oder beim Quadrieren a'*a':  2*(e-exponent(a))>exp_mid-exp_low-1
  28 //      d.h. exponent(b)-exponent(a)>floor((exp_mid-exp_low-1)/2) ).
  29 //  b':=b/2^e bzw. 0.0 bei Underflowmöglichkeit (beim Skalieren b':=b/2^e
  30 //      oder beim Quadrieren b'*b':  2*(e-exponent(b))>exp_mid-exp_low-1
  31 //      d.h. exponent(a)-exponent(b)>floor((exp_mid-exp_low-1)/2) ).
  32 //  c':=a'*a'+b'*b', c':=sqrt(c'), liefere 2^e*c'.
  33         var sintL a_exp;
  34         var sintL b_exp;
  35         {
  36                 // Exponenten von a holen:
  37                 var uintL uexp = SF_uexp(a);
  38                 if (uexp == 0)
  39                         // a=0.0 -> liefere (abs b) :
  40                         return (minusp_inline(b) ? -b : b);
  41                 a_exp = (sintL)(uexp - SF_exp_mid);
  42         }
  43         {
  44                 // Exponenten von b holen:
  45                 var uintL uexp = SF_uexp(b);
  46                 if (uexp == 0)
  47                         // b=0.0 -> liefere (abs a) :
  48                         return (minusp_inline(a) ? -a : a);
  49                 b_exp = (sintL)(uexp - SF_exp_mid);
  50         }
  51         // Nun a_exp = float_exponent(a), b_exp = float_exponent(b).
  52         var sintL e = (a_exp > b_exp ? a_exp : b_exp); // Maximum der Exponenten
  53         // a und b durch 2^e dividieren:
  54         var cl_SF na = (b_exp-a_exp > floor(SF_exp_mid-SF_exp_low-1,2) ? SF_0 : scale_float(a,-e));
  55         var cl_SF nb = (a_exp-b_exp > floor(SF_exp_mid-SF_exp_low-1,2) ? SF_0 : scale_float(b,-e));
  56         // c' := a'*a'+b'*b' berechnen:
  57         var cl_SF nc = square(na) + square(nb);
  58         return scale_float(sqrt(nc),e); // c' := sqrt(c'), 2^e*c'
  59 }
  60
  61 }  // namespace cln