Changeset 543

Timestamp:

01/15/08 03:12:46 (11 months ago)

Author:

Don Clugston

Message:

Merging improvements to Tango.Math into Phobos:
* the low level routines which operate on IEEE floating-point numbers now have versions for 64 and 80 bit reals, and for both big-endian and little-endian systems.
* internal functions isPosZero(), isNegZero() removed in favour of the more generally useful isIdentical().
* asm versions of functions which were not implemented by DMD Windows: scalb, lrint.
* added creal expi(real y) which is useful for simultaneous calculation of sin + cos.
* added comments noting the buggy/unimplemented nature of functions fma(), nan() and poly().

Files:

• candidate/phobos/std/math.d (modified) (26 diffs)

candidate/phobos/std/math.d

@@ -67 +67 @@
-private import std.c.math;
-private import std.traits;
+
+/* 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.
+ * Only three 'real' ABIs are currently supported:
+ * 64 bit Big-endian    (eg PowerPC)
+ * 64 bit Little-endian
+ * 80 bit Little-endian, with implied bit (eg x87, Itanium).
+ */
+version(LittleEndian) {
+    static assert(real.mant_dig == 53 || real.mant_dig==64,
+        "Only 64-bit and 80-bit reals are supported for LittleEndian CPUs");
+} else {
+    static assert(real.mant_dig == 53,
+     "Only 64-bit reals are supported for BigEndian CPUs.");
+}
+
-
-class NotImplemented : Error
@@ -145 +162 @@
-unittest
-{
-PosZero(abs(-0.0L)
+Identical(abs(-0.0L), 0.0L
-    assert(isnan(abs(real.nan)));
-    assert(abs(-real.infinity) == real.infinity);
@@ -470 +487 @@
-unittest
-{
-PosZero(asinh(0.0)
-NegZero(asinh(-0.0)
+Identical(asinh(0.0), 0.0
+Identical(asinh(-0.0), -0.0
-    assert(asinh(real.infinity) == real.infinity);
-    assert(asinh(-real.infinity) == -real.infinity);
@@ -502 +519 @@
-unittest
-{
-PosZero(atanh(0.0)
-NegZero(atanh(-0.0)
+Identical(atanh(0.0), 0.0
+Identical(atanh(-0.0),-0.0
-    assert(isnan(atanh(real.nan)));
-    assert(isnan(atanh(-real.infinity))); 
@@ -619 +636 @@
-real expm1(real x)      { return std.c.math.expm1l(x); }
-
+/**
+ * Calculate cos(y) + i sin(y).
+ *
+ * On many CPUs (such as x86), this is a very efficient operation;
+ * almost twice as fast as calculating sin(y) and cos(y)
+ * seperately, and is the preferred method when both are required.
+ */
+creal expi(real y)
+{
+    version(D_InlineAsm_X86)
+    {
+        asm
+        {
+            fld y;
+            fsincos;
+            fxch st(1), st(0);
+        }
+    }
+    else
+    {
+        return cos(y) + sin(y)*1i;
+    }
+}
+
+unittest
+{
+    assert(expi(1.3e5L) == cos(1.3e5L) + sin(1.3e5L) * 1i);
+    assert(expi(0.0L) == 1L + 0.0Li);
+}
-
-/*********************************************************************
@@ -638 +684 @@
- */
-
-
-real frexp(real value, out int exp)
-{
@@ -645 +690 @@
-    uint ex;
-
-
-
-
-
-
-
-
-
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
-        exp = int.min;
-
-
-
+
-        exp = int.min;
-
-
-
+
-        exp = int.max;
-        }
-
-
-
+
-        exp = ex - 0x3FFE;
-
-
-
-
-
+
+
+
-    // value is +-0.0
-    exp = 0;
-
-
-
+
+
-    int i = -0x3FFD;
-
-
-
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
-        i--;
-        *vl <<= 1;
-    } while (*vl > 0);
-    exp = i;
-
+
+
-    }
-    return value;
@@ -700 +781 @@
-    [-1.0,  -.5,    1],
-    [2.0,   .5, 2],
-
-
-
-
-
+
-    [real.infinity,real.infinity,int.max],
-    [-real.infinity,-real.infinity,int.min],
-    [real.nan,real.nan,int.min],
-    [-real.nan,-real.nan,int.min],
-
-    // Don't really support signalling nan's in D
-    //[real.nans,real.nan,int.min],
-    //[-real.nans,-real.nan,int.min],
-    ];
+
-    int i;
-
-
-
+
-    real x = vals[i][0];
-    real e = vals[i][1];
@@ -723 +796 @@
-    int eptr;
-    real v = frexp(x, eptr);
-
-
-
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
-    assert(exp == eptr);
-
-
-
+
+
+
+
-
-/******************************************
@@ -861 +952 @@
-    version (linux)
-    return std.c.math.scalbnl(x, n);
-
+
+
+
+
+
+
+
+
-    throw new NotImplemented("scalbn");
+}
+
+unittest {
+    assert(scalbn(-real.infinity, 5) == -real.infinity);
-}
-
@@ -1113 +1215 @@
- * Rounds x to the nearest integer value, using the current rounding
- * mode.
+ *
+ * This is generally the fastest method to convert a floating-point number
+ * to an integer. Note that the results from this function
+ * depend on the rounding mode, if the fractional part of x is exactly 0.5.
+ * If using the default rounding mode (ties round to even integers)
+ * lrint(4.5) == 4, lrint(5.5)==6.
- */
-long lrint(real x)
@@ -1118 +1226 @@
-    version (linux)
-    return std.c.math.llrintl(x);
+    else version(D_InlineAsm_X86)
+    {
+        long n;
+        asm
+        {
+            fld x;
+            fistp n;
+        }
+        return n;
+    }
-    else
-    throw new NotImplemented("lrint");
@@ -1184 +1302 @@
- */
-
-
-
-
-
+
+
+
+
+
+
+
+
+
+
-
-    return (pe[4] & 0x7FFF) == 0x7FFF &&
-        *ps & 0x7FFFFFFFFFFFFFFF;
+}
-}
-
@@ -1209 +1334 @@
-int isfinite(real e)
-{
+    static if (real.mant_dig == double.mant_dig) {
+        return isfinite(cast(double)e);
+    } else {
-    ushort* pe = cast(ushort *)&e;
-
-    return (pe[4] & 0x7FFF) != 0x7FFF;
+}
-}
-
@@ -1253 +1382 @@
-/// ditto
-
-
-
-
-
+
+
+
+
+
+
+
-
-    return (pe[4] & 0x7FFF) != 0x7FFF && *ps < 0;
+}
-}
-
@@ -1318 +1451 @@
-int issubnormal(real e)
-{
+    static if (real.mant_dig == double.mant_dig) {
+        return isSubnormal(cast(double)e);
+    } else {
-    ushort* pe = cast(ushort *)&e;
-    long*   ps = cast(long *)&e;
@@ -1323 +1459 @@
-    return (pe[4] & 0x7FFF) == 0 && *ps > 0;
-}
+}
-
-unittest
@@ -1338 +1475 @@
-int isinf(real e)
-{
+    static if (real.mant_dig == double.mant_dig) {
+        return ((*cast(ulong *)&x)&0x7FFF_FFFF_FFFF_FFFF) == 0x7FF8_0000_0000_0000;
+    } else {
-    ushort* pe = cast(ushort *)&e;
-    ulong*  ps = cast(ulong *)&e;
-
-    return (pe[4] & 0x7FFF) == 0x7FFF &&
-
+
+
-}
-
@@ -1356 +1497 @@
-
-/*********************************
+ * Is the binary representation of x identical to y?
+ *
+ * Same as ==, except that positive and negative zero are not identical,
+ * and two $(NAN)s are identical if they have the same 'payload'.
+ */
+
+bool isIdentical(real x, real y)
+{
+    long*   pxs = cast(long *)&x;
+    long*   pys = cast(long *)&y;
+  static if (real.mant_dig == double.mant_dig){
+    return pxs[0] == pys[0];
+  } else {
+    ushort* pxe = cast(ushort *)&x;
+    ushort* pye = cast(ushort *)&y;
+    return pxe[4] == pye[4] && pxs[0] == pys[0];
+  }
+}
+
+/*********************************
- * Return 1 if sign bit of e is set, 0 if not.
- */
-
-
-
-
-
-
+
+
+
+
+
+
-    return (pe[9] & 0x80) != 0;
+}
-}
-
@@ -1415 +1578 @@
-/******************************************
- * Creates a quiet NAN with the information from tagp[] embedded in it.
+ *
+ * BUGS: DMD always returns real.nan, ignoring the payload.
- */
-real nan(char[] tagp) { return std.c.math.nanl(toStringz(tagp)); }
@@ -1485 +1650 @@
- * Returns (x * y) + z, rounding only once according to the
- * current rounding mode.
+ *
+ * BUGS: Not currently implemented - rounds twice.
- */
-real fma(real x, real y, real z) { return (x * y) + z; }
@@ -1704 +1871 @@
-}
-
-// Returns true if x is +0.0 (This function is used in unit tests)
-bool isPosZero(real x)
-{
-    return (x == 0) && (signbit(x) == 0);
-}
-
-// Returns true if x is -0.0 (This function is used in unit tests)
-bool isNegZero(real x)
-{
-    return (x == 0) && signbit(x);
-}
-
-/**************************************
- * To what precision is x equal to y?
@@ -1736 +1891 @@
-    /* Public Domain. Author: Don Clugston, 18 Aug 2005.
-     */
-
-    if (x == y)
-    return real.mant_dig; // ensure diff!=0, cope with INF.
@@ -1755 +1909 @@
-    // always 1 lower than we want, except that if bitsdiff==0,
-    // they could have 0 or 1 bits in common.
-
-
-
+
+
+
+
+
+
+
+
+
+
+
+
+
+
-    {   // Difference is denormal
-    // For denormals, we need to add the number of zeros that
-   // lie at the start of diff's mantissa
+        // lie at the start of diff's significand
-    // We do this by multiplying by 2^real.mant_dig
-   diff *= 0x1p+63
-   return bitsdiff + real.mant_dig - pd[4
+        diff *= POW2MANTDIG
+        return bitsdiff + real.mant_dig - pd[EXPONENTPOS
-    }
-
@@ -1770 +1935 @@
-
-    // Avoid out-by-1 errors when factor is almost 2.
-
+
+
+
+
+
+
+
-}
-
@@ -1834 +2005 @@
-    version (Windows)
-    {
+        // BUG: This code assumes a frame pointer in EBP.
-        asm // assembler by W. Bright
-        {