divexact_1 and bdiv_q_1

bodrato at mail.dm.unipi.it bodrato at mail.dm.unipi.it
Thu Jan 20 10:45:13 CET 2011

Ciao Niels,

On Wed, January 19, 2011 7:34 pm, Niels Möller wrote:
> If I remember correctly, it should work like this. Inputs are U (n
> limbs) and d (single limb, and odd!). Then we should get back Q (n

Both divexact_1 and bdiv_q_1 accept an even divisor. They remove any even
factor, then implement (or call) a pi1_bdiv_q_1 routine.

> limbs), defined by Q = U d^{-1} (mod B^n) and and r (single limb) such
> that
>   U = Q * d - B^n r

A comment from mpn/generic/dive1_c reads:
   mpn_divexact_1c could be created, accepting and returning c.  This would
   let a long calculation be done piece by piece.  Currently there's no
   particular need for that, and not returning c means that a final umul can
   be skipped.
As far as I can understand, returning the correct r needs some more work,
but is never used...

> (This is a bit different from the interface of general bdiv, where an
> n/1 division would produce a quotient of n-1 limbs and a remainder

Well... bdiv_q does not return any reminder, and I do not expect any
returned by bdiv_q_1.

> So it seems that at least r == 0 if and only if d divides U, which makes

This sounds reasonable... but there currently are two implementations.
The C version in mpn/generic/bdiv_q_1.c containing:

mpn_bdiv_q_1 (mp_ptr rp, mp_srcptr up, mp_size_t n, mp_limb_t d)
  return mpn_pi1_bdiv_q_1 (rp, up, n, d, di, shift);

mpn_pi1_bdiv_q_1 (mp_ptr rp, mp_srcptr up, mp_size_t n, mp_limb_t d,
		  mp_limb_t di, int shift)
      c = 0;
      for (i = 1; i < n; i++) /* modifies c, but is skipped if n==1 */
   return c;

I mean, if size (n) is one, the returned value is always zero, regardless
of divisibility.

The asm version in mpn/x86_64/bdiv_q_1.asm contains:
	C rdi	rp end
	mov	%rax, (%rdi)
	pop	%rbx

The returned value is the most significant limb of the result, right?

My mpn/x86/bdiv_q_1.asm implementation should return the most significant
limb too. This can be useful for the typical
  cy = bdiv_q_1 ( ..., size, ...);
  size -= (cy == 0);
This could also be extended to bdiv_q, so that one could write
  cy = bdiv_q ( ..., num_size, ..., den_size, ...);
  quotient_size = num_size - den_size  + (cy != 0);
if we desire a coherent interface.



More information about the gmp-devel mailing list