dead code in div_q.c?

[paul zimmermann]

       Niels,

I fully agree with your analysis. Small comments below.

> From: nisse at lysator.liu.se (Niels Möller)
> Date: Thu, 26 Apr 2018 22:12:09 +0200
> 
> paul zimmermann <Paul.Zimmermann at inria.fr> writes:
> 
> > 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.
> 
> You're right again... I thought we had more margin.
> 
> According to comments, quotient from sb and dc divappr is correct or
> one too large, while quotient from mu divappr can be at most 4 too
> large.
> 
> And the limit case, N = B^n/2 - 1, D = B^k/2 is 
> 
>    floor ((B^n/2 - 1) / (B^k/2)) = floor(B^{n-k} - B^{-k}/2) 
>                                  = B^{n-k} - 1 
> 
> So even we only a single unit of error, we could get qh == 1.
> 
> I wonder what actual error can be produced by the divappr functions in
> the corner cases that error could spill over into qh?
> 
> In the cy != 0 case, it looks like we we make a balanced divappr call,
> (2qn+2) / (qn+1). Lets put k = qn+1. For divappr we then have
> 
>   N <= B^{2k}/2 - 1, 
>   D >= B^k/2, lets say D = B/2 + e, e >= 0

D = B/2 + e should read D = B^k/2 + e

> To have any possibility to get qh > 0 with from the approximation error,
> we must have Q >= B^k - 4, where Q denotes the correct quotient. So
> assume that is the case. Now,
> 
>   D <= N/Q <= (B^{2k}/2 - 1) / (B^k - 4)
> 
> which implies
> 
>   e <= N/Q - B^k/2 <= (2 B^k - 1) / (B^k - 4) < 3
> 
> So we need to consider only D = B^k/2 + {0, 1, 2}. Then, any inverse not
> taking the very least significant D limb into account (including mu inverse)
> must be *all* ones, right?

yes if you compute the inverse as floor((B^(2k-2)-1)/floor(D/B))

> It would be good to either prove that qh > 0 can't happen, or add a test
> that will exercise the code handling that case.

yes more generally it would be nice to know for the divappr functions, if the
exact quotient fits on qn limbs (i.e., the exact qh is 0), can the
"approximate" qh be 1?

> > beware that in the case cy <> 0, one divides {new_np, nn+1} by {new_dp, dn}
> > thus filling {qp, nn-dn+1}. 
> 
> Thanks for pointing this out. Now I understand better.
> 
> Regards,
> /Niels
> 
> -- 
> Niels Möller. PGP-encrypted email is preferred. Keyid 368C6677.
> Internet email is subject to wholesale government surveillance.

Paul