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