[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