1) The overflow detection in mpd_qln() has a surprising number of case splits.

   List all of them in the comment.

2) Use the recently stated relative error of _mpd_qln() to generate the
   interval for the exact value of ln(x). See also the comment in mpd_qexp().

1 files changed, 19 insertions, 7 deletions

diff --git a/Modules/_decimal/libmpdec/mpdecimal.c b/Modules/_decimal/libmpdec/mpdecimal.c
index 346e5ee..4b76fcd 100644
--- a/Modules/_decimal/libmpdec/mpdecimal.c
+++ b/Modules/_decimal/libmpdec/mpdecimal.c
@@ -4632,14 +4632,26 @@ mpd_qln(mpd_t *result, const mpd_t *a, const mpd_context_t *ctx,
         _settriple(result, MPD_POS, 0, 0);
         return;
     }
-    /* Check if the result will overflow.
+    /*
+     * Check if the result will overflow (0 < x, x != 1):
+     *   1) log10(x) < 0 iff adjexp(x) < 0
+     *   2) 0 < x /\ x <= y ==> adjexp(x) <= adjexp(y)
+     *   3) 0 < x /\ x != 1 ==> 2 * abs(log10(x)) < abs(log(x))
+     *   4) adjexp(x) <= log10(x) < adjexp(x) + 1
      *
-     * 1) adjexp(a) + 1 > log10(a) >= adjexp(a)
+     * Case adjexp(x) >= 0:
+     *   5) 2 * adjexp(x) < abs(log(x))
+     *   Case adjexp(x) > 0:
+     *     6) adjexp(2 * adjexp(x)) <= adjexp(abs(log(x)))
+     *   Case adjexp(x) == 0:
+     *     mpd_exp_digits(t)-1 == 0 <= emax (the shortcut is not triggered)
      *
-     * 2) |log10(a)| >= adjexp(a), if adjexp(a) >= 0
-     *    |log10(a)| > -adjexp(a)-1, if adjexp(a) < 0
-     *
-     * 3) |log(a)| > 2*|log10(a)|
+     * Case adjexp(x) < 0:
+     *   7) 2 * (-adjexp(x) - 1) < abs(log(x))
+     *   Case adjexp(x) < -1:
+     *     8) adjexp(2 * (-adjexp(x) - 1)) <= adjexp(abs(log(x)))
+     *   Case adjexp(x) == -1:
+     *     mpd_exp_digits(t)-1 == 0 <= emax (the shortcut is not triggered)
      */
     adjexp = mpd_adjexp(a);
     t = (adjexp < 0) ? -adjexp-1 : adjexp;
@@ -4674,7 +4686,7 @@ mpd_qln(mpd_t *result, const mpd_t *a, const mpd_context_t *ctx,
             workctx.prec = prec;
             _mpd_qln(result, a, &workctx, status);
             _ssettriple(&ulp, MPD_POS, 1,
-                        result->exp + result->digits-workctx.prec-1);
+                        result->exp + result->digits-workctx.prec);
 
             workctx.prec = ctx->prec;
             mpd_qadd(&t1, result, &ulp, &workctx, &workctx.status);