mpn_sqrtrem1

[Torbjörn Granlund]

paul zimmermann <Paul.Zimmermann at inria.fr> writes:

  in my current effort of optimizing MPFR for small precisions (1-2 limbs),
  I am currently reviewing mpn_sqrtrem1() and mpn_sqrtrem2(), since we need
  such functions to speed up mpfr_sqrt.
  
  In mpn_sqrtrem1, I wonder what the 0x30000 constant is used for. By exhaustive
  search on all possible values of a1 (there are only 3/4*2^33) on a 64-bit computer,
  the maximal error on the 16-bit approximation x0 is the same without this constant,
  more precisely -0.000023 <= x0/2^32 - (a0/2^64)^(-1/2) <= 0.
  
  The worst error is even smaller (in absolute value) without the 0x30000 constant,
  since with - 0x30000 it is attained for a1=2248147012 (error -0.0000231867), and
  without 0x30000 it is attained for a1=2248147013 (error -0.0000231863).
  
  Also "make check" passes without the 0x30000 constant. If the 0x30000 constant is
  really needed, it would be nice to add an example that breaks make check if we
  remove it.
  
I am not proud of the mpn_sqrtrem1 code.  It was made with some analysis
and lots of testing.  The many undocumented magic constants are ugly.
Perhaps I should never have checked in this code.

Furthermore, for a real performance win one would need to rewrite
mpn_sqrtrem2.

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