Anomaly in mpn_sqrtrem and mpn_rottrem

[Torbjörn Granlund]

nisse at lysator.liu.se (Niels Möller) writes:

  tg at gmplib.org (Torbjörn Granlund) writes:
  
  > Surely, we could make mpn_rootrem run faster, in particular for small
  > arguments.  But also for large arguments, 2x slowdown seems like a lot.
  
  I've had a quick look. Both mpn_dc_sqrtrem and mpn_rootrem_internal
  (which seem to be the work horses for larger operands) do a division in
  the loop. The latter is organized as a newton iteration, the former
  isn't (but I guess the computation is more or less equivalent to a
  Newton iteration).
  
  The code is a bit too complex for me to really compare them...
  
Ahum, same sensation here...

  For larger operands, I imagine a rewrite to use a division-free
  Newton-iteration for an approximation of the inverse root would be
  faster.
  
I played some years ago with a loop computing A^(-0.5), which naturally
becomes division-free.  An alternative would be to keep the current
algorithm, but use the fact that the previous divisor is a damn good
staring approximation to the new one, and use Newton to maintain an
inverse to it.

For mpn_rootrem, i.e., a^(1/k), I suppose we could make the work take
O(M(log(a^(1/k)))) as opposed to the current O(M(log(a))), where M(n) is
the time for an n-bit multiply.  We'd need to make use of
mpn_pow_1_highpart.  This will be a lot of work.



-- 
Torbjörn
Please encrypt, key id 0xC8601622