mirrors/cpython
1
0
Fork 0
mirror of https://github.com/python/cpython.git synced 2025-08-02 16:13:13 +00:00
Code Issues Projects Releases Packages Wiki Activity Actions 6
cpython/Modules/_decimal/libmpdec/basearith.c
Stefan Krah c35a8e5c98 Proactive reliability fix for broken FPUs: The base conversion functions 
use log10() to calculate the size of the output array. The current code
has been tested on x86/amd64 (and to a lesser extent on qemu-mips qemu-sparc)
and produces sufficiently large values for all inputs tested so far (coefficient
sizes of 10**18 - 1 are hard to test exhaustively).

The new code does not rely on the correctness of log10() and resizes
the output arrays if the allocated space is insufficient.
2012-06-30 18:05:33 +02:00

658 lines
17 KiB
C
Raw Blame History

/*
* Copyright (c) 2008-2012 Stefan Krah. All rights reserved.
*
* Redistribution and use in source and binary forms, with or without
* modification, are permitted provided that the following conditions
* are met:
*
* 1. Redistributions of source code must retain the above copyright
* notice, this list of conditions and the following disclaimer.
*
* 2. Redistributions in binary form must reproduce the above copyright
* notice, this list of conditions and the following disclaimer in the
* documentation and/or other materials provided with the distribution.
*
* THIS SOFTWARE IS PROVIDED BY THE AUTHOR AND CONTRIBUTORS "AS IS" AND
* ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
* IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
* ARE DISCLAIMED. IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE
* FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
* DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
* OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
* HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
* LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
* OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
* SUCH DAMAGE.
*/
#include "mpdecimal.h"
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <assert.h>
#include "constants.h"
#include "memory.h"
#include "typearith.h"
#include "basearith.h"
/*********************************************************************/
/* Calculations in base MPD_RADIX */
/*********************************************************************/
/*
* Knuth, TAOCP, Volume 2, 4.3.1:
* w := sum of u (len m) and v (len n)
* n > 0 and m >= n
* The calling function has to handle a possible final carry.
*/
mpd_uint_t
_mpd_baseadd(mpd_uint_t *w, const mpd_uint_t *u, const mpd_uint_t *v,
mpd_size_t m, mpd_size_t n)
{
mpd_uint_t s;
mpd_uint_t carry = 0;
mpd_size_t i;
assert(n > 0 && m >= n);
/* add n members of u and v */
for (i = 0; i < n; i++) {
s = u[i] + (v[i] + carry);
carry = (s < u[i]) | (s >= MPD_RADIX);
w[i] = carry ? s-MPD_RADIX : s;
}
/* if there is a carry, propagate it */
for (; carry && i < m; i++) {
s = u[i] + carry;
carry = (s == MPD_RADIX);
w[i] = carry ? 0 : s;
}
/* copy the rest of u */
for (; i < m; i++) {
w[i] = u[i];
}
return carry;
}
/*
* Add the contents of u to w. Carries are propagated further. The caller
* has to make sure that w is big enough.
*/
void
_mpd_baseaddto(mpd_uint_t *w, const mpd_uint_t *u, mpd_size_t n)
{
mpd_uint_t s;
mpd_uint_t carry = 0;
mpd_size_t i;
if (n == 0) return;
/* add n members of u to w */
for (i = 0; i < n; i++) {
s = w[i] + (u[i] + carry);
carry = (s < w[i]) | (s >= MPD_RADIX);
w[i] = carry ? s-MPD_RADIX : s;
}
/* if there is a carry, propagate it */
for (; carry; i++) {
s = w[i] + carry;
carry = (s == MPD_RADIX);
w[i] = carry ? 0 : s;
}
}
/*
* Add v to w (len m). The calling function has to handle a possible
* final carry. Assumption: m > 0.
*/
mpd_uint_t
_mpd_shortadd(mpd_uint_t *w, mpd_size_t m, mpd_uint_t v)
{
mpd_uint_t s;
mpd_uint_t carry;
mpd_size_t i;
assert(m > 0);
/* add v to w */
s = w[0] + v;
carry = (s < v) | (s >= MPD_RADIX);
w[0] = carry ? s-MPD_RADIX : s;
/* if there is a carry, propagate it */
for (i = 1; carry && i < m; i++) {
s = w[i] + carry;
carry = (s == MPD_RADIX);
w[i] = carry ? 0 : s;
}
return carry;
}
/* Increment u. The calling function has to handle a possible carry. */
mpd_uint_t
_mpd_baseincr(mpd_uint_t *u, mpd_size_t n)
{
mpd_uint_t s;
mpd_uint_t carry = 1;
mpd_size_t i;
assert(n > 0);
/* if there is a carry, propagate it */
for (i = 0; carry && i < n; i++) {
s = u[i] + carry;
carry = (s == MPD_RADIX);
u[i] = carry ? 0 : s;
}
return carry;
}
/*
* Knuth, TAOCP, Volume 2, 4.3.1:
* w := difference of u (len m) and v (len n).
* number in u >= number in v;
*/
void
_mpd_basesub(mpd_uint_t *w, const mpd_uint_t *u, const mpd_uint_t *v,
mpd_size_t m, mpd_size_t n)
{
mpd_uint_t d;
mpd_uint_t borrow = 0;
mpd_size_t i;
assert(m > 0 && n > 0);
/* subtract n members of v from u */
for (i = 0; i < n; i++) {
d = u[i] - (v[i] + borrow);
borrow = (u[i] < d);
w[i] = borrow ? d + MPD_RADIX : d;
}
/* if there is a borrow, propagate it */
for (; borrow && i < m; i++) {
d = u[i] - borrow;
borrow = (u[i] == 0);
w[i] = borrow ? MPD_RADIX-1 : d;
}
/* copy the rest of u */
for (; i < m; i++) {
w[i] = u[i];
}
}
/*
* Subtract the contents of u from w. w is larger than u. Borrows are
* propagated further, but eventually w can absorb the final borrow.
*/
void
_mpd_basesubfrom(mpd_uint_t *w, const mpd_uint_t *u, mpd_size_t n)
{
mpd_uint_t d;
mpd_uint_t borrow = 0;
mpd_size_t i;
if (n == 0) return;
/* subtract n members of u from w */
for (i = 0; i < n; i++) {
d = w[i] - (u[i] + borrow);
borrow = (w[i] < d);
w[i] = borrow ? d + MPD_RADIX : d;
}
/* if there is a borrow, propagate it */
for (; borrow; i++) {
d = w[i] - borrow;
borrow = (w[i] == 0);
w[i] = borrow ? MPD_RADIX-1 : d;
}
}
/* w := product of u (len n) and v (single word) */
void
_mpd_shortmul(mpd_uint_t *w, const mpd_uint_t *u, mpd_size_t n, mpd_uint_t v)
{
mpd_uint_t hi, lo;
mpd_uint_t carry = 0;
mpd_size_t i;
assert(n > 0);
for (i=0; i < n; i++) {
_mpd_mul_words(&hi, &lo, u[i], v);
lo = carry + lo;
if (lo < carry) hi++;
_mpd_div_words_r(&carry, &w[i], hi, lo);
}
w[i] = carry;
}
/*
* Knuth, TAOCP, Volume 2, 4.3.1:
* w := product of u (len m) and v (len n)
* w must be initialized to zero
*/
void
_mpd_basemul(mpd_uint_t *w, const mpd_uint_t *u, const mpd_uint_t *v,
mpd_size_t m, mpd_size_t n)
{
mpd_uint_t hi, lo;
mpd_uint_t carry;
mpd_size_t i, j;
assert(m > 0 && n > 0);
for (j=0; j < n; j++) {
carry = 0;
for (i=0; i < m; i++) {
_mpd_mul_words(&hi, &lo, u[i], v[j]);
lo = w[i+j] + lo;
if (lo < w[i+j]) hi++;
lo = carry + lo;
if (lo < carry) hi++;
_mpd_div_words_r(&carry, &w[i+j], hi, lo);
}
w[j+m] = carry;
}
}
/*
* Knuth, TAOCP Volume 2, 4.3.1, exercise 16:
* w := quotient of u (len n) divided by a single word v
*/
mpd_uint_t
_mpd_shortdiv(mpd_uint_t *w, const mpd_uint_t *u, mpd_size_t n, mpd_uint_t v)
{
mpd_uint_t hi, lo;
mpd_uint_t rem = 0;
mpd_size_t i;
assert(n > 0);
for (i=n-1; i != MPD_SIZE_MAX; i--) {
_mpd_mul_words(&hi, &lo, rem, MPD_RADIX);
lo = u[i] + lo;
if (lo < u[i]) hi++;
_mpd_div_words(&w[i], &rem, hi, lo, v);
}
return rem;
}
/*
* Knuth, TAOCP Volume 2, 4.3.1:
* q, r := quotient and remainder of uconst (len nplusm)
* divided by vconst (len n)
* nplusm >= n
*
* If r is not NULL, r will contain the remainder. If r is NULL, the
* return value indicates if there is a remainder: 1 for true, 0 for
* false. A return value of -1 indicates an error.
*/
int
_mpd_basedivmod(mpd_uint_t *q, mpd_uint_t *r,
const mpd_uint_t *uconst, const mpd_uint_t *vconst,
mpd_size_t nplusm, mpd_size_t n)
{
mpd_uint_t ustatic[MPD_MINALLOC_MAX];
mpd_uint_t vstatic[MPD_MINALLOC_MAX];
mpd_uint_t *u = ustatic;
mpd_uint_t *v = vstatic;
mpd_uint_t d, qhat, rhat, w2[2];
mpd_uint_t hi, lo, x;
mpd_uint_t carry;
mpd_size_t i, j, m;
int retval = 0;
assert(n > 1 && nplusm >= n);
m = sub_size_t(nplusm, n);
/* D1: normalize */
d = MPD_RADIX / (vconst[n-1] + 1);
if (nplusm >= MPD_MINALLOC_MAX) {
if ((u = mpd_alloc(nplusm+1, sizeof *u)) == NULL) {
return -1;
}
}
if (n >= MPD_MINALLOC_MAX) {
if ((v = mpd_alloc(n+1, sizeof *v)) == NULL) {
mpd_free(u);
return -1;
}
}
_mpd_shortmul(u, uconst, nplusm, d);
_mpd_shortmul(v, vconst, n, d);
/* D2: loop */
for (j=m; j != MPD_SIZE_MAX; j--) {
/* D3: calculate qhat and rhat */
rhat = _mpd_shortdiv(w2, u+j+n-1, 2, v[n-1]);
qhat = w2[1] * MPD_RADIX + w2[0];
while (1) {
if (qhat < MPD_RADIX) {
_mpd_singlemul(w2, qhat, v[n-2]);
if (w2[1] <= rhat) {
if (w2[1] != rhat || w2[0] <= u[j+n-2]) {
break;
}
}
}
qhat -= 1;
rhat += v[n-1];
if (rhat < v[n-1] || rhat >= MPD_RADIX) {
break;
}
}
/* D4: multiply and subtract */
carry = 0;
for (i=0; i <= n; i++) {
_mpd_mul_words(&hi, &lo, qhat, v[i]);
lo = carry + lo;
if (lo < carry) hi++;
_mpd_div_words_r(&hi, &lo, hi, lo);
x = u[i+j] - lo;
carry = (u[i+j] < x);
u[i+j] = carry ? x+MPD_RADIX : x;
carry += hi;
}
q[j] = qhat;
/* D5: test remainder */
if (carry) {
q[j] -= 1;
/* D6: add back */
(void)_mpd_baseadd(u+j, u+j, v, n+1, n);
}
}
/* D8: unnormalize */
if (r != NULL) {
_mpd_shortdiv(r, u, n, d);
/* we are not interested in the return value here */
retval = 0;
}
else {
retval = !_mpd_isallzero(u, n);
}
if (u != ustatic) mpd_free(u);
if (v != vstatic) mpd_free(v);
return retval;
}
/*
* Left shift of src by 'shift' digits; src may equal dest.
*
* dest := area of n mpd_uint_t with space for srcdigits+shift digits.
* src := coefficient with length m.
*
* The case splits in the function are non-obvious. The following
* equations might help:
*
* Let msdigits denote the number of digits in the most significant
* word of src. Then 1 <= msdigits <= rdigits.
*
* 1) shift = q * rdigits + r
* 2) srcdigits = qsrc * rdigits + msdigits
* 3) destdigits = shift + srcdigits
* = q * rdigits + r + qsrc * rdigits + msdigits
* = q * rdigits + (qsrc * rdigits + (r + msdigits))
*
* The result has q zero words, followed by the coefficient that
* is left-shifted by r. The case r == 0 is trivial. For r > 0, it
* is important to keep in mind that we always read m source words,
* but write m+1 destination words if r + msdigits > rdigits, m words
* otherwise.
*/
void
_mpd_baseshiftl(mpd_uint_t *dest, mpd_uint_t *src, mpd_size_t n, mpd_size_t m,
mpd_size_t shift)
{
#if defined(__GNUC__) && !defined(__INTEL_COMPILER) && !defined(__clang__)
/* spurious uninitialized warnings */
mpd_uint_t l=l, lprev=lprev, h=h;
#else
mpd_uint_t l, lprev, h;
#endif
mpd_uint_t q, r;
mpd_uint_t ph;
assert(m > 0 && n >= m);
_mpd_div_word(&q, &r, (mpd_uint_t)shift, MPD_RDIGITS);
if (r != 0) {
ph = mpd_pow10[r];
--m; --n;
_mpd_divmod_pow10(&h, &lprev, src[m--], MPD_RDIGITS-r);
if (h != 0) { /* r + msdigits > rdigits <==> h != 0 */
dest[n--] = h;
}
/* write m-1 shifted words */
for (; m != MPD_SIZE_MAX; m--,n--) {
_mpd_divmod_pow10(&h, &l, src[m], MPD_RDIGITS-r);
dest[n] = ph * lprev + h;
lprev = l;
}
/* write least significant word */
dest[q] = ph * lprev;
}
else {
while (--m != MPD_SIZE_MAX) {
dest[m+q] = src[m];
}
}
mpd_uint_zero(dest, q);
}
/*
* Right shift of src by 'shift' digits; src may equal dest.
* Assumption: srcdigits-shift > 0.
*
* dest := area with space for srcdigits-shift digits.
* src := coefficient with length 'slen'.
*
* The case splits in the function rely on the following equations:
*
* Let msdigits denote the number of digits in the most significant
* word of src. Then 1 <= msdigits <= rdigits.
*
* 1) shift = q * rdigits + r
* 2) srcdigits = qsrc * rdigits + msdigits
* 3) destdigits = srcdigits - shift
* = qsrc * rdigits + msdigits - (q * rdigits + r)
* = (qsrc - q) * rdigits + msdigits - r
*
* Since destdigits > 0 and 1 <= msdigits <= rdigits:
*
* 4) qsrc >= q
* 5) qsrc == q ==> msdigits > r
*
* The result has slen-q words if msdigits > r, slen-q-1 words otherwise.
*/
mpd_uint_t
_mpd_baseshiftr(mpd_uint_t *dest, mpd_uint_t *src, mpd_size_t slen,
mpd_size_t shift)
{
#if defined(__GNUC__) && !defined(__INTEL_COMPILER) && !defined(__clang__)
/* spurious uninitialized warnings */
mpd_uint_t l=l, h=h, hprev=hprev; /* low, high, previous high */
#else
mpd_uint_t l, h, hprev; /* low, high, previous high */
#endif
mpd_uint_t rnd, rest; /* rounding digit, rest */
mpd_uint_t q, r;
mpd_size_t i, j;
mpd_uint_t ph;
assert(slen > 0);
_mpd_div_word(&q, &r, (mpd_uint_t)shift, MPD_RDIGITS);
rnd = rest = 0;
if (r != 0) {
ph = mpd_pow10[MPD_RDIGITS-r];
_mpd_divmod_pow10(&hprev, &rest, src[q], r);
_mpd_divmod_pow10(&rnd, &rest, rest, r-1);
if (rest == 0 && q > 0) {
rest = !_mpd_isallzero(src, q);
}
/* write slen-q-1 words */
for (j=0,i=q+1; i<slen; i++,j++) {
_mpd_divmod_pow10(&h, &l, src[i], r);
dest[j] = ph * l + hprev;
hprev = h;
}
/* write most significant word */
if (hprev != 0) { /* always the case if slen==q-1 */
dest[j] = hprev;
}
}
else {
if (q > 0) {
_mpd_divmod_pow10(&rnd, &rest, src[q-1], MPD_RDIGITS-1);
/* is there any non-zero digit below rnd? */
if (rest == 0) rest = !_mpd_isallzero(src, q-1);
}
for (j = 0; j < slen-q; j++) {
dest[j] = src[q+j];
}
}
/* 0-4 ==> rnd+rest < 0.5 */
/* 5 ==> rnd+rest == 0.5 */
/* 6-9 ==> rnd+rest > 0.5 */
return (rnd == 0 || rnd == 5) ? rnd + !!rest : rnd;
}
/*********************************************************************/
/* Calculations in base b */
/*********************************************************************/
/*
* Add v to w (len m). The calling function has to handle a possible
* final carry. Assumption: m > 0.
*/
mpd_uint_t
_mpd_shortadd_b(mpd_uint_t *w, mpd_size_t m, mpd_uint_t v, mpd_uint_t b)
{
mpd_uint_t s;
mpd_uint_t carry;
mpd_size_t i;
assert(m > 0);
/* add v to w */
s = w[0] + v;
carry = (s < v) | (s >= b);
w[0] = carry ? s-b : s;
/* if there is a carry, propagate it */
for (i = 1; carry && i < m; i++) {
s = w[i] + carry;
carry = (s == b);
w[i] = carry ? 0 : s;
}
return carry;
}
/* w := product of u (len n) and v (single word). Return carry. */
mpd_uint_t
_mpd_shortmul_c(mpd_uint_t *w, const mpd_uint_t *u, mpd_size_t n, mpd_uint_t v)
{
mpd_uint_t hi, lo;
mpd_uint_t carry = 0;
mpd_size_t i;
assert(n > 0);
for (i=0; i < n; i++) {
_mpd_mul_words(&hi, &lo, u[i], v);
lo = carry + lo;
if (lo < carry) hi++;
_mpd_div_words_r(&carry, &w[i], hi, lo);
}
return carry;
}
/* w := product of u (len n) and v (single word) */
mpd_uint_t
_mpd_shortmul_b(mpd_uint_t *w, const mpd_uint_t *u, mpd_size_t n,
mpd_uint_t v, mpd_uint_t b)
{
mpd_uint_t hi, lo;
mpd_uint_t carry = 0;
mpd_size_t i;
assert(n > 0);
for (i=0; i < n; i++) {
_mpd_mul_words(&hi, &lo, u[i], v);
lo = carry + lo;
if (lo < carry) hi++;
_mpd_div_words(&carry, &w[i], hi, lo, b);
}
return carry;
}
/*
* Knuth, TAOCP Volume 2, 4.3.1, exercise 16:
* w := quotient of u (len n) divided by a single word v
*/
mpd_uint_t
_mpd_shortdiv_b(mpd_uint_t *w, const mpd_uint_t *u, mpd_size_t n,
mpd_uint_t v, mpd_uint_t b)
{
mpd_uint_t hi, lo;
mpd_uint_t rem = 0;
mpd_size_t i;
assert(n > 0);
for (i=n-1; i != MPD_SIZE_MAX; i--) {
_mpd_mul_words(&hi, &lo, rem, b);
lo = u[i] + lo;
if (lo < u[i]) hi++;
_mpd_div_words(&w[i], &rem, hi, lo, v);
}
return rem;
}
Reference in a new issue View git blame Copy permalink