[patrick.pelissier@gmail.com: Fwd: Optimizing 1/x]

[Ashod Nakashian]

Thanks Patrick Pe'lissier!

Paul Zimmermann wrote:
> Patrick Pe'lissier forwarded me the following mail (I'm not on gmp-discuss).
> Two remarks:
> 
> * what GMP needs is a fast inverse at the mpn level, which would compute,
>   given a n-limb normalized input d, a n-limb quotient q (plus possibly one
>   carry bit), and a n-limb remainder 0 <= r < d, such that
> 
>   2^(2*n*mp_bits_per_limb) = q * d + r.
> 
>   With such a function, computing 1/x in mpf would become easy.
> 

I agree; an mpn-level implementation would probably be the way to go 
(vs. a higher-level).

> * Newton's iteration isn't currently implemented in GMP, even for the
>   integer division. The main reason is that you need to carefully analyze
>   the roundoff errors to get a correct answer (as in the above proposal).
>   I've implemented several variants which are quite fast; however the
>   fastest known algorithm is currently under a patent from HP, and I failed
>   so far to obtain a license for use within GMP...

Well, I assume then that the benefits of such a method does *at least* 
pay-off. The next question would be how slower is the next fastest 
non-patented algorithm compared to HP's? Also, although I'm anything but 
a lawyer, I'd be surprised if law forbids using a patented algorithm for 
educational/non-profit purposes. So I assume that one could argue that 
making patch (with the said HP algorithm) available would benefit many 
and the burden of following the patent laws would be on the user, not 
the provider of the code.

Have you considered making contributions with your implementation to GMP 
or seperate from GMP in form of patches to GMP?

Out of curiosity, how much faster is the HP algorithm compared to the 
current GMP division implementation. I assume it's sub linear (or is it 
not?), but in practice is it close to log2(n) + C? Its complexity is 
probably a function in the multiplication algorithm's complexity.

-Ash

> 
> Paul Zimmermann
> 
> - ---------- Forwarded message ----------
> From: Ashod Nakashian <saghmos at xter.net>
> Date: Sep 11, 2005 11:18 AM
> Subject: Optimizing 1/x
> To: gmp-discuss at swox.com
> 
> 
> Hi,
> 
> Since GMP is highly optimized for the most common (and sometimes not so
> common) calculations, I assumed that there would be an optimized
> function for calculating the multiplicative inverse of a real (1/'mpf_t').
> 
> I don't seem to find such a function, so I assumed that mpf_ui_div might
> have a specialized path for values of 1 (i.e. inversion), but there
> seems to be no such path.
> 
> As I understand it, the code currently creates an mpf_t from the 'ui'
> value. This value is normalized such that the division is done using the
> mpf divisor as is and corrected by the exponent of the divisor.
> 
> This seems like a very simplistic approach that is sound (and probably
> optimal) for arbitrary numerators. But the dire question is:
> 
> Why not use the Newton Iterations to computer 1/x?
> 
> And if it's not worth it (i.e. it won't be faster), why not simply
> optimize the current algorithm using the fact that the numerator is a
> simple 1.
> 
> It is hard to imagine that the GMP user base (in deed the developers as
> well) didn't come across the need of a faster inversion calculation.
> 
> Am I missing something? Is it assumed since Newton's iterative method is
> fairly easy to implement, one will do just that? (That still doesn't
> sound like th GMP philosophy to me.)
> 
> So what do/did you guys do when you need(ed) to find 1/x quite rapidly?
> 
> - -Ash
> _______________________________________________
> 
>