[net.lang.c] Request for multiply/divide routines.

baba@tymix.UUCP (Duane Hentrich) (03/14/86)

Greetings net people.

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.

Pointers to algorithms, flow charts, and even actual C code are all 
appreciated.  Email is preferred.  Will post responses if there is 
enough interest.

Thank you for your kind attention,

Baba Rum Dudu	(aka Duane Hentrich)

doug@terak.UUCP (Doug Pardee) (03/17/86)

> 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.

NS32000 Series assembler:

	MEID	OPND32,R0		;Result in R0 and R1

What?  C???  Oh.  I forgot that C compilers always produce code that is
at least as good as you can write by hand in assembler  :-)
-- 
Doug Pardee -- CalComp -- {elrond,savax,seismo,decvax,ihnp4}!terak!doug

jsdy@hadron.UUCP (Joseph S. D. Yao) (03/19/86)

In article <695@tymix.UUCP> baba@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@seismo.{CSS.GOV,ARPA,UUCP}

gwyn@brl-smoke.ARPA (Doug Gwyn ) (03/20/86)

Knuth Vol. 2 contains algorithms for multiple-precision arithmetic.