divexact_1 and bdiv_q_1

[bodrato at mail.dm.unipi.it]

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:

mp_limb_t
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);
}

mp_limb_t
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
	ret

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.

Regards,
Marco

-- 
http://bodrato.it/papers/