src/float/ffloat/elem/cl_FF_div.cc

   1 // binary operator /
   2
   3 // General includes.
   4 #include "cl_sysdep.h"
   5
   6 // Specification.
   7 #include "cln/ffloat.h"
   8
   9
  10 // Implementation.
  11
  12 #include "cl_FF.h"
  13 #include "cl_N.h"
  14 #include "cl_F.h"
  15 #include "cl_low.h"
  16 #include "cl_ieee.h"
  17
  18 #undef MAYBE_INLINE
  19 #define MAYBE_INLINE inline
  20 #include "cl_FF_zerop.cc"
  21
  22 namespace cln {
  23
  24 NEED_IEEE_FLOATS()
  25
  26 const cl_FF operator/ (const cl_FF& x1, const cl_FF& x2)
  27 {
  28 // Methode:
  29 // x2 = 0.0 -> Error
  30 // x1 = 0.0 -> Ergebnis 0.0
  31 // Sonst:
  32 // Ergebnis-Vorzeichen = xor der beiden Vorzeichen von x1 und x2
  33 // Ergebnis-Exponent = Differenz der beiden Exponenten von x1 und x2
  34 // Ergebnis-Mantisse = Mantisse mant1 / Mantisse mant2, gerundet.
  35 //   mant1/mant2 > 1/2, mant1/mant2 < 2;
  36 //   nach Rundung mant1/mant2 >=1/2, <=2*mant1<2.
  37 //   Bei mant1/mant2 >=1 brauche 23 Nachkommabits,
  38 //   bei mant1/mant2 <1 brauche 24 Nachkommabits.
  39 //   Fürs Runden: brauche ein Rundungsbit (Rest gibt an, ob exakt).
  40 //   Brauche daher insgesamt 25 Nachkommabits von mant1/mant2.
  41 //   Dividiere daher (als Unsigned Integers) 2^25*(2^24*mant1) durch (2^24*mant2).
  42 //   Falls der Quotient >=2^25 ist, runde die letzten zwei Bits weg und
  43 //     erhöhe den Exponenten um 1.
  44 //   Falls der Quotient <2^25 ist, runde das letzte Bit weg. Bei rounding
  45 //     overflow schiebe um ein weiteres Bit nach rechts, incr. Exponenten.
  46   #if defined(FAST_FLOAT) && !defined(__i386__)
  47       float_to_FF(FF_to_float(x1) / FF_to_float(x2), return ,
  48                   TRUE, TRUE, // Overflow und subnormale Zahl abfangen
  49                   !zerop(x1), // ein Ergebnis +/- 0.0
  50                               // ist genau dann in Wirklichkeit ein Underflow
  51                   zerop(x2), // Division durch Null abfangen
  52                   FALSE // kein NaN als Ergebnis möglich
  53                  );
  54   #else
  55       // x1,x2 entpacken:
  56       var cl_signean sign1;
  57       var sintL exp1;
  58       var uintL mant1;
  59       var cl_signean sign2;
  60       var sintL exp2;
  61       var uintL mant2;
  62       FF_decode(x2, { throw division_by_0_exception(); }, sign2=,exp2=,mant2=);
  63       FF_decode(x1, { return x1; }, sign1=,exp1=,mant1=);
  64       exp1 = exp1 - exp2; // Differenz der Exponenten
  65       sign1 = sign1 ^ sign2; // Ergebnis-Vorzeichen
  66       // Dividiere 2^25*mant1 durch mant2 oder (äquivalent)
  67       // 2^i*2^25*mant1 durch 2^i*mant2 für irgendein i mit 0 <= i <= 32-24 :
  68       var uintL mant;
  69       var uintL rest;
  70       // wähle i = 32-(FF_mant_len+1), also i+(FF_mant_len+2) = 33.
  71       divu_6432_3232(mant1<<1,0, mant2<<(32-(FF_mant_len+1)), mant=,rest=);
  72       if (mant >= bit(FF_mant_len+2))
  73         // Quotient >=2^25 -> 2 Bits wegrunden
  74         { var uintL rounding_bits = mant & (bit(2)-1);
  75           exp1 += 1; // Exponenten incrementieren
  76           mant = mant >> 2;
  77           if ( (rounding_bits < bit(1)) // 00,01 werden abgerundet
  78                || ( (rounding_bits == bit(1)) // 10
  79                     && (rest == 0) // und genau halbzahlig
  80                     && ((mant & bit(0)) ==0) // -> round-to-even
  81              )    )
  82             // abrunden
  83             {}
  84             else
  85             // aufrunden
  86             { mant += 1; }
  87         }
  88         else
  89         // Quotient <2^25 -> 1 Bit wegrunden
  90         { var uintL rounding_bit = mant & bit(0);
  91           mant = mant >> 1;
  92           if ( (rounding_bit == 0) // 0 wird abgerundet
  93                || ( (rest == 0) // genau halbzahlig
  94                     && ((mant & bit(0)) ==0) // -> round-to-even
  95              )    )
  96             // abrunden
  97             {}
  98             else
  99             // aufrunden
 100             { mant += 1;
 101               if (mant >= bit(FF_mant_len+1)) // rounding overflow?
 102                 { mant = mant>>1; exp1 = exp1+1; }
 103         }   }
 104       return encode_FF(sign1,exp1,mant);
 105   #endif
 106 }
 107
 108 }  // namespace cln