[pypy-commit] pypy decimal-libmpdec: Add libmpdec source code, and start a _decimal module.
amauryfa
noreply at buildbot.pypy.org
Sun May 11 00:27:31 CEST 2014
Author: Amaury Forgeot d'Arc <amauryfa at gmail.com>
Branch: decimal-libmpdec
Changeset: r71445:ffe162a71565
Date: 2014-05-01 21:01 +0200
http://bitbucket.org/pypy/pypy/changeset/ffe162a71565/
Log: Add libmpdec source code, and start a _decimal module.
diff too long, truncating to 2000 out of 16459 lines
diff --git a/pypy/module/_decimal/__init__.py b/pypy/module/_decimal/__init__.py
new file mode 100644
--- /dev/null
+++ b/pypy/module/_decimal/__init__.py
@@ -0,0 +1,9 @@
+from pypy.interpreter.mixedmodule import MixedModule
+
+class Module(MixedModule):
+ appleveldefs = {
+ }
+
+ interpleveldefs = {
+ 'IEEE_CONTEXT_MAX_BITS': 'space.wrap(interp_decimal.IEEE_CONTEXT_MAX_BITS)',
+ }
diff --git a/pypy/module/_decimal/interp_decimal.py b/pypy/module/_decimal/interp_decimal.py
new file mode 100644
--- /dev/null
+++ b/pypy/module/_decimal/interp_decimal.py
@@ -0,0 +1,3 @@
+from rpython.rlib import rmpdec
+
+IEEE_CONTEXT_MAX_BITS = rmpdec.MPD_IEEE_CONTEXT_MAX_BITS
diff --git a/pypy/module/_decimal/test/test_module.py b/pypy/module/_decimal/test/test_module.py
new file mode 100644
--- /dev/null
+++ b/pypy/module/_decimal/test/test_module.py
@@ -0,0 +1,6 @@
+class AppTestDecimalModule:
+ spaceconfig = dict(usemodules=('_decimal',))
+
+ def test_constants(self):
+ import _decimal
+ assert _decimal.IEEE_CONTEXT_MAX_BITS > 3
diff --git a/rpython/rlib/rmpdec.py b/rpython/rlib/rmpdec.py
new file mode 100644
--- /dev/null
+++ b/rpython/rlib/rmpdec.py
@@ -0,0 +1,43 @@
+import py
+import sys
+
+from rpython.translator.tool.cbuild import ExternalCompilationInfo
+from rpython.rtyper.tool import rffi_platform as platform
+from rpython.conftest import cdir
+
+libdir = py.path.local(cdir).join('src', 'libmpdec')
+
+compile_extra = []
+if sys.maxsize > 1<<32:
+ compile_extra.append("-DCONFIG_64")
+else:
+ compile_extra.append("-DCONFIG_32")
+
+eci = ExternalCompilationInfo(
+ includes=['src/libmpdec/mpdecimal.h'],
+ include_dirs=[cdir],
+ separate_module_files=[libdir.join('mpdecimal.c'),
+ libdir.join('basearith.c'),
+ libdir.join('convolute.c'),
+ libdir.join('constants.c'),
+ libdir.join('context.c'),
+ libdir.join('fourstep.c'),
+ libdir.join('sixstep.c'),
+ libdir.join('transpose.c'),
+ libdir.join('difradix2.c'),
+ libdir.join('numbertheory.c'),
+ libdir.join('fnt.c'),
+ libdir.join('crt.c'),
+ libdir.join('memory.c'),
+ ],
+ compile_extra=compile_extra,
+ libraries=['m'],
+ )
+
+class CConfig:
+ _compilation_info_ = eci
+
+ MPD_IEEE_CONTEXT_MAX_BITS = platform.ConstantInteger(
+ 'MPD_IEEE_CONTEXT_MAX_BITS')
+
+globals().update(platform.configure(CConfig))
diff --git a/rpython/rlib/test/test_rmpdec.py b/rpython/rlib/test/test_rmpdec.py
new file mode 100644
--- /dev/null
+++ b/rpython/rlib/test/test_rmpdec.py
@@ -0,0 +1,1 @@
+from rpython.rlib import rmpdec
diff --git a/rpython/translator/c/src/libmpdec/README-pypy.txt b/rpython/translator/c/src/libmpdec/README-pypy.txt
new file mode 100644
--- /dev/null
+++ b/rpython/translator/c/src/libmpdec/README-pypy.txt
@@ -0,0 +1,3 @@
+This libmpdec directory was directly copied from CPython.
+
+pyconfig.h was added, with a default configuration which works on Linux.
diff --git a/rpython/translator/c/src/libmpdec/README.txt b/rpython/translator/c/src/libmpdec/README.txt
new file mode 100644
--- /dev/null
+++ b/rpython/translator/c/src/libmpdec/README.txt
@@ -0,0 +1,90 @@
+
+
+libmpdec
+========
+
+libmpdec is a fast C/C++ library for correctly-rounded arbitrary precision
+decimal floating point arithmetic. It is a complete implementation of
+Mike Cowlishaw/IBM's General Decimal Arithmetic Specification.
+
+
+Files required for the Python _decimal module
+=============================================
+
+ Core files for small and medium precision arithmetic
+ ----------------------------------------------------
+
+ basearith.{c,h} -> Core arithmetic in base 10**9 or 10**19.
+ bits.h -> Portable detection of least/most significant one-bit.
+ constants.{c,h} -> Constants that are used in multiple files.
+ context.c -> Context functions.
+ io.{c,h} -> Conversions between mpd_t and ASCII strings,
+ mpd_t formatting (allows UTF-8 fill character).
+ memory.{c,h} -> Allocation handlers with overflow detection
+ and functions for switching between static
+ and dynamic mpd_t.
+ mpdecimal.{c,h} -> All (quiet) functions of the specification.
+ typearith.h -> Fast primitives for double word multiplication,
+ division etc.
+
+ Visual Studio only:
+ ~~~~~~~~~~~~~~~~~~~
+ vccompat.h -> snprintf <==> sprintf_s and similar things.
+ vcstdint.h -> stdint.h (included in VS 2010 but not in VS 2008).
+ vcdiv64.asm -> Double word division used in typearith.h. VS 2008 does
+ not allow inline asm for x64. Also, it does not provide
+ an intrinsic for double word division.
+
+ Files for bignum arithmetic:
+ ----------------------------
+
+ The following files implement the Fast Number Theoretic Transform
+ used for multiplying coefficients with more than 1024 words (see
+ mpdecimal.c: _mpd_fntmul()).
+
+ umodarith.h -> Fast low level routines for unsigned modular arithmetic.
+ numbertheory.{c,h} -> Routines for setting up the Number Theoretic Transform.
+ difradix2.{c,h} -> Decimation in frequency transform, used as the
+ "base case" by the following three files:
+
+ fnt.{c,h} -> Transform arrays up to 4096 words.
+ sixstep.{c,h} -> Transform larger arrays of length 2**n.
+ fourstep.{c,h} -> Transform larger arrays of length 3 * 2**n.
+
+ convolute.{c,h} -> Fast convolution using one of the three transform
+ functions.
+ transpose.{c,h} -> Transpositions needed for the sixstep algorithm.
+ crt.{c,h} -> Chinese Remainder Theorem: use information from three
+ transforms modulo three different primes to get the
+ final result.
+
+
+Pointers to literature, proofs and more
+=======================================
+
+ literature/
+ -----------
+
+ REFERENCES.txt -> List of relevant papers.
+ bignum.txt -> Explanation of the Fast Number Theoretic Transform (FNT).
+ fnt.py -> Verify constants used in the FNT; Python demo for the
+ O(N**2) discrete transform.
+
+ matrix-transform.txt -> Proof for the Matrix Fourier Transform used in
+ fourstep.c.
+ six-step.txt -> Show that the algorithm used in sixstep.c is
+ a variant of the Matrix Fourier Transform.
+ mulmod-64.txt -> Proof for the mulmod64 algorithm from
+ umodarith.h.
+ mulmod-ppro.txt -> Proof for the x87 FPU modular multiplication
+ from umodarith.h.
+ umodarith.lisp -> ACL2 proofs for many functions from umodarith.h.
+
+
+Library Author
+==============
+
+ Stefan Krah <skrah at bytereef.org>
+
+
+
diff --git a/rpython/translator/c/src/libmpdec/basearith.c b/rpython/translator/c/src/libmpdec/basearith.c
new file mode 100644
--- /dev/null
+++ b/rpython/translator/c/src/libmpdec/basearith.c
@@ -0,0 +1,658 @@
+/*
+ * 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;
+}
+
+
+
diff --git a/rpython/translator/c/src/libmpdec/basearith.h b/rpython/translator/c/src/libmpdec/basearith.h
new file mode 100644
--- /dev/null
+++ b/rpython/translator/c/src/libmpdec/basearith.h
@@ -0,0 +1,222 @@
+/*
+ * 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.
+ */
+
+
+#ifndef BASEARITH_H
+#define BASEARITH_H
+
+
+#include "mpdecimal.h"
+#include <stdio.h>
+#include "typearith.h"
+
+
+/* Internal header file: all symbols have local scope in the DSO */
+MPD_PRAGMA(MPD_HIDE_SYMBOLS_START)
+
+
+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);
+void _mpd_baseaddto(mpd_uint_t *w, const mpd_uint_t *u, mpd_size_t n);
+mpd_uint_t _mpd_shortadd(mpd_uint_t *w, mpd_size_t m, mpd_uint_t v);
+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 _mpd_baseincr(mpd_uint_t *u, mpd_size_t n);
+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);
+void _mpd_basesubfrom(mpd_uint_t *w, const mpd_uint_t *u, mpd_size_t n);
+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);
+void _mpd_shortmul(mpd_uint_t *w, const mpd_uint_t *u, mpd_size_t n,
+ mpd_uint_t v);
+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 _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 _mpd_shortdiv(mpd_uint_t *w, const mpd_uint_t *u, mpd_size_t n,
+ mpd_uint_t 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);
+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);
+void _mpd_baseshiftl(mpd_uint_t *dest, mpd_uint_t *src, mpd_size_t n,
+ mpd_size_t m, mpd_size_t shift);
+mpd_uint_t _mpd_baseshiftr(mpd_uint_t *dest, mpd_uint_t *src, mpd_size_t slen,
+ mpd_size_t shift);
+
+
+
+#ifdef CONFIG_64
+extern const mpd_uint_t mprime_rdx;
+
+/*
+ * Algorithm from: Division by Invariant Integers using Multiplication,
+ * T. Granlund and P. L. Montgomery, Proceedings of the SIGPLAN '94
+ * Conference on Programming Language Design and Implementation.
+ *
+ * http://gmplib.org/~tege/divcnst-pldi94.pdf
+ *
+ * Variables from the paper and their translations (See section 8):
+ *
+ * N := 64
+ * d := MPD_RADIX
+ * l := 64
+ * m' := floor((2**(64+64) - 1)/MPD_RADIX) - 2**64
+ *
+ * Since N-l == 0:
+ *
+ * dnorm := d
+ * n2 := hi
+ * n10 := lo
+ *
+ * ACL2 proof: mpd-div-words-r-correct
+ */
+static inline void
+_mpd_div_words_r(mpd_uint_t *q, mpd_uint_t *r, mpd_uint_t hi, mpd_uint_t lo)
+{
+ mpd_uint_t n_adj, h, l, t;
+ mpd_uint_t n1_neg;
+
+ /* n1_neg = if lo >= 2**63 then MPD_UINT_MAX else 0 */
+ n1_neg = (lo & (1ULL<<63)) ? MPD_UINT_MAX : 0;
+ /* n_adj = if lo >= 2**63 then lo+MPD_RADIX else lo */
+ n_adj = lo + (n1_neg & MPD_RADIX);
+
+ /* (h, l) = if lo >= 2**63 then m'*(hi+1) else m'*hi */
+ _mpd_mul_words(&h, &l, mprime_rdx, hi-n1_neg);
+ l = l + n_adj;
+ if (l < n_adj) h++;
+ t = h + hi;
+ /* At this point t == qest, with q == qest or q == qest+1:
+ * 1) 0 <= 2**64*hi + lo - qest*MPD_RADIX < 2*MPD_RADIX
+ */
+
+ /* t = 2**64-1 - qest = 2**64 - (qest+1) */
+ t = MPD_UINT_MAX - t;
+
+ /* (h, l) = 2**64*MPD_RADIX - (qest+1)*MPD_RADIX */
+ _mpd_mul_words(&h, &l, t, MPD_RADIX);
+ l = l + lo;
+ if (l < lo) h++;
+ h += hi;
+ h -= MPD_RADIX;
+ /* (h, l) = 2**64*hi + lo - (qest+1)*MPD_RADIX (mod 2**128)
+ * Case q == qest+1:
+ * a) h == 0, l == r
+ * b) q := h - t == qest+1
+ * c) r := l
+ * Case q == qest:
+ * a) h == MPD_UINT_MAX, l == 2**64-(MPD_RADIX-r)
+ * b) q := h - t == qest
+ * c) r := l + MPD_RADIX = r
+ */
+
+ *q = (h - t);
+ *r = l + (MPD_RADIX & h);
+}
+#else
+static inline void
+_mpd_div_words_r(mpd_uint_t *q, mpd_uint_t *r, mpd_uint_t hi, mpd_uint_t lo)
+{
+ _mpd_div_words(q, r, hi, lo, MPD_RADIX);
+}
+#endif
+
+
+/* Multiply two single base MPD_RADIX words, store result in array w[2]. */
+static inline void
+_mpd_singlemul(mpd_uint_t w[2], mpd_uint_t u, mpd_uint_t v)
+{
+ mpd_uint_t hi, lo;
+
+ _mpd_mul_words(&hi, &lo, u, v);
+ _mpd_div_words_r(&w[1], &w[0], hi, lo);
+}
+
+/* Multiply u (len 2) and v (len m, 1 <= m <= 2). */
+static inline void
+_mpd_mul_2_le2(mpd_uint_t w[4], mpd_uint_t u[2], mpd_uint_t v[2], mpd_ssize_t m)
+{
+ mpd_uint_t hi, lo;
+
+ _mpd_mul_words(&hi, &lo, u[0], v[0]);
+ _mpd_div_words_r(&w[1], &w[0], hi, lo);
+
+ _mpd_mul_words(&hi, &lo, u[1], v[0]);
+ lo = w[1] + lo;
+ if (lo < w[1]) hi++;
+ _mpd_div_words_r(&w[2], &w[1], hi, lo);
+ if (m == 1) return;
+
+ _mpd_mul_words(&hi, &lo, u[0], v[1]);
+ lo = w[1] + lo;
+ if (lo < w[1]) hi++;
+ _mpd_div_words_r(&w[3], &w[1], hi, lo);
+
+ _mpd_mul_words(&hi, &lo, u[1], v[1]);
+ lo = w[2] + lo;
+ if (lo < w[2]) hi++;
+ lo = w[3] + lo;
+ if (lo < w[3]) hi++;
+ _mpd_div_words_r(&w[3], &w[2], hi, lo);
+}
+
+
+/*
+ * Test if all words from data[len-1] to data[0] are zero. If len is 0, nothing
+ * is tested and the coefficient is regarded as "all zero".
+ */
+static inline int
+_mpd_isallzero(const mpd_uint_t *data, mpd_ssize_t len)
+{
+ while (--len >= 0) {
+ if (data[len] != 0) return 0;
+ }
+ return 1;
+}
+
+/*
+ * Test if all full words from data[len-1] to data[0] are MPD_RADIX-1
+ * (all nines). Return true if len == 0.
+ */
+static inline int
+_mpd_isallnine(const mpd_uint_t *data, mpd_ssize_t len)
+{
+ while (--len >= 0) {
+ if (data[len] != MPD_RADIX-1) return 0;
+ }
+ return 1;
+}
+
+
+MPD_PRAGMA(MPD_HIDE_SYMBOLS_END) /* restore previous scope rules */
+
+
+#endif /* BASEARITH_H */
+
+
+
diff --git a/rpython/translator/c/src/libmpdec/bits.h b/rpython/translator/c/src/libmpdec/bits.h
new file mode 100644
--- /dev/null
+++ b/rpython/translator/c/src/libmpdec/bits.h
@@ -0,0 +1,192 @@
+/*
+ * 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.
+ */
+
+
+#ifndef BITS_H
+#define BITS_H
+
+
+#include "mpdecimal.h"
+#include <stdio.h>
+
+
+/* Check if n is a power of 2. */
+static inline int
+ispower2(mpd_size_t n)
+{
+ return n != 0 && (n & (n-1)) == 0;
+}
+
+#if defined(ANSI)
+/*
+ * Return the most significant bit position of n from 0 to 31 (63).
+ * Assumptions: n != 0.
+ */
+static inline int
+mpd_bsr(mpd_size_t n)
+{
+ int pos = 0;
+ mpd_size_t tmp;
+
+#ifdef CONFIG_64
+ tmp = n >> 32;
+ if (tmp != 0) { n = tmp; pos += 32; }
+#endif
+ tmp = n >> 16;
+ if (tmp != 0) { n = tmp; pos += 16; }
+ tmp = n >> 8;
+ if (tmp != 0) { n = tmp; pos += 8; }
+ tmp = n >> 4;
+ if (tmp != 0) { n = tmp; pos += 4; }
+ tmp = n >> 2;
+ if (tmp != 0) { n = tmp; pos += 2; }
+ tmp = n >> 1;
+ if (tmp != 0) { n = tmp; pos += 1; }
+
+ return pos + (int)n - 1;
+}
+
+/*
+ * Return the least significant bit position of n from 0 to 31 (63).
+ * Assumptions: n != 0.
+ */
+static inline int
+mpd_bsf(mpd_size_t n)
+{
+ int pos;
+
+#ifdef CONFIG_64
+ pos = 63;
+ if (n & 0x00000000FFFFFFFFULL) { pos -= 32; } else { n >>= 32; }
+ if (n & 0x000000000000FFFFULL) { pos -= 16; } else { n >>= 16; }
+ if (n & 0x00000000000000FFULL) { pos -= 8; } else { n >>= 8; }
+ if (n & 0x000000000000000FULL) { pos -= 4; } else { n >>= 4; }
+ if (n & 0x0000000000000003ULL) { pos -= 2; } else { n >>= 2; }
+ if (n & 0x0000000000000001ULL) { pos -= 1; }
+#else
+ pos = 31;
+ if (n & 0x000000000000FFFFUL) { pos -= 16; } else { n >>= 16; }
+ if (n & 0x00000000000000FFUL) { pos -= 8; } else { n >>= 8; }
+ if (n & 0x000000000000000FUL) { pos -= 4; } else { n >>= 4; }
+ if (n & 0x0000000000000003UL) { pos -= 2; } else { n >>= 2; }
+ if (n & 0x0000000000000001UL) { pos -= 1; }
+#endif
+ return pos;
+}
+/* END ANSI */
+
+#elif defined(ASM)
+/*
+ * Bit scan reverse. Assumptions: a != 0.
+ */
+static inline int
+mpd_bsr(mpd_size_t a)
+{
+ mpd_size_t retval;
+
+ __asm__ (
+#ifdef CONFIG_64
+ "bsrq %1, %0\n\t"
+#else
+ "bsr %1, %0\n\t"
+#endif
+ :"=r" (retval)
+ :"r" (a)
+ :"cc"
+ );
+
+ return (int)retval;
+}
+
+/*
+ * Bit scan forward. Assumptions: a != 0.
+ */
+static inline int
+mpd_bsf(mpd_size_t a)
+{
+ mpd_size_t retval;
+
+ __asm__ (
+#ifdef CONFIG_64
+ "bsfq %1, %0\n\t"
+#else
+ "bsf %1, %0\n\t"
+#endif
+ :"=r" (retval)
+ :"r" (a)
+ :"cc"
+ );
+
+ return (int)retval;
+}
+/* END ASM */
+
+#elif defined(MASM)
+#include <intrin.h>
+/*
+ * Bit scan reverse. Assumptions: a != 0.
+ */
+static inline int __cdecl
+mpd_bsr(mpd_size_t a)
+{
+ unsigned long retval;
+
+#ifdef CONFIG_64
+ _BitScanReverse64(&retval, a);
+#else
+ _BitScanReverse(&retval, a);
+#endif
+
+ return (int)retval;
+}
+
+/*
+ * Bit scan forward. Assumptions: a != 0.
+ */
+static inline int __cdecl
+mpd_bsf(mpd_size_t a)
+{
+ unsigned long retval;
+
+#ifdef CONFIG_64
+ _BitScanForward64(&retval, a);
+#else
+ _BitScanForward(&retval, a);
+#endif
+
+ return (int)retval;
+}
+/* END MASM (_MSC_VER) */
+#else
+ #error "missing preprocessor definitions"
+#endif /* BSR/BSF */
+
+
+#endif /* BITS_H */
+
+
+
diff --git a/rpython/translator/c/src/libmpdec/constants.c b/rpython/translator/c/src/libmpdec/constants.c
new file mode 100644
--- /dev/null
+++ b/rpython/translator/c/src/libmpdec/constants.c
@@ -0,0 +1,132 @@
+/*
+ * 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 "constants.h"
+
+
+#if defined(CONFIG_64)
+
+ /* number-theory.c */
+ const mpd_uint_t mpd_moduli[3] = {
+ 18446744069414584321ULL, 18446744056529682433ULL, 18446742974197923841ULL
+ };
+ const mpd_uint_t mpd_roots[3] = {7ULL, 10ULL, 19ULL};
+
+ /* crt.c */
+ const mpd_uint_t INV_P1_MOD_P2 = 18446744055098026669ULL;
+ const mpd_uint_t INV_P1P2_MOD_P3 = 287064143708160ULL;
+ const mpd_uint_t LH_P1P2 = 18446744052234715137ULL; /* (P1*P2) % 2^64 */
+ const mpd_uint_t UH_P1P2 = 18446744052234715141ULL; /* (P1*P2) / 2^64 */
+
+ /* transpose.c */
+ const mpd_size_t mpd_bits[64] = {
+ 1, 2, 4, 8, 16, 32, 64, 128, 256, 512, 1024, 2048, 4096, 8192, 16384,
+ 32768, 65536, 131072, 262144, 524288, 1048576, 2097152, 4194304, 8388608,
+ 16777216, 33554432, 67108864, 134217728, 268435456, 536870912, 1073741824,
+ 2147483648ULL, 4294967296ULL, 8589934592ULL, 17179869184ULL, 34359738368ULL,
+ 68719476736ULL, 137438953472ULL, 274877906944ULL, 549755813888ULL,
+ 1099511627776ULL, 2199023255552ULL, 4398046511104, 8796093022208ULL,
+ 17592186044416ULL, 35184372088832ULL, 70368744177664ULL, 140737488355328ULL,
+ 281474976710656ULL, 562949953421312ULL, 1125899906842624ULL,
+ 2251799813685248ULL, 4503599627370496ULL, 9007199254740992ULL,
+ 18014398509481984ULL, 36028797018963968ULL, 72057594037927936ULL,
+ 144115188075855872ULL, 288230376151711744ULL, 576460752303423488ULL,
+ 1152921504606846976ULL, 2305843009213693952ULL, 4611686018427387904ULL,
+ 9223372036854775808ULL
+ };
+
+ /* mpdecimal.c */
+ const mpd_uint_t mpd_pow10[MPD_RDIGITS+1] = {
+ 1,10,100,1000,10000,100000,1000000,10000000,100000000,1000000000,
+ 10000000000ULL,100000000000ULL,1000000000000ULL,10000000000000ULL,
+ 100000000000000ULL,1000000000000000ULL,10000000000000000ULL,
+ 100000000000000000ULL,1000000000000000000ULL,10000000000000000000ULL
+ };
+
+ /* magic number for constant division by MPD_RADIX */
+ const mpd_uint_t mprime_rdx = 15581492618384294730ULL;
+
+#elif defined(CONFIG_32)
+
+ /* number-theory.c */
+ const mpd_uint_t mpd_moduli[3] = {2113929217UL, 2013265921UL, 1811939329UL};
+ const mpd_uint_t mpd_roots[3] = {5UL, 31UL, 13UL};
+
+ /* PentiumPro modular multiplication: These constants have to be loaded as
+ * 80 bit long doubles, which are not supported by certain compilers. */
+ const uint32_t mpd_invmoduli[3][3] = {
+ {4293885170U, 2181570688U, 16352U}, /* ((long double) 1 / 2113929217UL) */
+ {1698898177U, 2290649223U, 16352U}, /* ((long double) 1 / 2013265921UL) */
+ {2716021846U, 2545165803U, 16352U} /* ((long double) 1 / 1811939329UL) */
+ };
+
+ const float MPD_TWO63 = 9223372036854775808.0; /* 2^63 */
+
+ /* crt.c */
+ const mpd_uint_t INV_P1_MOD_P2 = 2013265901UL;
+ const mpd_uint_t INV_P1P2_MOD_P3 = 54UL;
+ const mpd_uint_t LH_P1P2 = 4127195137UL; /* (P1*P2) % 2^32 */
+ const mpd_uint_t UH_P1P2 = 990904320UL; /* (P1*P2) / 2^32 */
+
+ /* transpose.c */
+ const mpd_size_t mpd_bits[32] = {
+ 1, 2, 4, 8, 16, 32, 64, 128, 256, 512, 1024, 2048, 4096, 8192, 16384,
+ 32768, 65536, 131072, 262144, 524288, 1048576, 2097152, 4194304, 8388608,
+ 16777216, 33554432, 67108864, 134217728, 268435456, 536870912, 1073741824,
+ 2147483648UL
+ };
+
+ /* mpdecimal.c */
+ const mpd_uint_t mpd_pow10[MPD_RDIGITS+1] = {
+ 1,10,100,1000,10000,100000,1000000,10000000,100000000,1000000000
+ };
+
+#else
+ #error "CONFIG_64 or CONFIG_32 must be defined."
+#endif
+
+const char *mpd_round_string[MPD_ROUND_GUARD] = {
+ "ROUND_UP", /* round away from 0 */
+ "ROUND_DOWN", /* round toward 0 (truncate) */
+ "ROUND_CEILING", /* round toward +infinity */
+ "ROUND_FLOOR", /* round toward -infinity */
+ "ROUND_HALF_UP", /* 0.5 is rounded up */
+ "ROUND_HALF_DOWN", /* 0.5 is rounded down */
+ "ROUND_HALF_EVEN", /* 0.5 is rounded to even */
+ "ROUND_05UP", /* round zero or five away from 0 */
+ "ROUND_TRUNC", /* truncate, but set infinity */
+};
+
+const char *mpd_clamp_string[MPD_CLAMP_GUARD] = {
+ "CLAMP_DEFAULT",
+ "CLAMP_IEEE_754"
+};
+
+
diff --git a/rpython/translator/c/src/libmpdec/constants.h b/rpython/translator/c/src/libmpdec/constants.h
new file mode 100644
--- /dev/null
+++ b/rpython/translator/c/src/libmpdec/constants.h
@@ -0,0 +1,90 @@
+/*
+ * 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.
+ */
+
+
+#ifndef CONSTANTS_H
+#define CONSTANTS_H
+
+
+#include "mpdecimal.h"
+
+
+/* Internal header file: all symbols have local scope in the DSO */
+MPD_PRAGMA(MPD_HIDE_SYMBOLS_START)
+
+
+/* choice of optimized functions */
+#if defined(CONFIG_64)
+/* x64 */
+ #define MULMOD(a, b) x64_mulmod(a, b, umod)
+ #define MULMOD2C(a0, a1, w) x64_mulmod2c(a0, a1, w, umod)
+ #define MULMOD2(a0, b0, a1, b1) x64_mulmod2(a0, b0, a1, b1, umod)
+ #define POWMOD(base, exp) x64_powmod(base, exp, umod)
+ #define SETMODULUS(modnum) std_setmodulus(modnum, &umod)
+ #define SIZE3_NTT(x0, x1, x2, w3table) std_size3_ntt(x0, x1, x2, w3table, umod)
+#elif defined(PPRO)
+/* PentiumPro (or later) gcc inline asm */
+ #define MULMOD(a, b) ppro_mulmod(a, b, &dmod, dinvmod)
+ #define MULMOD2C(a0, a1, w) ppro_mulmod2c(a0, a1, w, &dmod, dinvmod)
+ #define MULMOD2(a0, b0, a1, b1) ppro_mulmod2(a0, b0, a1, b1, &dmod, dinvmod)
+ #define POWMOD(base, exp) ppro_powmod(base, exp, &dmod, dinvmod)
+ #define SETMODULUS(modnum) ppro_setmodulus(modnum, &umod, &dmod, dinvmod)
+ #define SIZE3_NTT(x0, x1, x2, w3table) ppro_size3_ntt(x0, x1, x2, w3table, umod, &dmod, dinvmod)
+#else
+ /* ANSI C99 */
+ #define MULMOD(a, b) std_mulmod(a, b, umod)
+ #define MULMOD2C(a0, a1, w) std_mulmod2c(a0, a1, w, umod)
+ #define MULMOD2(a0, b0, a1, b1) std_mulmod2(a0, b0, a1, b1, umod)
+ #define POWMOD(base, exp) std_powmod(base, exp, umod)
+ #define SETMODULUS(modnum) std_setmodulus(modnum, &umod)
+ #define SIZE3_NTT(x0, x1, x2, w3table) std_size3_ntt(x0, x1, x2, w3table, umod)
+#endif
+
+/* PentiumPro (or later) gcc inline asm */
+extern const float MPD_TWO63;
+extern const uint32_t mpd_invmoduli[3][3];
+
+enum {P1, P2, P3};
+
+extern const mpd_uint_t mpd_moduli[];
+extern const mpd_uint_t mpd_roots[];
+extern const mpd_size_t mpd_bits[];
+extern const mpd_uint_t mpd_pow10[];
+
+extern const mpd_uint_t INV_P1_MOD_P2;
+extern const mpd_uint_t INV_P1P2_MOD_P3;
+extern const mpd_uint_t LH_P1P2;
+extern const mpd_uint_t UH_P1P2;
+
+
+MPD_PRAGMA(MPD_HIDE_SYMBOLS_END) /* restore previous scope rules */
+
+
+#endif /* CONSTANTS_H */
+
+
+
diff --git a/rpython/translator/c/src/libmpdec/context.c b/rpython/translator/c/src/libmpdec/context.c
new file mode 100644
--- /dev/null
+++ b/rpython/translator/c/src/libmpdec/context.c
@@ -0,0 +1,286 @@
+/*
+ * 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 <string.h>
+#include <signal.h>
+
+
+void
+mpd_dflt_traphandler(mpd_context_t *ctx UNUSED)
+{
+ raise(SIGFPE);
+}
+
+void (* mpd_traphandler)(mpd_context_t *) = mpd_dflt_traphandler;
+
+
+/* Set guaranteed minimum number of coefficient words. The function may
+ be used once at program start. Setting MPD_MINALLOC to out-of-bounds
+ values is a catastrophic error, so in that case the function exits rather
+ than relying on the user to check a return value. */
+void
+mpd_setminalloc(mpd_ssize_t n)
+{
+ static int minalloc_is_set = 0;
+
+ if (minalloc_is_set) {
+ mpd_err_warn("mpd_setminalloc: ignoring request to set "
+ "MPD_MINALLOC a second time\n");
+ return;
+ }
+ if (n < MPD_MINALLOC_MIN || n > MPD_MINALLOC_MAX) {
+ mpd_err_fatal("illegal value for MPD_MINALLOC"); /* GCOV_NOT_REACHED */
+ }
+ MPD_MINALLOC = n;
+ minalloc_is_set = 1;
+}
+
+void
+mpd_init(mpd_context_t *ctx, mpd_ssize_t prec)
+{
+ mpd_ssize_t ideal_minalloc;
+
+ mpd_defaultcontext(ctx);
+
+ if (!mpd_qsetprec(ctx, prec)) {
+ mpd_addstatus_raise(ctx, MPD_Invalid_context);
+ return;
+ }
+
+ ideal_minalloc = 2 * ((prec+MPD_RDIGITS-1) / MPD_RDIGITS);
+ if (ideal_minalloc < MPD_MINALLOC_MIN) ideal_minalloc = MPD_MINALLOC_MIN;
+ if (ideal_minalloc > MPD_MINALLOC_MAX) ideal_minalloc = MPD_MINALLOC_MAX;
+
+ mpd_setminalloc(ideal_minalloc);
+}
+
+void
+mpd_maxcontext(mpd_context_t *ctx)
+{
+ ctx->prec=MPD_MAX_PREC;
+ ctx->emax=MPD_MAX_EMAX;
+ ctx->emin=MPD_MIN_EMIN;
+ ctx->round=MPD_ROUND_HALF_EVEN;
+ ctx->traps=MPD_Traps;
+ ctx->status=0;
+ ctx->newtrap=0;
+ ctx->clamp=0;
+ ctx->allcr=1;
+}
+
+void
+mpd_defaultcontext(mpd_context_t *ctx)
+{
+ ctx->prec=2*MPD_RDIGITS;
+ ctx->emax=MPD_MAX_EMAX;
+ ctx->emin=MPD_MIN_EMIN;
+ ctx->round=MPD_ROUND_HALF_UP;
+ ctx->traps=MPD_Traps;
+ ctx->status=0;
+ ctx->newtrap=0;
+ ctx->clamp=0;
+ ctx->allcr=1;
+}
+
+void
+mpd_basiccontext(mpd_context_t *ctx)
+{
+ ctx->prec=9;
+ ctx->emax=MPD_MAX_EMAX;
+ ctx->emin=MPD_MIN_EMIN;
+ ctx->round=MPD_ROUND_HALF_UP;
+ ctx->traps=MPD_Traps|MPD_Clamped;
+ ctx->status=0;
+ ctx->newtrap=0;
+ ctx->clamp=0;
+ ctx->allcr=1;
+}
+
+int
+mpd_ieee_context(mpd_context_t *ctx, int bits)
+{
+ if (bits <= 0 || bits > MPD_IEEE_CONTEXT_MAX_BITS || bits % 32) {
+ return -1;
+ }
+
+ ctx->prec = 9 * (bits/32) - 2;
+ ctx->emax = 3 * ((mpd_ssize_t)1<<(bits/16+3));
+ ctx->emin = 1 - ctx->emax;
+ ctx->round=MPD_ROUND_HALF_EVEN;
+ ctx->traps=0;
+ ctx->status=0;
+ ctx->newtrap=0;
+ ctx->clamp=1;
+ ctx->allcr=1;
+
+ return 0;
+}
+
+mpd_ssize_t
+mpd_getprec(const mpd_context_t *ctx)
+{
+ return ctx->prec;
+}
+
+mpd_ssize_t
+mpd_getemax(const mpd_context_t *ctx)
+{
+ return ctx->emax;
+}
+
+mpd_ssize_t
+mpd_getemin(const mpd_context_t *ctx)
+{
+ return ctx->emin;
+}
+
+int
+mpd_getround(const mpd_context_t *ctx)
+{
+ return ctx->round;
+}
+
+uint32_t
+mpd_gettraps(const mpd_context_t *ctx)
+{
+ return ctx->traps;
+}
+
+uint32_t
+mpd_getstatus(const mpd_context_t *ctx)
+{
+ return ctx->status;
+}
+
+int
+mpd_getclamp(const mpd_context_t *ctx)
+{
+ return ctx->clamp;
+}
+
+int
+mpd_getcr(const mpd_context_t *ctx)
+{
+ return ctx->allcr;
+}
+
+
+int
+mpd_qsetprec(mpd_context_t *ctx, mpd_ssize_t prec)
+{
+ if (prec <= 0 || prec > MPD_MAX_PREC) {
+ return 0;
+ }
+ ctx->prec = prec;
+ return 1;
+}
+
+int
+mpd_qsetemax(mpd_context_t *ctx, mpd_ssize_t emax)
+{
+ if (emax < 0 || emax > MPD_MAX_EMAX) {
+ return 0;
+ }
+ ctx->emax = emax;
+ return 1;
+}
+
+int
+mpd_qsetemin(mpd_context_t *ctx, mpd_ssize_t emin)
+{
+ if (emin > 0 || emin < MPD_MIN_EMIN) {
+ return 0;
+ }
+ ctx->emin = emin;
+ return 1;
+}
+
+int
+mpd_qsetround(mpd_context_t *ctx, int round)
+{
+ if (!(0 <= round && round < MPD_ROUND_GUARD)) {
+ return 0;
+ }
+ ctx->round = round;
+ return 1;
+}
+
+int
+mpd_qsettraps(mpd_context_t *ctx, uint32_t traps)
+{
+ if (traps > MPD_Max_status) {
+ return 0;
+ }
+ ctx->traps = traps;
+ return 1;
+}
+
+int
+mpd_qsetstatus(mpd_context_t *ctx, uint32_t flags)
+{
+ if (flags > MPD_Max_status) {
+ return 0;
+ }
+ ctx->status = flags;
+ return 1;
+}
+
+int
+mpd_qsetclamp(mpd_context_t *ctx, int c)
+{
+ if (c != 0 && c != 1) {
+ return 0;
+ }
+ ctx->clamp = c;
+ return 1;
+}
+
+int
+mpd_qsetcr(mpd_context_t *ctx, int c)
+{
+ if (c != 0 && c != 1) {
+ return 0;
+ }
+ ctx->allcr = c;
+ return 1;
+}
+
+
+void
+mpd_addstatus_raise(mpd_context_t *ctx, uint32_t flags)
+{
+ ctx->status |= flags;
+ if (flags&ctx->traps) {
+ ctx->newtrap = (flags&ctx->traps);
+ mpd_traphandler(ctx);
+ }
+}
+
+
diff --git a/rpython/translator/c/src/libmpdec/convolute.c b/rpython/translator/c/src/libmpdec/convolute.c
new file mode 100644
--- /dev/null
+++ b/rpython/translator/c/src/libmpdec/convolute.c
@@ -0,0 +1,174 @@
+/*
+ * 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 "bits.h"
+#include "constants.h"
+#include "fnt.h"
+#include "fourstep.h"
+#include "numbertheory.h"
+#include "sixstep.h"
+#include "umodarith.h"
+#include "convolute.h"
+
+
+/* Bignum: Fast convolution using the Number Theoretic Transform. Used for
+ the multiplication of very large coefficients. */
+
+
+/* Convolute the data in c1 and c2. Result is in c1. */
+int
+fnt_convolute(mpd_uint_t *c1, mpd_uint_t *c2, mpd_size_t n, int modnum)
+{
+ int (*fnt)(mpd_uint_t *, mpd_size_t, int);
+ int (*inv_fnt)(mpd_uint_t *, mpd_size_t, int);
+#ifdef PPRO
+ double dmod;
+ uint32_t dinvmod[3];
+#endif
+ mpd_uint_t n_inv, umod;
+ mpd_size_t i;
+
+
+ SETMODULUS(modnum);
+ n_inv = POWMOD(n, (umod-2));
+
+ if (ispower2(n)) {
+ if (n > SIX_STEP_THRESHOLD) {
+ fnt = six_step_fnt;
+ inv_fnt = inv_six_step_fnt;
+ }
+ else {
+ fnt = std_fnt;
+ inv_fnt = std_inv_fnt;
+ }
+ }
+ else {
+ fnt = four_step_fnt;
+ inv_fnt = inv_four_step_fnt;
+ }
+
+ if (!fnt(c1, n, modnum)) {
+ return 0;
+ }
+ if (!fnt(c2, n, modnum)) {
+ return 0;
+ }
+ for (i = 0; i < n-1; i += 2) {
+ mpd_uint_t x0 = c1[i];
+ mpd_uint_t y0 = c2[i];
+ mpd_uint_t x1 = c1[i+1];
+ mpd_uint_t y1 = c2[i+1];
+ MULMOD2(&x0, y0, &x1, y1);
+ c1[i] = x0;
+ c1[i+1] = x1;
+ }
+
+ if (!inv_fnt(c1, n, modnum)) {
+ return 0;
+ }
+ for (i = 0; i < n-3; i += 4) {
+ mpd_uint_t x0 = c1[i];
+ mpd_uint_t x1 = c1[i+1];
+ mpd_uint_t x2 = c1[i+2];
+ mpd_uint_t x3 = c1[i+3];
+ MULMOD2C(&x0, &x1, n_inv);
+ MULMOD2C(&x2, &x3, n_inv);
+ c1[i] = x0;
+ c1[i+1] = x1;
+ c1[i+2] = x2;
+ c1[i+3] = x3;
+ }
+
+ return 1;
+}
+
+/* Autoconvolute the data in c1. Result is in c1. */
+int
+fnt_autoconvolute(mpd_uint_t *c1, mpd_size_t n, int modnum)
+{
+ int (*fnt)(mpd_uint_t *, mpd_size_t, int);
+ int (*inv_fnt)(mpd_uint_t *, mpd_size_t, int);
+#ifdef PPRO
+ double dmod;
+ uint32_t dinvmod[3];
+#endif
+ mpd_uint_t n_inv, umod;
+ mpd_size_t i;
+
+
+ SETMODULUS(modnum);
+ n_inv = POWMOD(n, (umod-2));
+
+ if (ispower2(n)) {
+ if (n > SIX_STEP_THRESHOLD) {
+ fnt = six_step_fnt;
+ inv_fnt = inv_six_step_fnt;
+ }
+ else {
+ fnt = std_fnt;
+ inv_fnt = std_inv_fnt;
+ }
+ }
+ else {
+ fnt = four_step_fnt;
+ inv_fnt = inv_four_step_fnt;
+ }
+
+ if (!fnt(c1, n, modnum)) {
+ return 0;
+ }
+ for (i = 0; i < n-1; i += 2) {
+ mpd_uint_t x0 = c1[i];
+ mpd_uint_t x1 = c1[i+1];
+ MULMOD2(&x0, x0, &x1, x1);
+ c1[i] = x0;
+ c1[i+1] = x1;
+ }
+
+ if (!inv_fnt(c1, n, modnum)) {
+ return 0;
+ }
+ for (i = 0; i < n-3; i += 4) {
+ mpd_uint_t x0 = c1[i];
+ mpd_uint_t x1 = c1[i+1];
+ mpd_uint_t x2 = c1[i+2];
+ mpd_uint_t x3 = c1[i+3];
+ MULMOD2C(&x0, &x1, n_inv);
+ MULMOD2C(&x2, &x3, n_inv);
+ c1[i] = x0;
+ c1[i+1] = x1;
+ c1[i+2] = x2;
+ c1[i+3] = x3;
+ }
+
+ return 1;
+}
+
+
diff --git a/rpython/translator/c/src/libmpdec/convolute.h b/rpython/translator/c/src/libmpdec/convolute.h
new file mode 100644
--- /dev/null
+++ b/rpython/translator/c/src/libmpdec/convolute.h
@@ -0,0 +1,50 @@
+/*
+ * 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
More information about the pypy-commit
mailing list