mpdecimal-cpython.patch

   1 From 90930e65455f60216f09d175586139242dbba260 Mon Sep 17 00:00:00 2001
   2 From: Stefan Krah <skrah@bytereef.org>
   3 Date: Fri, 21 Feb 2020 01:52:47 +0100
   4 Subject: [PATCH] bpo-39576: Prevent memory error for overly optimistic
   5  precisions (GH-18581)
   6
   7 ---
   8  Lib/test/test_decimal.py              |  35 +++++++
   9  Modules/_decimal/libmpdec/mpdecimal.c |  77 +++++++++++++-
  10  Modules/_decimal/tests/deccheck.py    | 139 +++++++++++++++++++++++++-
  11  3 files changed, 245 insertions(+), 6 deletions(-)
  12
  13 #diff --git a/Modules/_decimal/libmpdec/mpdecimal.c b/Modules/_decimal/libmpdec/mpdecimal.c
  14 --- libmpdec/mpdecimal.c
  15 +++ libmpdec/mpdecimal.c
  16 @@ -3781,6 +3781,43 @@ mpd_qdiv(mpd_t *q, const mpd_t *a, const mpd_t *b,
  17           const mpd_context_t *ctx, uint32_t *status)
  18  {
  19      _mpd_qdiv(SET_IDEAL_EXP, q, a, b, ctx, status);
  20 +
  21 +    if (*status & MPD_Malloc_error) {
  22 +        /* Inexact quotients (the usual case) fill the entire context precision,
  23 +         * which can lead to malloc() failures for very high precisions. Retry
  24 +         * the operation with a lower precision in case the result is exact.
  25 +         *
  26 +         * We need an upper bound for the number of digits of a_coeff / b_coeff
  27 +         * when the result is exact.  If a_coeff' * 1 / b_coeff' is in lowest
  28 +         * terms, then maxdigits(a_coeff') + maxdigits(1 / b_coeff') is a suitable
  29 +         * bound.
  30 +         *
  31 +         * 1 / b_coeff' is exact iff b_coeff' exclusively has prime factors 2 or 5.
  32 +         * The largest amount of digits is generated if b_coeff' is a power of 2 or
  33 +         * a power of 5 and is less than or equal to log5(b_coeff') <= log2(b_coeff').
  34 +         *
  35 +         * We arrive at a total upper bound:
  36 +         *
  37 +         *   maxdigits(a_coeff') + maxdigits(1 / b_coeff') <=
  38 +         *   a->digits + log2(b_coeff) =
  39 +         *   a->digits + log10(b_coeff) / log10(2) <=
  40 +         *   a->digits + b->digits * 4;
  41 +         */
  42 +        uint32_t workstatus = 0;
  43 +        mpd_context_t workctx = *ctx;
  44 +        workctx.prec = a->digits + b->digits * 4;
  45 +        if (workctx.prec >= ctx->prec) {
  46 +            return;  /* No point in retrying, keep the original error. */
  47 +        }
  48 +
  49 +        _mpd_qdiv(SET_IDEAL_EXP, q, a, b, &workctx, &workstatus);
  50 +        if (workstatus == 0) { /* The result is exact, unrounded, normal etc. */
  51 +            *status = 0;
  52 +            return;
  53 +        }
  54 +
  55 +        mpd_seterror(q, *status, status);
  56 +    }
  57  }
  58
  59  /* Internal function. */
  60 @@ -7702,9 +7739,9 @@ mpd_qinvroot(mpd_t *result, const mpd_t *a, const mpd_context_t *ctx,
  61  /* END LIBMPDEC_ONLY */
  62
  63  /* Algorithm from decimal.py */
  64 -void
  65 -mpd_qsqrt(mpd_t *result, const mpd_t *a, const mpd_context_t *ctx,
  66 -          uint32_t *status)
  67 +static void
  68 +_mpd_qsqrt(mpd_t *result, const mpd_t *a, const mpd_context_t *ctx,
  69 +           uint32_t *status)
  70  {
  71      mpd_context_t maxcontext;
  72      MPD_NEW_STATIC(c,0,0,0,0);
  73 @@ -7836,6 +7873,40 @@ mpd_qsqrt(mpd_t *result, const mpd_t *a, const mpd_context_t *ctx,
  74      goto out;
  75  }
  76
  77 +void
  78 +mpd_qsqrt(mpd_t *result, const mpd_t *a, const mpd_context_t *ctx,
  79 +          uint32_t *status)
  80 +{
  81 +    _mpd_qsqrt(result, a, ctx, status);
  82 +
  83 +    if (*status & (MPD_Malloc_error|MPD_Division_impossible)) {
  84 +        /* The above conditions can occur at very high context precisions
  85 +         * if intermediate values get too large. Retry the operation with
  86 +         * a lower context precision in case the result is exact.
  87 +         *
  88 +         * If the result is exact, an upper bound for the number of digits
  89 +         * is the number of digits in the input.
  90 +         *
  91 +         * NOTE: sqrt(40e9) = 2.0e+5 /\ digits(40e9) = digits(2.0e+5) = 2
  92 +         */
  93 +        uint32_t workstatus = 0;
  94 +        mpd_context_t workctx = *ctx;
  95 +        workctx.prec = a->digits;
  96 +
  97 +        if (workctx.prec >= ctx->prec) {
  98 +            return; /* No point in repeating this, keep the original error. */
  99 +        }
 100 +
 101 +        _mpd_qsqrt(result, a, &workctx, &workstatus);
 102 +        if (workstatus == 0) {
 103 +            *status = 0;
 104 +            return;
 105 +        }
 106 +
 107 +        mpd_seterror(result, *status, status);
 108 +    }
 109 +}
 110 +
 111
 112  /******************************************************************************/
 113  /*                              Base conversions                              */