mod to symmetric range

[Torbjorn Granlund]

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