diff options
Diffstat (limited to 'runtime/c/i64_udivmod.c')
-rw-r--r-- | runtime/c/i64_udivmod.c | 158 |
1 files changed, 158 insertions, 0 deletions
diff --git a/runtime/c/i64_udivmod.c b/runtime/c/i64_udivmod.c new file mode 100644 index 00000000..d8f5073a --- /dev/null +++ b/runtime/c/i64_udivmod.c @@ -0,0 +1,158 @@ +/***************************************************************** + * + * The Compcert verified compiler + * + * Xavier Leroy, INRIA Paris-Rocquencourt + * + * Copyright (c) 2013 Institut National de Recherche en Informatique et + * en Automatique. + * + * Redistribution and use in source and binary forms, with or without + * modification, are permitted provided that the following conditions are met: + * * Redistributions of source code must retain the above copyright + * notice, this list of conditions and the following disclaimer. + * * 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. + * * Neither the name of the <organization> nor the + * names of its contributors may be used to endorse or promote products + * derived from this software without specific prior written permission. + * + * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS 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 <COPYRIGHT + * HOLDER> 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. + * + **********************************************************************/ + +/* Helper functions for 64-bit integer arithmetic. Reference C implementation */ + +#include <stddef.h> +#include "i64.h" + +static unsigned __i64_udiv6432(unsigned u1, unsigned u0, + unsigned v, unsigned *r); +static int __i64_nlz(unsigned x); + +/* Unsigned division and remainder */ + +unsigned long long __i64_udivmod(unsigned long long n, + unsigned long long d, + unsigned long long * rp) +{ + unsigned dh = d >> 32; + if (dh == 0) { + unsigned nh = n >> 32; + if (nh == 0) { + /* Special case 32 / 32 */ + unsigned nl = n; + unsigned dl = d; + *rp = (unsigned long long) (nl % dl); + return (unsigned long long) (nl / dl); + } else { + /* Special case 64 / 32 */ + unsigned nl = n; + unsigned dl = d; + unsigned qh = nh / dl; + unsigned rl; + unsigned ql = __i64_udiv6432(nh % dl, nl, dl, &rl); + *rp = (unsigned long long) rl; /* high word of remainder is 0 */ + return ((unsigned long long) qh) << 32 | ql; + } + } else { + /* General case 64 / 64 */ + unsigned dl = d; + /* Scale N and D down, giving N' and D' such that 2^31 <= D' < 2^32 */ + int s = 32 - __i64_nlz(dh); /* shift amount, between 1 and 32 */ + unsigned long long np = __i64_shr(n, s); + unsigned dp = (unsigned) __i64_shr(d, s); + /* Divide N' by D' to get an approximate quotient Q */ + unsigned q = __i64_udiv6432(np >> 32, np, dp, NULL); + again: ; + /* Tentative quotient Q is either correct or one too high */ + /* Compute Q * D, checking for overflow */ + unsigned long long p1 = (unsigned long long) dl * q; + unsigned long long p2 = (unsigned long long) (dh * q) << 32; + unsigned long long p = p1 + p2; + if (p < p1) { + /* Overflow occurred: adjust Q down by 1 and redo check */ + q = q - 1; goto again; + } + /* Compute remainder R */ + unsigned long long r = n - p; + if (n < p) { + /* Q is one too high, adjust Q down by 1 and R up by D */ + q = q - 1; r = r + d; + } + *rp = r; + return (unsigned long long) q; + } +} + +/* Unsigned division and remainder for 64 bits divided by 32 bits. */ +/* This is algorithm "divlu" from _Hacker's Delight_, fig 9.3 */ + +static unsigned __i64_udiv6432(unsigned u1, unsigned u0, + unsigned v, unsigned *r) +{ + const unsigned b = 65536; // Number base (16 bits). + unsigned un1, un0, // Norm. dividend LSD's. + vn1, vn0, // Norm. divisor digits. + q1, q0, // Quotient digits. + un32, un21, un10,// Divided digit pairs. + rhat; // A remainder. + int s; // Shift amount for norm. + + if (u1 >= v) { // If overflow, + if (r != NULL) *r = 0xFFFFFFFFU; // set rem to an impossible value, + return 0xFFFFFFFFU; // and return largest possible quotient. + } + s = __i64_nlz(v); // 0 <= s <= 31. + v = v << s; // Normalize divisor. + vn1 = v >> 16; // Break divisor up into + vn0 = v & 0xFFFF; // two 16-bit digits. + un32 = (u1 << s) | ((u0 >> (32 - s)) & (-s >> 31)); + un10 = u0 << s; // Shift dividend left. + un1 = un10 >> 16; // Break right half of + un0 = un10 & 0xFFFF; // dividend into two digits. + q1 = un32/vn1; // Compute the first quotient digit, q1. + rhat = un32 - q1*vn1; + again1: + if (q1 >= b || q1*vn0 > b*rhat + un1) { + q1 = q1 - 1; + rhat = rhat + vn1; + if (rhat < b) goto again1; + } + un21 = un32*b + un1 - q1*v; // Multiply and subtract. + q0 = un21/vn1; // Compute the second quotient digit, q0 + rhat = un21 - q0*vn1; +again2: + if (q0 >= b || q0*vn0 > b*rhat + un0) { + q0 = q0 - 1; + rhat = rhat + vn1; + if (rhat < b) goto again2; + } + if (r != NULL) *r = (un21*b + un0 - q0*v) >> s; + return q1*b + q0; +} + +/* Number of leading zeroes */ + +static int __i64_nlz(unsigned x) +{ + if (x == 0) return 32; + int n = 0; + if (x <= 0x0000FFFF) { n += 16; x <<= 16; } + if (x <= 0x00FFFFFF) { n += 8; x <<= 8; } + if (x <= 0x0FFFFFFF) { n += 4; x <<= 4; } + if (x <= 0x3FFFFFFF) { n += 2; x <<= 2; } + if (x <= 0x7FFFFFFF) { n += 1; x <<= 1; } + return n; +} |