Request for multiply/divide routines.
Joseph S. D. Yao
jsdy at hadron.UUCP
Wed Mar 19 15:22:28 AEST 1986
In article <695 at tymix.UUCP> baba at tymix.UUCP (Duane Hentrich) writes:
>I am looking for two routines/algorithms:
>1. multiply two 32-bit unsigned integers giving a 64-bit result
> stored in two 32-bit variables.
>2. divide a 64-bit integer, contained in two 32-bit variables,
> by a 32-bit divisor giving a 32-bit signed remainder and a 32-bit
> quotient.
>Baba Rum Dudu (aka Duane Hentrich)
To: seismo!hplabs!oliveb!tymix!baba
Here is a rough sketch of the algorithms I mentioned:
long long [64-bit, in r0:r1] llmul(mulr, muld)
long [32-bit] int mulr;
long int muld;
{
mask r2, r3, r4, r5;
r0:r1 = 0LL; /* Future product */
r3 = mulr; /* Extended multiplier */
r2 = (r3 < 0 ? -1 : 0);
(* i.e. r2:r3 = mulr *)
r4 = muld; /* Non-extended multiplicand */
r5 = 1;
/*
** This will test bits in r4 (multiplicand). For each
** bit that is set, an r2:r3 shifted left that many bits
** is added into the product. This algorithm works
** regardless of the signs of the operands.
*/
while (r5 != 0) {
if (r4 & r5) { /* is the muld bit set? */
/* if so, add in mulr shifted just so. */
r1 += r3;
addc r0;
r0 += r2;
(* i.e. r0:r1 += r2:r3 *)
}
(r2:r3) <<= 1; /* here mulr is shifted. */
r5 <<= 1; /* also shift bit. */
/*
** ASSUMPTION: bit shifts left leaving no trailing
** bits, and disappears after 32d shift.
*/
}
return(r0:r1);
}
lldiv/llrem is a little mite treackier. Not much, though.
typedef char bool;
typedef int boolean; /* for args & retvals. */
/*
** hide the results of last op: there is just enough repeat
** to make this worth while, given an "expensive" divide.
*/
static long long hide-numerator = 0LL;
static long int hide-denominator = 0L;
static long int hide-quotient = 0L;
static long int hide-remainder = 0L;
/*
** Note that denom could be ll, but that's not what you asked for.
** Also note that, even with ll/l, both q & r could be ll; but
** again, that's not what you asked for. The following can be
** generalised, perhaps not entirely easily, to a G.P. MPA (multi-
** precision arithmetic) schema.
*/
long int entry llrem(num, denom)
long long num;
long int denom;
{
mask r2, r3, r4, r5, r6, r7;
bool is_div = FALSE;
bool is_neg, num_neg;
goto main_body; /* GASP! */
long int entry lldiv(num, denom)
long long num;
long int denom;
{
mask r2, r3, r4, r5, r6, r7;
bool is_div = TRUE;
bool is_neg, num_neg;
main_body:
is_neg = FALSE;
num_neg = FALSE;
if (num == hide-numerator && denom == hide-denominator)
goto give-answer; /* GASP! GASP! */
/*
** for the benefit of flamers: this is a pseudo-code
** transcription of old PDP-11 assembler code. The
** goto's are in a structured tradition (cf. Knuth,
** "Structured Programming with Goto's," and consider
** how to do it without them). Constructive (!!!)
** suggestions welcomed.
*/
r2 = 0; /* Denominator: gets subtracted */
r3 = denom;
(* r2:r3 = denom; *)
/* Correct and record a negative denominator */
if (r3 < 0L) {
r3 = -r3;
is_neg = !is_neg;
}
/*
** MACHINE-DEPENDENT!
** In the one case where a negative negative is still
** negative, we could get caught with high numbers.
*/
if (r3 < 0L) { /* Only 0x8000 does this. */
hide-remainder = r1 & 0x7fff;
hide-quotient = r0 << 1;
if (r1 & 0x8000)
hide-quotient |= 1;
goto give-answer;
}
/*
** And, of course, division by zero.
*/
if (r3 == 0L) { /* Handle divide-by-zero */
#ifdef TRAP_IT
r0:r1 = r0:r1 / r3; anyway - arithmetic trap
# else
hide-remainder = 0L;
hide-quotient = num < 0 ? NEG_INF : POS_INF;
goto give-answer;
#endif/*TRAP_IT*/
}
r0:r1 = num; /* Numerator: subtract from here */
if (r0:r1 < 0LL) {
r0:r1 = -r0:r1;
is_neg = !is_neg;
num_neg = TRUE;
}
r4:r5 = 1LL;
r6 = 0L;
/*
** Shift both denominator and '1' left until
** denominator >= numerator.
*/
while (r2:r3 < r0:r1) {
r2:r3 <<= 1;
r4:r5 <<= 1;
}
/*
** Now shift back, as long as there is a '1' bit;
** and whenever shifted-denom <= num, subtract it out
** and add the low-order 32 bits of the shifted '1' bit
** to the quotient. Stop also if numerator becomes 0LL.
*/
while (r0:r1 != 0LL && r4:r5 != 0LL) {
/* Will this iteration contribute to the quotient? */
if (r2:r3 <= r0:r1) {
r0:r1 -= r2:r3;
r6 += r5;
}
/* Shift denom and '1' back towards "normal". */
r2:r3 >>= 1;
(*
clip any carry bits: e.g.,
r2 &= 0x7fff;
if r2 was zero last time (store in r7?)
r3 &= 0x7fff;
*)
r4:r5 >>= 1;
(* and clip carries. easier: only one bit aloud. *)
}
/*
** An alternative version of these last two loops dumps
** r4:r5 altogether. It increments r6 (re-named r4) by
** one only, shifts it left each loop, and terminates
** the second loop when r2:r3 <= abs(denom). No test for
** r0:r1 == 0. Not sure which is "better."
*/
ASSERT(r0:r1 < denom);
ASSERT((r4:r5 == 0 && r2:r3 == denom/2) || r0:r1 == 0);
/* but probably meaningless. */
ASSERT(r6 is abs(the quotient) && r1 is abs(the remainder));
/*
** Now we do your favourite sign conversions for remainder
** and quotient. It happens that always-positive is modulus,
** not remainder. But there are folks on the net who w a n t
** modulus instead of remainder. They say it's more mathe-
** matical. As a graduate mathematician, I say it's a matter
** of definition; as a (now) software engineer, I say it's a
** matter of ALU architecture.
**
** Example: mine one of several alternatives
** 13/3 4R1 oh, really. never 5R-2.
** -13/3 -4R-1 -5R2
** 13/-3 -4R1 -5R-2 (i don't think anyone wants.)
** -13/-3 4R-1 5R2
*/
if (num_neg) {
/* get a negative remainder. */
r1 = -r1;
}
/*
** Oh, all right. Here's an alternative code to the above:
** if (num_neg) {
** / * get a negative remainder * /
** r1 = -r1;
** / * produce a positive modulus * /
** if (is_neg)
** r1 += denom;
** else
** r1 -= denom;
** / * adjust quotient appropriately * /
** ++r6;
** }
** Any others you will have to figure out on your own.
*/
if (is_neg) /* Little or no question on this. */
r6 = -r6;
hide-quotient = r6;
hide-remainder = r1;
give-answer:
return(is_div ? hide-quotient : hide-remainder);
}
--
Joe Yao hadron!jsdy at seismo.{CSS.GOV,ARPA,UUCP}
More information about the Comp.lang.c
mailing list