dead code in div_q.c?

[paul zimmermann]

       Niels,

> From: nisse at lysator.liu.se (Niels Möller)
> Cc: gmp-devel at gmplib.org,  raphael.rieu-helft at lri.fr
> Date: Wed, 25 Apr 2018 22:22:41 +0200
> User-Agent: Gnus/5.13 (Gnus v5.13) Emacs/25.3 (berkeley-unix)
> 
> nisse at lysator.liu.se (Niels Möller) writes:
> 
> > paul zimmermann <Paul.Zimmermann at inria.fr> writes:
> >
> >> together with Raphaël Rieu-Hleft (in cc), we believe we have found some dead code in
> >> mpn/generic/div_q.c around lines 173-182:
> >>
> >>           else if (UNLIKELY (qh != 0))
> >>             {
> >>               /* This happens only when the quotient is close to B^n and        
> >>                  mpn_*_divappr_q returned B^n.  */
> >
> > I think your right. This comment is wrong in two ways: First, divappr is
> > not called (directly) for this branch of code, and second, the quotient
> > can be at most close to B^n/2.
> 
> I'm trying the appended patch. Besides the analysis, I've ran 
> 
>   while GMP_CHECK_RANDOMIZE=0 ./tests/mpn/t-div ; do : ; done
> 
> for at least half an hour. So I think we can drop the code. 

for the first case qh != 0 (lines 174-183) we are sure this cannot happen,
since cnt < GMP_NUMB_BITS, thus new_dp[new_nn-1] has its most significant
bit 0, and since new_dp[dn-1] has its most significant set, we cannot have
a carry in the division.

In the second case (lines 267-276) this is less obvious since you use an
approximate division (divappr), and I don't know what is the maximal allowed
difference between divappr and div. If divappr always return a quotient
less or equal to the true quotient, then qh cannot be zero in this case.
But if qh can be larger than the true quotient, even by one bit, then it
is possible that qh is 1, in the case {new_np, new_nn} = 0111...111 and
{new_dp, dn} = 1000...000.

> Question is,
> how much of comments and/or ASSERT are useful? And should we make the
> assignments
> 
>   qp[qn - 1] = qh;
> 
> and 
> 
>   tp[qn] = qh;
> 
> unconditional? As far as I see, the areas are always large enough.

beware that in the case cy <> 0, one divides {new_np, nn+1} by {new_dp, dn}
thus filling {qp, nn-dn+1}. If you write qp[qn-1]=qh afterwards, you will
erase the high limb of the quotient. (By the way, this is another way to
prove that qh=0 in this case, since we know the quotient fits on nn-dn+1
limbs.)

If the division routines already check for zero high limbs of the dividend,
then it would be simpler to have new_nn = nn + 1 in all cases, and we would
not need to write qp[qn-1].

Paul

> But it's a unclear to me if and why it's sometimes permissible to omit
> writing the top limbs? Reads of the top limbs look unconditional. 
> 
> I've tried adding 
> 
> --- a/tests/mpn/t-div.c	Wed Apr 25 07:38:14 2018 +0200
> +++ b/tests/mpn/t-div.c	Wed Apr 25 22:19:17 2018 +0200
> @@ -445,6 +445,7 @@ main (int argc, char **argv)
>  	      alloc = itch + 1;
>  	    }
>  	  scratch[itch] = ran;
> +	  MPN_COPY (qp, junkp, nn - dn + 1);
>  	  mpn_div_q (qp, np, nn, dup, dn, scratch);
>  	  ASSERT_ALWAYS (ran == scratch[itch]);
>  	  ASSERT_ALWAYS (qp[-1] == qran0);  ASSERT_ALWAYS (qp[nn - dn + 1] == qran1);
> 
> to the test, to ensure that limbs aren't correct thanks to the previous
> call. Tests still pass. So I'm a bit puzzled.
> 
> Regards,
> /Niels
> 
> *** /tmp/extdiff.0aBNfO/gmp.765c2c27523b/mpn/generic/div_q.c	2018-04-25 21:55:51.457769871 +0200
> --- /home/nisse/hack/gmp/mpn/generic/div_q.c	2018-04-25 21:18:01.516501993 +0200
> *************** mpn_div_q (mp_ptr qp,
> *** 171,184 ****
>   	  if (cy == 0)
>   	    qp[qn - 1] = qh;
> ! 	  else if (UNLIKELY (qh != 0))
> ! 	    {
> ! 	      /* This happens only when the quotient is close to B^n and
> ! 		 mpn_*_divappr_q returned B^n.  */
> ! 	      mp_size_t i, n;
> ! 	      n = new_nn - dn;
> ! 	      for (i = 0; i < n; i++)
> ! 		qp[i] = GMP_NUMB_MAX;
> ! 	      qh = 0;		/* currently ignored */
> ! 	    }
>   	}
>         else  /* divisor is already normalised */
> --- 171,176 ----
>   	  if (cy == 0)
>   	    qp[qn - 1] = qh;
> ! 	  else
> ! 	    ASSERT_ALWAYS (qh == 0);
>   	}
>         else  /* divisor is already normalised */
> *************** mpn_div_q (mp_ptr qp,
> *** 263,276 ****
>   	  if (cy == 0)
>   	    tp[qn] = qh;
> ! 	  else if (UNLIKELY (qh != 0))
> ! 	    {
> ! 	      /* This happens only when the quotient is close to B^n and
> ! 		 mpn_*_divappr_q returned B^n.  */
> ! 	      mp_size_t i, n;
> ! 	      n = new_nn - (qn + 1);
> ! 	      for (i = 0; i < n; i++)
> ! 		tp[i] = GMP_NUMB_MAX;
> ! 	      qh = 0;		/* currently ignored */
> ! 	    }
>   	}
>         else  /* divisor is already normalised */
> --- 255,260 ----
>   	  if (cy == 0)
>   	    tp[qn] = qh;
> ! 	  else
> ! 	    ASSERT_ALWAYS (qh == 0);
>   	}
>         else  /* divisor is already normalised */
> 
> 
> -- 
> Niels Möller. PGP-encrypted email is preferred. Keyid 368C6677.
> Internet email is subject to wholesale government surveillance.
>