7e6cabd06e
Renumber cluase 4 to 3, per what everybody else did when BSD granted them permission to remove clause 3. My insistance on keeping the same numbering for legal reasons is too pedantic, so give up on that point. Submitted by: Jan Schaumann <jschauma@stevens.edu> Pull Request: https://github.com/freebsd/freebsd/pull/96
349 lines
8.5 KiB
C
349 lines
8.5 KiB
C
/*-
|
|
* Copyright (c) 1992, 1993
|
|
* The Regents of the University of California. All rights reserved.
|
|
*
|
|
* This software was developed by the Computer Systems Engineering group
|
|
* at Lawrence Berkeley Laboratory under DARPA contract BG 91-66 and
|
|
* contributed to Berkeley.
|
|
*
|
|
* 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.
|
|
* 3. Neither the name of the University 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 REGENTS 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 REGENTS 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.
|
|
*
|
|
* From: Id: qdivrem.c,v 1.7 1997/11/07 09:20:40 phk Exp
|
|
*/
|
|
|
|
#include <sys/cdefs.h>
|
|
__FBSDID("$FreeBSD$");
|
|
|
|
/*
|
|
* Multiprecision divide. This algorithm is from Knuth vol. 2 (2nd ed),
|
|
* section 4.3.1, pp. 257--259.
|
|
*/
|
|
|
|
#include "quad.h"
|
|
|
|
#define B (1 << HALF_BITS) /* digit base */
|
|
|
|
/* Combine two `digits' to make a single two-digit number. */
|
|
#define COMBINE(a, b) (((u_int)(a) << HALF_BITS) | (b))
|
|
|
|
_Static_assert(sizeof(int) / 2 == sizeof(short),
|
|
"Bitwise functions in libstand are broken on this architecture\n");
|
|
|
|
/* select a type for digits in base B: use unsigned short if they fit */
|
|
typedef unsigned short digit;
|
|
|
|
/*
|
|
* Shift p[0]..p[len] left `sh' bits, ignoring any bits that
|
|
* `fall out' the left (there never will be any such anyway).
|
|
* We may assume len >= 0. NOTE THAT THIS WRITES len+1 DIGITS.
|
|
*/
|
|
static void
|
|
shl(digit *p, int len, int sh)
|
|
{
|
|
int i;
|
|
|
|
for (i = 0; i < len; i++)
|
|
p[i] = LHALF(p[i] << sh) | (p[i + 1] >> (HALF_BITS - sh));
|
|
p[i] = LHALF(p[i] << sh);
|
|
}
|
|
|
|
/*
|
|
* __qdivrem(u, v, rem) returns u/v and, optionally, sets *rem to u%v.
|
|
*
|
|
* We do this in base 2-sup-HALF_BITS, so that all intermediate products
|
|
* fit within u_int. As a consequence, the maximum length dividend and
|
|
* divisor are 4 `digits' in this base (they are shorter if they have
|
|
* leading zeros).
|
|
*/
|
|
u_quad_t
|
|
__qdivrem(uq, vq, arq)
|
|
u_quad_t uq, vq, *arq;
|
|
{
|
|
union uu tmp;
|
|
digit *u, *v, *q;
|
|
digit v1, v2;
|
|
u_int qhat, rhat, t;
|
|
int m, n, d, j, i;
|
|
digit uspace[5], vspace[5], qspace[5];
|
|
|
|
/*
|
|
* Take care of special cases: divide by zero, and u < v.
|
|
*/
|
|
if (vq == 0) {
|
|
/* divide by zero. */
|
|
static volatile const unsigned int zero = 0;
|
|
|
|
tmp.ul[H] = tmp.ul[L] = 1 / zero;
|
|
if (arq)
|
|
*arq = uq;
|
|
return (tmp.q);
|
|
}
|
|
if (uq < vq) {
|
|
if (arq)
|
|
*arq = uq;
|
|
return (0);
|
|
}
|
|
u = &uspace[0];
|
|
v = &vspace[0];
|
|
q = &qspace[0];
|
|
|
|
/*
|
|
* Break dividend and divisor into digits in base B, then
|
|
* count leading zeros to determine m and n. When done, we
|
|
* will have:
|
|
* u = (u[1]u[2]...u[m+n]) sub B
|
|
* v = (v[1]v[2]...v[n]) sub B
|
|
* v[1] != 0
|
|
* 1 < n <= 4 (if n = 1, we use a different division algorithm)
|
|
* m >= 0 (otherwise u < v, which we already checked)
|
|
* m + n = 4
|
|
* and thus
|
|
* m = 4 - n <= 2
|
|
*/
|
|
tmp.uq = uq;
|
|
u[0] = 0;
|
|
u[1] = HHALF(tmp.ul[H]);
|
|
u[2] = LHALF(tmp.ul[H]);
|
|
u[3] = HHALF(tmp.ul[L]);
|
|
u[4] = LHALF(tmp.ul[L]);
|
|
tmp.uq = vq;
|
|
v[1] = HHALF(tmp.ul[H]);
|
|
v[2] = LHALF(tmp.ul[H]);
|
|
v[3] = HHALF(tmp.ul[L]);
|
|
v[4] = LHALF(tmp.ul[L]);
|
|
for (n = 4; v[1] == 0; v++) {
|
|
if (--n == 1) {
|
|
u_int rbj; /* r*B+u[j] (not root boy jim) */
|
|
digit q1, q2, q3, q4;
|
|
|
|
/*
|
|
* Change of plan, per exercise 16.
|
|
* r = 0;
|
|
* for j = 1..4:
|
|
* q[j] = floor((r*B + u[j]) / v),
|
|
* r = (r*B + u[j]) % v;
|
|
* We unroll this completely here.
|
|
*/
|
|
t = v[2]; /* nonzero, by definition */
|
|
q1 = u[1] / t;
|
|
rbj = COMBINE(u[1] % t, u[2]);
|
|
q2 = rbj / t;
|
|
rbj = COMBINE(rbj % t, u[3]);
|
|
q3 = rbj / t;
|
|
rbj = COMBINE(rbj % t, u[4]);
|
|
q4 = rbj / t;
|
|
if (arq)
|
|
*arq = rbj % t;
|
|
tmp.ul[H] = COMBINE(q1, q2);
|
|
tmp.ul[L] = COMBINE(q3, q4);
|
|
return (tmp.q);
|
|
}
|
|
}
|
|
|
|
/*
|
|
* By adjusting q once we determine m, we can guarantee that
|
|
* there is a complete four-digit quotient at &qspace[1] when
|
|
* we finally stop.
|
|
*/
|
|
for (m = 4 - n; u[1] == 0; u++)
|
|
m--;
|
|
for (i = 4 - m; --i >= 0;)
|
|
q[i] = 0;
|
|
q += 4 - m;
|
|
|
|
/*
|
|
* Here we run Program D, translated from MIX to C and acquiring
|
|
* a few minor changes.
|
|
*
|
|
* D1: choose multiplier 1 << d to ensure v[1] >= B/2.
|
|
*/
|
|
d = 0;
|
|
for (t = v[1]; t < B / 2; t <<= 1)
|
|
d++;
|
|
if (d > 0) {
|
|
shl(&u[0], m + n, d); /* u <<= d */
|
|
shl(&v[1], n - 1, d); /* v <<= d */
|
|
}
|
|
/*
|
|
* D2: j = 0.
|
|
*/
|
|
j = 0;
|
|
v1 = v[1]; /* for D3 -- note that v[1..n] are constant */
|
|
v2 = v[2]; /* for D3 */
|
|
do {
|
|
digit uj0, uj1, uj2;
|
|
|
|
/*
|
|
* D3: Calculate qhat (\^q, in TeX notation).
|
|
* Let qhat = min((u[j]*B + u[j+1])/v[1], B-1), and
|
|
* let rhat = (u[j]*B + u[j+1]) mod v[1].
|
|
* While rhat < B and v[2]*qhat > rhat*B+u[j+2],
|
|
* decrement qhat and increase rhat correspondingly.
|
|
* Note that if rhat >= B, v[2]*qhat < rhat*B.
|
|
*/
|
|
uj0 = u[j + 0]; /* for D3 only -- note that u[j+...] change */
|
|
uj1 = u[j + 1]; /* for D3 only */
|
|
uj2 = u[j + 2]; /* for D3 only */
|
|
if (uj0 == v1) {
|
|
qhat = B;
|
|
rhat = uj1;
|
|
goto qhat_too_big;
|
|
} else {
|
|
u_int nn = COMBINE(uj0, uj1);
|
|
qhat = nn / v1;
|
|
rhat = nn % v1;
|
|
}
|
|
while (v2 * qhat > COMBINE(rhat, uj2)) {
|
|
qhat_too_big:
|
|
qhat--;
|
|
if ((rhat += v1) >= B)
|
|
break;
|
|
}
|
|
/*
|
|
* D4: Multiply and subtract.
|
|
* The variable `t' holds any borrows across the loop.
|
|
* We split this up so that we do not require v[0] = 0,
|
|
* and to eliminate a final special case.
|
|
*/
|
|
for (t = 0, i = n; i > 0; i--) {
|
|
t = u[i + j] - v[i] * qhat - t;
|
|
u[i + j] = LHALF(t);
|
|
t = (B - HHALF(t)) & (B - 1);
|
|
}
|
|
t = u[j] - t;
|
|
u[j] = LHALF(t);
|
|
/*
|
|
* D5: test remainder.
|
|
* There is a borrow if and only if HHALF(t) is nonzero;
|
|
* in that (rare) case, qhat was too large (by exactly 1).
|
|
* Fix it by adding v[1..n] to u[j..j+n].
|
|
*/
|
|
if (HHALF(t)) {
|
|
qhat--;
|
|
for (t = 0, i = n; i > 0; i--) { /* D6: add back. */
|
|
t += u[i + j] + v[i];
|
|
u[i + j] = LHALF(t);
|
|
t = HHALF(t);
|
|
}
|
|
u[j] = LHALF(u[j] + t);
|
|
}
|
|
q[j] = qhat;
|
|
} while (++j <= m); /* D7: loop on j. */
|
|
|
|
/*
|
|
* If caller wants the remainder, we have to calculate it as
|
|
* u[m..m+n] >> d (this is at most n digits and thus fits in
|
|
* u[m+1..m+n], but we may need more source digits).
|
|
*/
|
|
if (arq) {
|
|
if (d) {
|
|
for (i = m + n; i > m; --i)
|
|
u[i] = (u[i] >> d) |
|
|
LHALF(u[i - 1] << (HALF_BITS - d));
|
|
u[i] = 0;
|
|
}
|
|
tmp.ul[H] = COMBINE(uspace[1], uspace[2]);
|
|
tmp.ul[L] = COMBINE(uspace[3], uspace[4]);
|
|
*arq = tmp.q;
|
|
}
|
|
|
|
tmp.ul[H] = COMBINE(qspace[1], qspace[2]);
|
|
tmp.ul[L] = COMBINE(qspace[3], qspace[4]);
|
|
return (tmp.q);
|
|
}
|
|
|
|
/*
|
|
* Divide two unsigned quads.
|
|
*/
|
|
|
|
u_quad_t
|
|
__udivdi3(a, b)
|
|
u_quad_t a, b;
|
|
{
|
|
|
|
return (__qdivrem(a, b, (u_quad_t *)0));
|
|
}
|
|
|
|
/*
|
|
* Return remainder after dividing two unsigned quads.
|
|
*/
|
|
u_quad_t
|
|
__umoddi3(a, b)
|
|
u_quad_t a, b;
|
|
{
|
|
u_quad_t r;
|
|
|
|
(void)__qdivrem(a, b, &r);
|
|
return (r);
|
|
}
|
|
|
|
/*
|
|
* Divide two signed quads.
|
|
* ??? if -1/2 should produce -1 on this machine, this code is wrong
|
|
*/
|
|
quad_t
|
|
__divdi3(a, b)
|
|
quad_t a, b;
|
|
{
|
|
u_quad_t ua, ub, uq;
|
|
int neg;
|
|
|
|
if (a < 0)
|
|
ua = -(u_quad_t)a, neg = 1;
|
|
else
|
|
ua = a, neg = 0;
|
|
if (b < 0)
|
|
ub = -(u_quad_t)b, neg ^= 1;
|
|
else
|
|
ub = b;
|
|
uq = __qdivrem(ua, ub, (u_quad_t *)0);
|
|
return (neg ? -uq : uq);
|
|
}
|
|
|
|
/*
|
|
* Return remainder after dividing two signed quads.
|
|
*
|
|
* XXX
|
|
* If -1/2 should produce -1 on this machine, this code is wrong.
|
|
*/
|
|
quad_t
|
|
__moddi3(a, b)
|
|
quad_t a, b;
|
|
{
|
|
u_quad_t ua, ub, ur;
|
|
int neg;
|
|
|
|
if (a < 0)
|
|
ua = -(u_quad_t)a, neg = 1;
|
|
else
|
|
ua = a, neg = 0;
|
|
if (b < 0)
|
|
ub = -(u_quad_t)b;
|
|
else
|
|
ub = b;
|
|
(void)__qdivrem(ua, ub, &ur);
|
|
return (neg ? -ur : ur);
|
|
}
|