mod to symmetric range
Torbjorn Granlund
tg at gmplib.org
Tue Sep 7 17:07:27 CEST 2010
Paul Zimmermann <Paul.Zimmermann at loria.fr> writes:
I did skeptically review the code (but not test it). See the attached file.
My comments are prefixed by PZ.
Thanks!
mpz_add (rem, rem, divisor); /* PZ: if B is odd, then the symmetric mod
is fully symmetric, thus we could simply
mod |rem|, and change its sign
afterwards, avoiding an addition.
If B is even, the only case where this
would fail is when rem = -|B|/2, but
this can be detected below. */
I am afraid I don't understand this comment.
I don't understand "we could simply mod |rem|". The addition could be
avoided if we used mpn division, computing mod |divisor|. Perhaps
that's what you're saying?
rp = PTR(rem);
bp = PTR(divisor);
/* Adjust if R > floor(|B|/2) */
if (rn == bn)
{
bhl = bp[bn - 1] >> 1;
if (rp[rn - 1] > bhl)
{
/* R is surely too large, since if rh := rp[rn - 1] and
bh := bp[bn - 1] then R >= rh*beta^(n-1) >= (bh+1)/2*beta^(n-1)
> B/2 */
mpz_sub (rem, rem, divisor);
}
else if (rp[rn - 1] == bhl) /* PZ: since rn = SIZ(rem), we have
rp[rn - 1] >= 1, thus bp[bn - 1] >= 2,
and B/2 has still bn limbs. */
{
mpn_rshift (tp, bp, bn, 1); /* PZ: still bn limbs (see above) */
if (mpn_cmp (rp, tp, bn) > 0) /* PZ: we could compare only the bn-1
lower limbs, since rp[bn-1]=tp[bn-1]
*/
mpz_sub (rem, rem, divisor);
}
/* if rp[rn - 1] < bhl, then R < (rh+1)*beta^(n-1) <= bhl*beta^(n-1)
<= B/2 */
}
else if (rn == bn - 1 && bp[bn - 1] == 1)
{
bhl = GMP_NUMB_HIGHBIT | (bp[bn - 2] >> 1);
/* PZ: since the if-else code is the same as above, maybe we can merge
them? */
The size passed to mpn_cmp is different (bn vs rn), unfortunately. But
your cmp comment implies that we can always use size rn!
if (rp[rn - 1] > bhl)
{
/* R is surely too large, same analysis as above */
mpz_sub (rem, rem, divisor);
}
else if (rp[rn - 1] == bhl)
{
mpn_rshift (tp, bp, bn, 1);
if (mpn_cmp (rp, tp, rn) > 0) /* PZ: same remark as above */
I don't think this cmp can use fewer limbs.
Thanks for nice feedback!
--
Torbjörn
More information about the gmp-discuss
mailing list