src/float/sfloat/elem/cl_SF_div.cc

   1 // binary operator /
   2
   3 // General includes.
   4 #include "cl_sysdep.h"
   5
   6 // Specification.
   7 #include "cln/sfloat.h"
   8
   9
  10 // Implementation.
  11
  12 #include "cl_SF.h"
  13 #include "cl_N.h"
  14 #include "cl_low.h"
  15
  16 namespace cln {
  17
  18 const cl_SF operator/ (const cl_SF& x1, const cl_SF& x2)
  19 {
  20 // Methode:
  21 // x2 = 0.0 -> Error
  22 // x1 = 0.0 -> Ergebnis 0.0
  23 // Sonst:
  24 // Ergebnis-Vorzeichen = xor der beiden Vorzeichen von x1 und x2
  25 // Ergebnis-Exponent = Differenz der beiden Exponenten von x1 und x2
  26 // Ergebnis-Mantisse = Mantisse mant1 / Mantisse mant2, gerundet.
  27 //   mant1/mant2 > 1/2, mant1/mant2 < 2;
  28 //   nach Rundung mant1/mant2 >=1/2, <=2*mant1<2.
  29 //   Bei mant1/mant2 >=1 brauche 16 Nachkommabits,
  30 //   bei mant1/mant2 <1 brauche 17 Nachkommabits.
  31 //   Fürs Runden: brauche ein Rundungsbit (Rest gibt an, ob exakt).
  32 //   Brauche daher insgesamt 18 Nachkommabits von mant1/mant2.
  33 //   Dividiere daher (als Unsigned Integers) 2^18*(2^17*mant1) durch (2^17*mant2).
  34 //   Falls der Quotient >=2^18 ist, runde die letzten zwei Bits weg und
  35 //     erhöhe den Exponenten um 1.
  36 //   Falls der Quotient <2^18 ist, runde das letzte Bit weg. Bei rounding
  37 //     overflow schiebe um ein weiteres Bit nach rechts, incr. Exponenten.
  38       // x1,x2 entpacken:
  39       var cl_signean sign1;
  40       var sintL exp1;
  41       var uintL mant1;
  42       var cl_signean sign2;
  43       var sintL exp2;
  44       var uintL mant2;
  45       SF_decode(x2, { throw division_by_0_exception(); }, sign2=,exp2=,mant2=);
  46       SF_decode(x1, { return x1; }, sign1=,exp1=,mant1=);
  47       exp1 = exp1 - exp2; // Differenz der Exponenten
  48       sign1 = sign1 ^ sign2; // Ergebnis-Vorzeichen
  49       // Dividiere 2^18*mant1 durch mant2 oder (äquivalent)
  50       // 2^i*2^18*mant1 durch 2^i*mant2 für irgendein i mit 0 <= i <= 32-17 :
  51       var uintL mant;
  52       var uintL rest;
  53       // wähle i = 32-(SF_mant_len+1), also i+(SF_mant_len+2) = 33.
  54       divu_6432_3232(mant1<<1,0, mant2<<(32-(SF_mant_len+1)), mant=,rest=);
  55       if (mant >= bit(SF_mant_len+2))
  56         // Quotient >=2^18 -> 2 Bits wegrunden
  57         { var uintL rounding_bits = mant & (bit(2)-1);
  58           exp1 += 1; // Exponenten incrementieren
  59           mant = mant >> 2;
  60           if ( (rounding_bits < bit(1)) // 00,01 werden abgerundet
  61                || ( (rounding_bits == bit(1)) // 10
  62                     && (rest == 0) // und genau halbzahlig
  63                     && ((mant & bit(0)) ==0) // -> round-to-even
  64              )    )
  65             // abrunden
  66             {}
  67             else
  68             // aufrunden
  69             { mant += 1; }
  70         }
  71         else
  72         // Quotient <2^18 -> 1 Bit wegrunden
  73         { var uintL rounding_bit = mant & bit(0);
  74           mant = mant >> 1;
  75           if ( (rounding_bit == 0) // 0 wird abgerundet
  76                || ( (rest == 0) // genau halbzahlig
  77                     && ((mant & bit(0)) ==0) // -> round-to-even
  78              )    )
  79             // abrunden
  80             {}
  81             else
  82             // aufrunden
  83             { mant += 1;
  84               if (mant >= bit(SF_mant_len+1)) // rounding overflow?
  85                 { mant = mant>>1; exp1 = exp1+1; }
  86         }   }
  87       return encode_SF(sign1,exp1,mant);
  88 }
  89
  90 }  // namespace cln