mirrors/cosmopolitan
1
0
Fork 0
mirror of https://github.com/jart/cosmopolitan.git synced 2025-02-12 17:27:56 +00:00
Code Issues Projects Releases Packages Wiki Activity Actions 1
cosmopolitan/third_party/python/Modules/_decimal/libmpdec/basearith.c
Justine Tunney 9b29358511 Make whitespace changes 
Status lines for Emacs and Vim have been added to Python sources so
they'll be easier to edit using Python's preferred coding style.

Some DNS helper functions have been broken up into multiple files. It's
nice to have one function per file whenever possible, since that way we
don't need -ffunction-sections.  Another reason it's good to have small
source files, is because the build will be enforcing resource limits on
compilation and testing soon.
2021-08-13 03:20:45 -07:00

656 lines
18 KiB
C
Raw Blame History

/*-*- mode:c;indent-tabs-mode:nil;c-basic-offset:4;tab-width:8;coding:utf-8 -*-│
│vi: set net ft=c ts=4 sts=4 sw=4 fenc=utf-8 :vi│
╞══════════════════════════════════════════════════════════════════════════════╡
│ Copyright (c) 2008-2016 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 "third_party/python/Modules/_decimal/libmpdec/basearith.h"
#include "third_party/python/Modules/_decimal/libmpdec/constants.h"
#include "third_party/python/Modules/_decimal/libmpdec/mpdecimal.h"
#include "third_party/python/Modules/_decimal/libmpdec/typearith.h"
/* clang-format off */
asm(".ident\t\"\\n\\n\
libmpdec (BSD-2)\\n\
Copyright 2008-2016 Stefan Krah\"");
asm(".include \"libc/disclaimer.inc\"");
/*********************************************************************/
/* 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