src/complex/elem/division/cl_C_SF_recip.cc

   1 // cl_C_recip().
   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/sfloat.h"
  13 #include "cl_SF.h"
  14
  15 namespace cln {
  16
  17 const cl_C_SF cl_C_recip (const cl_SF& a, const cl_SF& b)
  18 {
  19 //  a=0.0 -> liefere die Komponenten a=0.0 und -1/b.
  20 //  b=0.0 -> liefere die Komponenten 1/a und b=0.0.
  21 //  e:=max(exponent(a),exponent(b)).
  22 //  a':=a/2^e bzw. 0.0 bei Underflowmöglichkeit (beim Skalieren a':=a/2^e
  23 //      oder beim Quadrieren a'*a':  2*(e-exponent(a))>exp_mid-exp_low-1
  24 //      d.h. exponent(b)-exponent(a)>floor((exp_mid-exp_low-1)/2) ).
  25 //  b':=b/2^e bzw. 0.0 bei Underflowmöglichkeit (beim Skalieren b':=b/2^e
  26 //      oder beim Quadrieren b'*b':  2*(e-exponent(b))>exp_mid-exp_low-1
  27 //      d.h. exponent(a)-exponent(b)>floor((exp_mid-exp_low-1)/2) ).
  28 //  c':=a'*a'+b'*b',
  29 //  liefere die beiden Komponenten 2^(-e)*a'/c' und -2^(-e)*b'/c'.
  30         var sintL a_exp;
  31         var sintL b_exp;
  32         {
  33                 // Exponenten von a holen:
  34                 var uintL uexp = SF_uexp(a);
  35                 if (uexp == 0)
  36                         // a=0.0 -> liefere (complex a (- (/ b))) :
  37                         return cl_C_SF(a,-recip(b));
  38                 a_exp = (sintL)(uexp - SF_exp_mid);
  39         }
  40         {
  41                 // Exponenten von b holen:
  42                 var uintL uexp = SF_uexp(b);
  43                 if (uexp == 0)
  44                         // b=0.0 -> liefere (complex (/ a) b) :
  45                         return cl_C_SF(recip(a),b);
  46                 b_exp = (sintL)(uexp - SF_exp_mid);
  47         }
  48         // Nun a_exp = float_exponent(a), b_exp = float_exponent(b).
  49         var sintL e = (a_exp > b_exp ? a_exp : b_exp); // Maximum der Exponenten
  50         // a und b durch 2^e dividieren:
  51         var cl_SF na = (b_exp-a_exp > floor(SF_exp_mid-SF_exp_low-1,2) ? SF_0 : scale_float(a,-e));
  52         var cl_SF nb = (a_exp-b_exp > floor(SF_exp_mid-SF_exp_low-1,2) ? SF_0 : scale_float(b,-e));
  53         // c' := a'*a'+b'*b' berechnen:
  54         var cl_SF nc = square(na) + square(nb);
  55         // 2^(-e)*a'/c' + i * -2^(-e)*b'/c'
  56         return cl_C_SF(scale_float(na/nc,-e), scale_float(-(nb/nc),-e));
  57 }
  58
  59 }  // namespace cln