* 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
* 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. All advertising materials mentioning features or use of this software
* must display the following acknowledgement:
* This product includes software developed by the University of
* California, Berkeley and its contributors.
* 4. 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
#if defined(LIBC_SCCS) && !defined(lint)
static char sccsid
[] = "@(#)qdivrem.c 8.1 (Berkeley) 6/4/93";
#endif /* LIBC_SCCS and not lint */
* Multiprecision divide. This algorithm is from Knuth vol. 2 (2nd ed),
* section 4.3.1, pp. 257--259.
#define B (1 << HALF_BITS) /* digit base */
/* Combine two `digits' to make a single two-digit number. */
#define COMBINE(a, b) (((u_long)(a) << HALF_BITS) | (b))
/* select a type for digits in base B: use unsigned short if they fit */
#if ULONG_MAX == 0xffffffff && USHRT_MAX >= 0xffff
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.
shl(register digit
*p
, register int len
, register int sh
)
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_long. As a consequence, the maximum length dividend and
* divisor are 4 `digits' in this base (they are shorter if they have
digit uspace
[5], vspace
[5], qspace
[5];
* Take care of special cases: divide by zero, and u < v.
static volatile const unsigned int zero
= 0;
tmp
.ul
[H
] = tmp
.ul
[L
] = 1 / zero
;
* Break dividend and divisor into digits in base B, then
* count leading zeros to determine m and n. When done, we
* u = (u[1]u[2]...u[m+n]) sub B
* v = (v[1]v[2]...v[n]) sub B
* 1 < n <= 4 (if n = 1, we use a different division algorithm)
* m >= 0 (otherwise u < v, which we already checked)
for (n
= 4; v
[1] == 0; v
++) {
u_long rbj
; /* r*B+u[j] (not root boy jim) */
* Change of plan, per exercise 16.
* q[j] = floor((r*B + u[j]) / v),
* We unroll this completely here.
t
= v
[2]; /* nonzero, by definition */
rbj
= COMBINE(u
[1] % t
, u
[2]);
rbj
= COMBINE(rbj
% t
, u
[3]);
rbj
= COMBINE(rbj
% t
, u
[4]);
tmp
.ul
[H
] = COMBINE(q1
, q2
);
tmp
.ul
[L
] = COMBINE(q3
, q4
);
* By adjusting q once we determine m, we can guarantee that
* there is a complete four-digit quotient at &qspace[1] when
for (m
= 4 - n
; u
[1] == 0; u
++)
for (i
= 4 - m
; --i
>= 0;)
* Here we run Program D, translated from MIX to C and acquiring
* D1: choose multiplier 1 << d to ensure v[1] >= B/2.
for (t
= v
[1]; t
< B
/ 2; t
<<= 1)
shl(&u
[0], m
+ n
, d
); /* u <<= d */
shl(&v
[1], n
- 1, d
); /* v <<= d */
v1
= v
[1]; /* for D3 -- note that v[1..n] are constant */
register 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 */
u_long n
= COMBINE(uj0
, uj1
);
while (v2
* qhat
> COMBINE(rhat
, uj2
)) {
* 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
;
t
= (B
- HHALF(t
)) & (B
- 1);
* 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].
for (t
= 0, i
= n
; i
> 0; i
--) { /* D6: add back. */
} 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).
for (i
= m
+ n
; i
> m
; --i
)
LHALF(u
[i
- 1] << (HALF_BITS
- d
));
tmp
.ul
[H
] = COMBINE(uspace
[1], uspace
[2]);
tmp
.ul
[L
] = COMBINE(uspace
[3], uspace
[4]);
tmp
.ul
[H
] = COMBINE(qspace
[1], qspace
[2]);
tmp
.ul
[L
] = COMBINE(qspace
[3], qspace
[4]);