src/float/dfloat/elem/cl_DF_div.cc

   1 // binary operator /
   2
   3 // General includes.
   4 #include "cl_sysdep.h"
   5
   6 // Specification.
   7 #include "cln/dfloat.h"
   8
   9
  10 // Implementation.
  11
  12 #include "cl_DF.h"
  13 #include "cl_N.h"
  14 #include "cl_F.h"
  15 #include "cl_low.h"
  16 #include "cl_DS.h"
  17 #include "cl_ieee.h"
  18
  19 #undef MAYBE_INLINE
  20 #define MAYBE_INLINE inline
  21 #include "cl_DF_zerop.cc"
  22
  23 namespace cln {
  24
  25 NEED_IEEE_FLOATS()
  26
  27 const cl_DF operator/ (const cl_DF& x1, const cl_DF& x2)
  28 {
  29 // Methode:
  30 // x2 = 0.0 -> Error
  31 // x1 = 0.0 -> Ergebnis 0.0
  32 // Sonst:
  33 // Ergebnis-Vorzeichen = xor der beiden Vorzeichen von x1 und x2
  34 // Ergebnis-Exponent = Differenz der beiden Exponenten von x1 und x2
  35 // Ergebnis-Mantisse = Mantisse mant1 / Mantisse mant2, gerundet.
  36 //   mant1/mant2 > 1/2, mant1/mant2 < 2;
  37 //   nach Rundung mant1/mant2 >=1/2, <=2*mant1<2.
  38 //   Bei mant1/mant2 >=1 brauche 52 Nachkommabits,
  39 //   bei mant1/mant2 <1 brauche 53 Nachkommabits.
  40 //   Fürs Runden: brauche ein Rundungsbit (Rest gibt an, ob exakt).
  41 //   Brauche daher insgesamt 54 Nachkommabits von mant1/mant2.
  42 //   Dividiere daher (als Unsigned Integers) 2^54*(2^53*mant1) durch (2^53*mant2).
  43 //   Falls der Quotient >=2^54 ist, runde die letzten zwei Bits weg und
  44 //     erhöhe den Exponenten um 1.
  45 //   Falls der Quotient <2^54 ist, runde das letzte Bit weg. Bei rounding
  46 //     overflow schiebe um ein weiteres Bit nach rechts, incr. Exponenten.
  47 #if defined(FAST_DOUBLE) && !defined(__i386__)
  48       double_to_DF(DF_to_double(x1) / DF_to_double(x2), return ,
  49                    TRUE, TRUE, // Overflow und subnormale Zahl abfangen
  50                    !zerop(x1), // ein Ergebnis +/- 0.0
  51                                // ist genau dann in Wirklichkeit ein Underflow
  52                    zerop(x2), // Division durch Null abfangen
  53                    FALSE // kein NaN als Ergebnis möglich
  54                   );
  55 #else
  56       // x1,x2 entpacken:
  57       var cl_signean sign1;
  58       var sintL exp1;
  59       #if (intDsize<=32)
  60       var uintL manthi1;
  61       var uintL mantlo1;
  62       #endif
  63       var cl_signean sign2;
  64       var sintL exp2;
  65       #if (intDsize<=32)
  66       var uintL manthi2;
  67       var uintL mantlo2;
  68       #endif
  69       #if (cl_word_size==64)
  70       var uint64 mantx1;
  71       var uint64 mantx2;
  72       DF_decode(x2, { cl_error_division_by_0(); }, sign2=,exp2=,mantx2=);
  73       DF_decode(x1, { return x1; }, sign1=,exp1=,mantx1=);
  74       #else
  75       DF_decode2(x2, { cl_error_division_by_0(); }, sign2=,exp2=,manthi2=,mantlo2=);
  76       DF_decode2(x1, { return x1; }, sign1=,exp1=,manthi1=,mantlo1=);
  77       #endif
  78       exp1 = exp1 - exp2; // Differenz der Exponenten
  79       sign1 = sign1 ^ sign2; // Ergebnis-Vorzeichen
  80       // Dividiere 2^54*mant1 durch mant2 oder (äquivalent)
  81       // 2^i*2^54*mant1 durch 2^i*mant2 für irgendein i mit 0 <= i <= 64-53 :
  82       // wähle i = 64-(DF_mant_len+1), also i+(DF_mant_len+2) = 65.
  83       #if (cl_word_size==64)
  84       mantx1 = mantx1 << 1;
  85       mantx2 = mantx2 << (64-(DF_mant_len+1));
  86       #if (intDsize<=32)
  87       manthi1 = high32(mantx1); mantlo1 = low32(mantx1);
  88       manthi2 = high32(mantx2); mantlo2 = low32(mantx2);
  89       #endif
  90       #else
  91       manthi1 = (manthi1 << 1) | (mantlo1 >> 31); mantlo1 = mantlo1 << 1;
  92       manthi2 = (manthi2 << (64-(DF_mant_len+1))) | (mantlo2 >> ((DF_mant_len+1)-32)); mantlo2 = mantlo2 << (64-(DF_mant_len+1));
  93       #endif
  94       var uintD mant1 [128/intDsize];
  95       var uintD mant2 [64/intDsize];
  96       #if (intDsize==64)
  97       arrayLSref(mant1,128/intDsize,1) = mantx1;
  98       arrayLSref(mant1,128/intDsize,0) = 0;
  99       arrayLSref(mant2,64/intDsize,0) = mantx2;
 100       #elif (intDsize==32) || (intDsize==16) || (intDsize==8)
 101       set_32_Dptr(arrayMSDptr(mant1,128/intDsize),manthi1);
 102       set_32_Dptr(arrayMSDptr(mant1,128/intDsize) mspop 32/intDsize,mantlo1);
 103       set_32_Dptr(arrayMSDptr(mant1,128/intDsize) mspop 2*32/intDsize,0);
 104       set_32_Dptr(arrayMSDptr(mant1,128/intDsize) mspop 3*32/intDsize,0);
 105       set_32_Dptr(arrayMSDptr(mant2,64/intDsize),manthi2);
 106       set_32_Dptr(arrayMSDptr(mant2,64/intDsize) mspop 32/intDsize,mantlo2);
 107       #else
 108       {var uintD* ptr;
 109        ptr = arrayLSDptr(mant1,128/intDsize);
 110        doconsttimes(64/intDsize, { lsprefnext(ptr) = 0; } );
 111        doconsttimes(32/intDsize, { lsprefnext(ptr) = (uintD)mantlo1; mantlo1 = mantlo1>>intDsize; } );
 112        doconsttimes(32/intDsize, { lsprefnext(ptr) = (uintD)manthi1; manthi1 = manthi1>>intDsize; } );
 113       }
 114       {var uintD* ptr;
 115        ptr = arrayLSDptr(mant2,64/intDsize);
 116        doconsttimes(32/intDsize, { lsprefnext(ptr) = (uintD)mantlo2; mantlo2 = mantlo2>>intDsize; } );
 117        doconsttimes(32/intDsize, { lsprefnext(ptr) = (uintD)manthi2; manthi2 = manthi2>>intDsize; } );
 118       }
 119       #endif
 120       #if (cl_word_size==64)
 121       var uint64 mantx;
 122       #endif
 123       #if (intDsize<=32)
 124       var uintL manthi;
 125       var uintL mantlo;
 126       #endif
 127       {CL_ALLOCA_STACK;
 128        var DS q;
 129        var DS r;
 130        UDS_divide(arrayMSDptr(mant1,128/intDsize),128/intDsize,arrayLSDptr(mant1,128/intDsize),
 131                   arrayMSDptr(mant2,64/intDsize),64/intDsize,arrayLSDptr(mant2,64/intDsize),
 132                   &q, &r
 133                  );
 134        // Es ist 2^53 <= q < 2^55, also q.len = ceiling(54/intDsize)=ceiling(55/intDsize),
 135        // und r=0 genau dann, wenn r.len=0.
 136        ASSERT(q.len==ceiling(54,intDsize))
 137        {var uintD* ptr = q.MSDptr;
 138         #if (intDsize==64)
 139         mantx = mspref(ptr,0);
 140         #else // (intDsize<=32)
 141         manthi = get_max32_Dptr(23,ptr);
 142         mantlo = get_32_Dptr(ptr mspop ceiling(23,intDsize));
 143         #endif
 144        }
 145        // q = 2^32*manthi+mantlo.
 146        #if (cl_word_size==64)
 147        #if (intDsize<=32)
 148        mantx = (manthi<<32) | (uint64)mantlo;
 149        #endif
 150        if (mantx >= bit(DF_mant_len+2))
 151          // Quotient >=2^54 -> 2 Bits wegrunden
 152          { var uint64 rounding_bits = mantx & (bit(2)-1);
 153            exp1 += 1; // Exponenten incrementieren
 154            mantx = mantx >> 2;
 155            if ( (rounding_bits < bit(1)) // 00,01 werden abgerundet
 156                 || ( (rounding_bits == bit(1)) // 10
 157                      && (r.len == 0) // und genau halbzahlig
 158                      && ((mantx & bit(0)) ==0) // -> round-to-even
 159               )    )
 160              // abrunden
 161              {}
 162              else
 163              // aufrunden
 164              { mantx += 1; }
 165          }
 166          else
 167          // Quotient <2^54 -> 1 Bit wegrunden
 168          { var uint64 rounding_bit = mantx & bit(0);
 169            mantx = mantx >> 1;
 170            if ( (rounding_bit == 0) // 0 wird abgerundet
 171                 || ( (r.len == 0) // genau halbzahlig
 172                      && ((mantx & bit(0)) ==0) // -> round-to-even
 173               )    )
 174              // abrunden
 175              {}
 176              else
 177              // aufrunden
 178              { mantx += 1;
 179                if (mantx >= bit(DF_mant_len+1)) // rounding overflow?
 180                  { mantx = mantx>>1; exp1 = exp1+1; }
 181          }   }
 182        #else
 183        if (manthi >= bit(DF_mant_len-32+2))
 184          // Quotient >=2^54 -> 2 Bits wegrunden
 185          { var uintL rounding_bits = mantlo & (bit(2)-1);
 186            exp1 += 1; // Exponenten incrementieren
 187            mantlo = (mantlo >> 2) | (manthi << 30); manthi = manthi >> 2;
 188            if ( (rounding_bits < bit(1)) // 00,01 werden abgerundet
 189                 || ( (rounding_bits == bit(1)) // 10
 190                      && (r.len == 0) // und genau halbzahlig
 191                      && ((mantlo & bit(0)) ==0) // -> round-to-even
 192               )    )
 193              // abrunden
 194              {}
 195              else
 196              // aufrunden
 197              { mantlo += 1; if (mantlo==0) { manthi += 1; } }
 198          }
 199          else
 200          // Quotient <2^54 -> 1 Bit wegrunden
 201          { var uintL rounding_bit = mantlo & bit(0);
 202            mantlo = (mantlo >> 1) | (manthi << 31); manthi = manthi >> 1;
 203            if ( (rounding_bit == 0) // 0 wird abgerundet
 204                 || ( (r.len == 0) // genau halbzahlig
 205                      && ((mantlo & bit(0)) ==0) // -> round-to-even
 206               )    )
 207              // abrunden
 208              {}
 209              else
 210              // aufrunden
 211              { mantlo += 1;
 212                if (mantlo==0)
 213                  { manthi += 1;
 214                    if (manthi >= bit(DF_mant_len-32+1)) // rounding overflow?
 215                      { manthi = manthi>>1; exp1 = exp1+1; }
 216          }   }   }
 217        #endif
 218       }
 219       #if (cl_word_size==64)
 220       return encode_DF(sign1,exp1,mantx);
 221       #else
 222       return encode_DF(sign1,exp1,manthi,mantlo);
 223       #endif
 224 #endif
 225 }
 226
 227 }  // namespace cln