cl_R_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/real.h"
  13 #include "cl_R.h"
  14 #include "cln/rational.h"
  15 #include "cl_RA.h"
  16 #include "cl_F.h"
  17 #include "cl_SF.h"
  18 #include "cl_FF.h"
  19 #include "cl_DF.h"
  20 #include "cl_LF.h"
  21
  22 namespace cln {
  23
  24 const cl_R cl_hypot (const cl_R& a, const cl_R& b)
  25 {
  26 // Methode:
  27 // Falls a=0: (abs b).
  28 // Falls b=0: (abs a).
  29 // Falls a und b beide rational sind:
  30 //   c:=a*a+b*b, liefere (sqrt c).
  31 // Falls a oder b Floats sind:
  32 //   Falls einer von beiden rational ist, runde ihn zum selben Float-Typ
  33 //     wie der andere und führe das UP durch.
  34 //   Falls beide Floats sind, erweitere auf den genaueren, führe das UP
  35 //     durch und runde wieder auf den ungenaueren.
  36 //   Das Ergebnis ist ein Float >=0.
  37 // UP: [a,b Floats vom selben Typ]
  38 //  a=0.0 -> liefere abs(b).
  39 //  b=0.0 -> liefere abs(a).
  40 //  e:=max(exponent(a),exponent(b)).
  41 //  a':=a/2^e bzw. 0.0 bei Underflowmöglichkeit (beim Skalieren a':=a/2^e
  42 //      oder beim Quadrieren a'*a':  2*(e-exponent(a))>exp_mid-exp_low-1
  43 //      d.h. exponent(b)-exponent(a)>floor((exp_mid-exp_low-1)/2) ).
  44 //  b':=b/2^e bzw. 0.0 bei Underflowmöglichkeit (beim Skalieren b':=b/2^e
  45 //      oder beim Quadrieren b'*b':  2*(e-exponent(b))>exp_mid-exp_low-1
  46 //      d.h. exponent(a)-exponent(b)>floor((exp_mid-exp_low-1)/2) ).
  47 //  c':=a'*a'+b'*b', c':=sqrt(c'), liefere 2^e*c'.
  48
  49         if (rationalp(a)) {
  50                 DeclareType(cl_RA,a);
  51                 if (eq(a,0)) // a=0 -> (abs b)
  52                         return abs(b);
  53                 if (rationalp(b)) {
  54                         DeclareType(cl_RA,b);
  55                         // a,b beide rational
  56                         return sqrt(square(a)+square(b));
  57                 } else {
  58                         DeclareType(cl_F,b);
  59                         // a rational, b Float
  60                         floatcase(b
  61                         ,       return cl_hypot(cl_RA_to_SF(a),b);
  62                         ,       return cl_hypot(cl_RA_to_FF(a),b);
  63                         ,       return cl_hypot(cl_RA_to_DF(a),b);
  64                         ,       return cl_hypot(cl_RA_to_LF(a,TheLfloat(b)->len),b);
  65                         );
  66                 }
  67         } else {
  68                 DeclareType(cl_F,a);
  69                 if (rationalp(b)) {
  70                         DeclareType(cl_RA,b);
  71                         // a Float, b rational
  72                         if (eq(b,0)) // b=0 -> (abs a)
  73                                 return abs(a);
  74                         floatcase(a
  75                         ,       return cl_hypot(a,cl_RA_to_SF(b));
  76                         ,       return cl_hypot(a,cl_RA_to_FF(b));
  77                         ,       return cl_hypot(a,cl_RA_to_DF(b));
  78                         ,       return cl_hypot(a,cl_RA_to_LF(b,TheLfloat(a)->len));
  79                         );
  80                 } else {
  81                         DeclareType(cl_F,b);
  82                         // a,b Floats
  83                         #ifndef CL_LF_PEDANTIC
  84                         GEN_F_OP2(a,b, cl_hypot, 1, 1, return);
  85                         #else
  86                         GEN_F_OP2(a,b, cl_hypot, 1, 0, return);
  87                         #endif
  88                 }
  89         }
  90 }
  91
  92 }  // namespace cln