Checking if A - B > D

[Niels Möller]

nisse at lysator.liu.se (Niels Möller) writes:

> When evaluating Jebelean's criterion for matching partial quotient
> sequences (this is in the context of gcd computation), I need to check if
> 
>   A - B < D
> 
> where A > B.
> 
> What's the best way to do this, using mpn functions? A and B are often
> much larger than D, and I'd prefer not using temporaries larger than
> D. In particular, I'd prefer not having to actually compute the
> complete difference A - B.

I ended up with the following function. It's a little too hairy for my
taste, so I wonder if there are any good ways to simplify it.

Regards,
/Niels

/* Checks if a - b < c. Overwrites c. */
int
nhgcd_diff_smaller_p(mp_srcptr ap, mp_size_t asize,
		     mp_srcptr bp, mp_size_t bsize,
		     mp_ptr cp, mp_size_t csize)
{
  mp_limb_t ch;
  mp_limb_t dh;
  unsigned i;

#if 1
  trace("nhgcd_diff_smaller_p: asize = %d, bsize = %d, csize = %d\n",
	asize, bsize, csize);
  trace("  a = %Nd\n", ap, asize);
  trace("  b = %Nd\n", bp, csize);
  trace("  c = %Nd\n", cp, csize);
#endif

  ASSERT(MPN_LEQ_P(bp, bsize, ap, asize));

  if (csize == 0)
    return 0;

  if (asize < csize)
    return 1;

  if (asize > bsize)
    {
      /* #(a - b) >= asize - 1 */
      if (asize > csize + 1)
	return 0;
      else
	{
	  ASSERT(bsize <= csize);
	  ch = mpn_add(cp, cp, csize, bp, bsize);
	  if (ch)
	    {
	      if (asize == csize)
		return 0;
	      ASSERT(asize == csize + 1);
	      if (ap[asize - 1] < ch)
		return 1;
	      else if (ap[asize - 1] > ch)
		return 0;

	      /* The high limb of a is canceled by ch. */
	      asize--;
	    }

	  return MPN_LESS_P(ap, asize, cp, csize);
	}
    }

  /* a and b are of the same size. Equal high limbs cancel out, so
     ignore them */
  while (asize && ap[asize - 1] == bp[asize - 1])
    asize--;

  if (asize < csize)
    return 1;

  /* Now asize = bsize >= csize */
  if (asize == csize)
    {
      ch = mpn_add_n(cp, cp, bp, csize);
      if (ch)
	return 1;

      return MPN_LESS_P(ap, csize, cp, csize);
    }

  /* asize == bsize > csize. */
  if (csize == 0)
    return 0;

  /* Subtract one from c, transforming the check to a - b <= c */
  ASSERT_NOCARRY(mpn_sub_1(cp, cp, csize, 1));
  MPN_NORMALIZE(cp, csize);

  /* Write a = ah W^k + al, b = bh W^k + bl, where c, al, bl < W^k.
     Then a - b <= c is equivalent to

     (ah - bh) W^k <= c + bl - al

     We always have (ah - bh) >= 0, and

     -W^k < c + bl - ah < 2 W^k
  */

  /* We know that ah - bh >= 1. Do we have ah - bh >= 2? */
  dh = ap[asize - 1] - bp[asize - 1];
  if (dh >= 2)
    return 0;

  /* Compute c += bl - al */
  ch = mpn_add_n(cp, cp, bp, csize);
  if (ch == 0)
    /* c + bl < W^k */
    return 0;

  ch -= mpn_sub_n(cp, cp, ap, csize);

  if (ch == 0)
    return 0;

  /* W^k <= c + bl - al < 2 W^k. Then (ah - bh) W^k is smaller iff ah
     = bh + 1, which happens iff ah ends with 0...0, and bh with
     f...f. */

  for (i = csize; i < asize - 1; i++)
    if (ap[i] != 0 || bp[i] != MP_LIMB_T_MAX)
      return 0;

  return 1;
}