GMP 4.4 remaining tasks
tg at gmplib.org
Tue Dec 8 08:36:47 CET 2009
bodrato at mail.dm.unipi.it writes:
We should try the "add a virtual limb" strategy... But if you have the
"exact or one unit too small" result and want to correct it... I suspect
that you need an expensive nxn multiplication...
Perhaps in the worst case, but in almost all cases a O(n) convolution
like product should suffice. (I'd like to see an example of the latter,
actually. This relates to the development of digits in a fraction; our
dividend is quite special.)
Now, I am not suggesting that we should necessarily correct the result
at any O(n) price. If we get the result right almost always, that's
sufficient. I believe we can get it right almost always by just
including an "extra" D limb during the Newton iterations, or perhaps
including an extra virtual limb at the end of D, or making sure the
penultimate approximation has at least ceil((n+1)/2) limbs.
> Getting mpn_invert to compute a slightly better value should not be
> hard. For example, including another limb from D (called A in MCA)
> should help. Also, one may avoid doubling the number of limbs in each
I did some tests, the "try mpn_inverse" test passed without any correction
with a code using "3/4" instead of "1/2" for the last iteration. Problems:
1) this does not mean that the code was correct for _any_ input (some
simulations I did with gp-pari suggest it was not)
2) that code was much _slower_, and this is obvious: the last step
required a (n*3n/4) and a (3n/4*3n/4) product instead of (n*n/2) and
Not a very good strategy if you ask me...
> iteration. The basecase will surely be rather large, judging from
> mpn_binvert it will be over 500 limbs, which means an extra limb will
> imply a very little relative cost.
I did not collect timings very carefully, but my first attempt suggest a
smaller threshold, somewhere around 20 limbs, compared to the mpn_invert
currently available in the repository.
That would be *extremely* impressive. You should swing your wand over
mpn_binvert next. :-)
More information about the gmp-devel