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