# Checking if A - B > D

Niels Möller nisse at lysator.liu.se
Sat Sep 13 22:39:15 CEST 2003

```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;
}
```