fix #10802: digamma maps +/-0 to -/+inf (#10825)

This adapts std.mathspecial.digamma so that it IEEE 754 way of handling
+/-0.0, i.e., it treats them as values that are infinitesimally
greater/less than zero. digamma(+0) = -inf and digamma(-0) = +inf.

The DDOC and unittests were also updated.

Author: tedgin

std/internal/math/gammafunction.d

@@ -1643,7 +1643,12 @@ real digamma(real x)
     negative = 0;
     nz = 0.0;
 
-    if ( x <= 0.0 )
+    if ( x == 0.0 )
+    {
+        return signbit(x) == 1 ? real.infinity : -real.infinity;
+    }
+
+    if ( x < 0.0 )
     {
         negative = 1;
         q = x;
@@ -1718,6 +1723,10 @@ done:
     assert(digamma(1.0)== -EULERGAMMA);
     assert(feqrel(digamma(0.25), -PI/2 - 3* LN2 - EULERGAMMA) >= real.mant_dig-7);
     assert(feqrel(digamma(1.0L/6), -PI/2 *sqrt(3.0L) - 2* LN2 -1.5*log(3.0L) - EULERGAMMA) >= real.mant_dig-7);
+    assert(digamma(-0.0) == real.infinity);
+    assert(!digamma(nextDown(-0.0)).isNaN());
+    assert(digamma(+0.0) == -real.infinity);
+    assert(!digamma(nextUp(+0.0)).isNaN());
     assert(digamma(-5.0).isNaN());
     assert(feqrel(digamma(2.5), -EULERGAMMA - 2*LN2 + 2.0 + 2.0L/3) >= real.mant_dig-9);
     assert(isIdentical(digamma(NaN(0xABC)), NaN(0xABC)));

std/mathspecial.d

@@ -29,6 +29,7 @@
  *      NAN = $(RED NAN)
  *      SUP = <span style="vertical-align:super;font-size:smaller">$0</span>
  *      GAMMA = &#915;
+ *      PSI = &Psi;
  *      THETA = &theta;
  *      INTEGRAL = &#8747;
  *      INTEGRATE = $(BIG &#8747;<sub>$(SMALL $1)</sub><sup>$2</sup>)
@@ -37,8 +38,10 @@
  *      BIGSUM = $(BIG &Sigma; <sup>$2</sup><sub>$(SMALL $1)</sub>)
  *      CHOOSE = $(BIG &#40;) <sup>$(SMALL $1)</sup><sub>$(SMALL $2)</sub> $(BIG &#41;)
  *      PLUSMN = &plusmn;
+ *      MNPLUS = &mnplus;
  *      INFIN = &infin;
  *      PLUSMNINF = &plusmn;&infin;
+ *      MNPLUSINF = &mnplus;&infin;
  *      PI = &pi;
  *      LT = &lt;
  *      GT = &gt;
@@ -166,19 +169,47 @@ real beta(real x, real y)
     assert(isIdentical(beta(2, NaN(0xABC)), NaN(0xABC)));
 }
 
-/** Digamma function
+/** Digamma function, $(PSI)(x)
  *
- *  The digamma function is the logarithmic derivative of the gamma function.
+ * $(PSI)(x), is the logarithmic derivative of the gamma function, $(GAMMA)(x).
  *
- *  digamma(x) = d/dx logGamma(x)
+ * $(PSI)(x) = $(SUP d)$(SUB /, dx) ln|$(GAMMA)(x)|  (the derivative of `logGamma(x)`)
  *
- *  See_Also: $(LREF logmdigamma), $(LREF logmdigammaInverse).
+ * Params:
+ *   x = the domain value
+ *
+ * Returns:
+ *   It returns $(PSI)(x).
+ *
+ * $(TABLE_SV
+ *   $(SVH x,             digamma(x)   )
+ *   $(SV  integer < 0,   $(NAN)       )
+ *   $(SV  $(PLUSMN)0.0,  $(MNPLUSINF) )
+ *   $(SV  +$(INFIN),     +$(INFIN)    )
+ *   $(SV  -$(INFIN),     $(NAN)       )
+ *   $(SV  $(NAN),        $(NAN)       )
+ * )
+ *
+ * See_Also: $(LREF logmdigamma), $(LREF logmdigammaInverse).
  */
 real digamma(real x)
 {
     return std.internal.math.gammafunction.digamma(x);
 }
 
+///
+@safe unittest
+{
+    const euler = 0.57721_56649_01532_86060_65121L;
+
+    assert(isClose(digamma(1), -euler));
+    assert(digamma(+0.) == -real.infinity);
+    assert(digamma(-0.) == +real.infinity);
+    assert(digamma(+real.infinity) == +real.infinity);
+    assert(isNaN(digamma(-1)));
+    assert(isNaN(digamma(-real.infinity)));
+}
+
 /** Log Minus Digamma function
  *
  *  logmdigamma(x) = log(x) - digamma(x)