I have resolved the issue with the asymptotically fast division code. Fortunately, our code was always returning the correct result. The original test code was not only broken, but it was testing the wrong thing!
The relevant result is in this paper: http://www.loria.fr/~zimmerma/papers/invert.pdf In particular, it guarantees only that AX < B^2n <= A(X + 2), where X is the inverse being computed. In our asymptotically fast division code, we need the more strict AX < B^2n <= A(X + 1), which is what the original test code was supposed to be testing for, but wasn't. With this assumption, by following the analysis of the paper, the value th in the original code is always zero, and cy is always zero. Thus the high bits are never altered by addition of cy, and thus remain correct. Whilst the original code does only guarantee AX < B^2n <= A(X + 2), we then add a further test which checks if the result satisfies the more restricted inequality or not. If not X is incremented by 1 so that the stricter inequality holds. So in short, none of the code in MPIR was returning the wrong results, despite the code and test code being wrong. However, I have now modified both. I removed the unnecessary test of th and all the code dealing with cy, which is always 0 in our case. I fixed the test code in both places that it occurs. All tests still pass. I've committed this to our repository. Bill. On 1 February 2015 at 15:41, Bill Hart <goodwillh...@googlemail.com> wrote: > I was just studying the asymptotically fast division code in MPIR (used > for division in the FFT range), and I have just spotted what I believe is a > serious bug. > > The value returned by mpn/generic/invert.c almost always returns the wrong > result. > > The associated test code for this code checks two inequalities. But after > checking the first inequality, it overwrites the return value (pass/fail) > with the result from the second inequality. Thus the first inequality is > effectively never checked. Therefore, the test code never picks up the fact > that the return value from invert.c is almost always too large. > > Presumably our division code doesn't actually return wrong answers because > of this. It treats the inverse as approximate and adjusts. But one imagines > there is the potential that it occasionally returns incorrect results or > runs more slowly than it ought. > > Of course we have tested the division code extensively, so the likelihood > of hitting such an incorrect value by chance is probably remote. But I felt > it best to disclose the bug in invert.c as soon as possible. > > Note that invert.c almost always returns incorrect results. But it is > presumably not being used directly by anyone. It's not part of the public > interface of MPIR. It's only used for our asymptotically fast division, and > for approximate division (which subsequently gets corrected in any actual > division code anyway). > > Bill. > -- You received this message because you are subscribed to the Google Groups "mpir-devel" group. To unsubscribe from this group and stop receiving emails from it, send an email to mpir-devel+unsubscr...@googlegroups.com. To post to this group, send email to mpir-devel@googlegroups.com. Visit this group at http://groups.google.com/group/mpir-devel. For more options, visit https://groups.google.com/d/optout.
