Changeset 3177

Timestamp:

02/13/08 06:03:42 (10 months ago)

Author:

Don Clugston

Message:

Refactoring to support both big-endian and little-endian quadruple reals. All functions in math.IEEE are now supported, although none have actually been tested...

Files:

• trunk/tango/math/IEEE.d (modified) (37 diffs)

trunk/tango/math/IEEE.d

@@ -103 +103 @@
-}
-
+private:
-/* Most of the functions depend on the format of the largest IEEE floating-point type.
- * These code will differ depending on whether 'real' is 64, 80, or 128 bits,
- * and whether it is a big-endian or little-endian architecture.
-
-
-
-
-
+
+
+
+
+
+
+
- * There is also an unsupported ABI which does not follow IEEE; several of its functions
- *  will generate run-time errors if used.
-(double-double, as used by GDC <= 0.23
+'doubledouble' (used by GDC <= 0.23 for PowerPC
- */
-
-version(LittleEndian) {
-    static assert(real.mant_dig == 53 || real.mant_dig==64 || real.mant_dig == 113,
-
+-bit
-} else {
-    static assert(real.mant_dig == 53 || real.mant_dig==106 || real.mant_dig == 113,
-
-
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
-
-/** IEEE exception status flags
@@ -367 +415 @@
-    uint ex;
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
+
+
-    if (ex) { // If exponent is non-zero
-        if (ex == EXPMASK) {   // infinity or NaN
@@ -387 +422 @@
-                *vl |= 0xC000000000000000;  // convert $(NAN)S to $(NAN)Q
-                exp = int.min;
-
+_SHORT
-                exp = int.min;
-            } else {   // positive infinity
@@ -394 +429 @@
-        } else {
-            exp = ex - 0x3FFE;
-] = cast(ushort)((0x8000 & vu[EXPONENTPOS
+_SHORT] = cast(ushort)((0x8000 & vu[EXPONENTPOS_SHORT
-        }
-    } else if (!*vl) {
@@ -407 +442 @@
-        } while (*vl > 0);
-        exp = i;
-] = cast(ushort)((0x8000 & vu[EXPONENTPOS
+_SHORT] = cast(ushort)((0x8000 & vu[EXPONENTPOS_SHORT
-    }
-  } else static if (real.mant_dig == 113) {
@@ -416 +451 @@
-                    vl[1] |= 0x0000_8000_0000_0000;  // convert $(NAN)S to $(NAN)Q
-                    exp = int.min;
-
+_SHORT
-                    exp = int.min;
-                } else {   // positive infinity
@@ -423 +458 @@
-            } else {
-                exp = ex - 0x3FFE;
-] = cast(ushort)((0x8000 & vu[EXPONENTPOS
+_SHORT] = cast(ushort)((0x8000 & vu[EXPONENTPOS_SHORT
-            }
-        } else if (*vl==0 && ((vl[1]&0x0000_FFFF_FFFF_FFFF)==0)) {
@@ -439 +474 @@
-        } while ((vl[1]&0x0000_8000_0000_0000)== 0);
-        exp = i;
-
-
-
-
+
+
+
-        assert(0, "Unsupported");
-
-
+
-    if (ex) { // If exponent is non-zero
-        if (ex == EXPMASK) {   // infinity or NaN
@@ -458 +491 @@
-        } else {
-            exp = (ex - 0x3FE0) >>> 4;
-] = (0x8000 & ve[EXPONENTPOS
+_SHORT] = (0x8000 & ve[EXPONENTPOS_SHORT
-        }
-    } else if (!(*vl & 0x7FFF_FFFF_FFFF_FFFF)) {
@@ -466 +499 @@
-        // denormal
-        ushort sgn;
-
+_SHORT
-        *vl &= 0x7FFF_FFFF_FFFF_FFFF;
-
@@ -475 +508 @@
-        } while (*vl > 0);
-        exp = i;
-
+_SHORT
-    }
-  }
@@ -587 +620 @@
-            return y;
-        } else static if (real.mant_dig==64) { // 80-bit reals
-0x7FFF
+EXPMASK
-            if (e == 0x7FFF) {
-                // BUG: should also set the invalid exception
@@ -816 +849 @@
-int isNaN(real x)
-{
-
-
+
-        ulong*  p = cast(ulong *)&x;
-
-
-
-
-
-
-
-
-
-
-
+
+
+
-        ulong*  ps = cast(ulong *)&x;
-        return e == 0x7FFF &&
-
+
+
+
+
+
+
+
+
+
+
+
-  } else {
-      return x!=x;
@@ -860 +894 @@
-{
-    uint *p = cast(uint *)&x;
-
-
-
-
+
+
-}
-
@@ -869 +901 @@
-int isNormal(double d)
-{
-
-
-
-
-
+
+
+
+
+
+
+
-}
-
@@ -879 +913 @@
-int isNormal(real x)
-{
-double.mant_dig
+53
-        return isNormal(cast(double)x);
+    } else static if(real.mant_dig==106) { // doubledouble
+    // doubledouble is normal if the least significant part is normal.
+        return isNormal((cast(double*)&x)[1]);
-    } else static if (real.mant_dig == 64) { // real80
-0x7FFF & (cast(ushort *)&x)[4
+EXPMASK & (cast(ushort *)&x)[EXPONENTPOS_SHORT
-        return (e != 0x7FFF && e!=0);
-    } else static if (real.mant_dig == 113) { // quadruple
-0x7FFF & (cast(ushort *)&x)[7
+EXPMASK & (cast(ushort *)&x)[EXPONENTPOS_SHORT
-        return (e != 0x7FFF && e!=0);
-    }
@@ -922 +959 @@
-    long*   pxs = cast(long *)&x;
-    long*   pys = cast(long *)&y;
-double.mant_dig
+53
-    return pxs[0] == pys[0];
-  } else static if (real.mant_dig == 113 || real.mant_dig==106) {
@@ -1006 +1043 @@
-/// ditto
-
-e
-
-double.mant_dig
-e
+x
+
+53
+x
-    } else static if (real.mant_dig == 113) { // quadruple        
-0x7FFF & (cast(ushort *)&x)[7
-e
-ps[0]!=0 || (ps[1]& 0x0000_FFFF_FFFF_FFFF
+EXPMASK & (cast(ushort *)&x)[EXPONENTPOS_SHORT
+x
+((ps[MANTISSA_LSB]|(ps[MANTISSA_MSB]& 0x0000_FFFF_FFFF_FFFF))
-    } else { // real80
-e
-e
-
-4] & 0x7FFF
+x
+x
+
+EXPONENTPOS_SHORT] & EXPMASK
-    }
-}
@@ -1037 +1074 @@
-int isZero(real x)
-{
-double.mant_dig) {
+53) { // double
-        return ((*cast(ulong *)&x) & 0x7FFF_FFFF_FFFF_FFFF) == 0;
-    } else static if (real.mant_dig == 113) { // quadruple   
-        long*   ps = cast(long *)&e;
-return (ps[0]==0 || (ps[1
-
+    return (ps[MANTISSA_LSB]==0 || (ps[MANTISSA_MSB
+ // real80
-        ushort* pe = cast(ushort *)&x;
-        ulong*  ps = cast(ulong  *)&x;
-4] & 0x7FFF
+EXPONENTPOS_SHORT] & EXPMASK
-    }
-}
@@ -1065 +1102 @@
-int isInfinity(real x)
-{
-double.mant_dig) {
+53) { // double
-        return ((*cast(ulong *)&x)&0x7FFF_FFFF_FFFF_FFFF) == 0x7FF8_0000_0000_0000;
+    } else static if(real.mant_dig == 106) { //doubledouble
+        return (((cast(ulong *)&x)[MANTISSA_MSB])&0x7FFF_FFFF_FFFF_FFFF) == 0x7FF8_0000_0000_0000;   
-    } else static if (real.mant_dig == 113) { // quadruple   
-        long*   ps = cast(long *)&x;
-
+
+
-    } else { // real80
-0x7FFF & (cast(ushort *)&x)[4
+EXPMASK & (cast(ushort *)&x)[EXPONENTPOS_SHORT
-        ulong*  ps = cast(ulong *)&x;
-
@@ -1111 +1151 @@
-real nextUp(real x)
-{
-double.mant_dig) {
+53) { // double
-        return nextDoubleUp(x);
-
-
-
-
-
+
+
+
+
-             return x; // +Inf and NaN are unchanged.
-        }     
-        ulong*   ps = cast(ulong *)&e;
-1
-0]==0 && ps[1
-*ps
-1
+MANTISSA_LSB
+MANTISSA_LSB]==0 && ps[MANTISSA_MSB
+ps[MANTISSA_LSB]
+MANTISSA_MSB
-                return x;
-            }
-            --*ps;
-*ps==0) --ps[1
+ps[MANTISSA_LSB]==0) --ps[MANTISSA_MSB
-        } else { // Positive number
-*ps
-*ps==0) ++ps[1
+ps[MANTISSA_LSB]
+ps[MANTISSA_LSB]==0) ++ps[MANTISSA_MSB
-        }
-        return x;
-          
-
+ // real80
-        // For 80-bit reals, the "implied bit" is a nuisance...
-        ushort *pe = cast(ushort *)&x;
-        ulong  *ps = cast(ulong  *)&x;
-
-4] & 0x7FFF
+EXPONENTPOS_SHORT] & EXPMASK
-            // First, deal with NANs and infinity
-            if (x == -real.infinity) return -real.max;
-            return x; // +Inf and NaN are unchanged.
-        }
-4
+EXPONENTPOS_SHORT
-            --*ps;
-            // Need to mask with 0x7FFF... so denormals are treated correctly.
-            if ((*ps & 0x7FFFFFFFFFFFFFFF) == 0x7FFFFFFFFFFFFFFF) {
-4
-4
+EXPONENTPOS_SHORT
+EXPONENTPOS_SHORT
-                    return x;
-                }
-4
-4
+EXPONENTPOS_SHORT
+EXPONENTPOS_SHORT
-                    return x; // it's become a denormal, implied bit stays low.
-                }
@@ -1165 +1204 @@
-            // Works automatically for positive zero.
-            ++*ps;
-FFFFFFFF
+_FFFF_FFFF_
-                // change in exponent
-4
-00000000
+EXPONENTPOS_SHORT
+_0000_0000_
-            }
-        }
@@ -1292 +1331 @@
-        uint *ps = cast(uint *)&x;
-        (*ps) &= 0xFFFF_FC00;
-double.mant_dig
+53
-        ulong *ps = cast(ulong *)&x;
-        (*ps) &= 0xFFFF_FFFF_FC00_0000;
@@ -1303 +1342 @@
-    } else static if (X.mant_dig==113) { // quadruple
-        ulong *ps = cast(ulong *)&x;
-(*ps)
+ps[MANTISSA_LSB]
-    }
-    //else static assert(0, "Unsupported size");
@@ -1430 +1469 @@
-    // they could have 0 or 1 bits in common.
-
-
-
-
-
-
-
-
-
-
-
-
-
+
+
+
+
-    {   // Difference is denormal
-        // For denormals, we need to add the number of zeros that
@@ -1447 +1478 @@
-        // We do this by multiplying by 2^real.mant_dig
-        diff *= POW2MANTDIG;
-
+_SHORT
-    }
-
@@ -1454 +1485 @@
-
-    // Avoid out-by-1 errors when factor is almost 2.
-
-
-
-
-
-
-
-
-
-
-
+
+
+
+
+
+
+
+
-    {   // Difference is denormal
-        // For denormals, we need to add the number of zeros that
-        // lie at the start of diff's significand.
-        // We do this by multiplying by 2^real.mant_dig
-0x1p+53
-
+POW2MANTDIG
+_SHORT
-    }
-
@@ -1477 +1505 @@
-
-    // Avoid out-by-1 errors when factor is almost 2.
-] ^ pb[EXPONENTPOS])&0x7FF0
+_SHORT] ^ pb[EXPONENTPOS_SHORT])&EXPMASK
-    else return 0;
- }
@@ -1535 +1563 @@
-int signbit(real x)
-{
-
-
-
-
-
-
+
-}
-
@@ -1562 +1585 @@
-real copysign(real to, real from)
-{
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
+
+
+
+
+
+
-}
-
@@ -1642 +1656 @@
-        ulong m = ((*xl) & 0x7FFF_FFFF_FFFF_FFFF) + ((*yl) & 0x7FFF_FFFF_FFFF_FFFF);
-
-4] & 0x7FFF) + (ye[4
+EXPONENTPOS_SHORT] & 0x7FFF) + (ye[EXPONENTPOS_SHORT
-        if (m & 0x8000_0000_0000_0000) {
-            ++e;
@@ -1655 +1669 @@
-        else *ul = m; // ... unless exponent is 0 (denormal or zero).
-        // Prevent a ridiculous warning (why does (ushort | ushort) get promoted to int???)
-
+
+
+
+
+
+
+
+
-    } else static if (T.mant_dig == double.mant_dig) {
-        ulong *ul = cast(ulong *)&u;
@@ -1711 +1732 @@
-real NaN(ulong payload)
-{
-double.mant_dig
+53
-      ulong v = 2; // no implied bit. quiet bit = 1
-    } else {
@@ -1734 +1755 @@
-    a >>=29;
-
-double.mant_dig
+53
-        v |=0x7FF0_0000_0000_0000;
-        real x;
@@ -1740 +1761 @@
-        return x;
-    } else {
-        // Extended real bits
-        v <<=11;
-        a &= 0x7FF;
-        v |= a;
-
-        real x = real.nan;
-
+
+
+
+
+
+
+
+
+
+
+
-        return x;
-    }
@@ -1764 +1793 @@
-{
-    assert(isNaN(x));
-ulong m = *cast(ulong *)(&x);
-static if (real.mant_dig == double.mant_dig) {
+static if (real.mant_dig == 53) {
+    ulong m = *cast(ulong *)(&x);
-        // Make it look like an 80-bit significand.
-        // Skip exponent, and quiet bit
-        m &= 0x0007_FFFF_FFFF_FFFF;
-        m <<= 10;
+    } else static if (real.mant_dig==113) { // quadruple
+        version(LittleEndian) {
+            ulong m = *cast(ulong*)(6+cast(ubyte*)(&x));
+        } else {
+            ulong m = *cast(ulong*)(2+cast(ubyte*)(&x));
+        }
+        m>>=1; // there's no implicit bit
+    } else {
+        ulong m = *cast(ulong *)(&x);
-    }
-    // ignore implicit bit and quiet bit
@@ -1782 +1820 @@
-unittest {
-  real nan4 = NaN(0x789_ABCD_EF12_3456);
-
+ || real.mant_dig==113
-      assert (getNaNPayload(nan4) == 0x789_ABCD_EF12_3456);
-  } else {