* mpf/sqrt.c: Change tsize calculation to get prec limbs result

	always, previously got prec+1 when exp was odd.

1 files changed, 29 insertions, 2 deletions

diff --git a/mpf/sqrt.c b/mpf/sqrt.c
index aea67204f..ccd5019f7 100644
--- a/mpf/sqrt.c
+++ b/mpf/sqrt.c
@@ -1,6 +1,6 @@
 /* mpf_sqrt -- Compute the square root of a float.
 
-Copyright 1993, 1994, 1996, 2000, 2001 Free Software Foundation, Inc.
+Copyright 1993, 1994, 1996, 2000, 2001, 2004 Free Software Foundation, Inc.
 
 This file is part of the GNU MP Library.
 
@@ -23,6 +23,33 @@ MA 02111-1307, USA. */
 #include "gmp.h"
 #include "gmp-impl.h"
 
+
+/* As usual, the aim is to produce PREC(r) limbs of result, with the high
+   limb non-zero.  This is accomplished by applying mpn_sqrtrem to either
+   2*prec or 2*prec-1 limbs, both such sizes resulting in prec limbs.
+
+   The choice between 2*prec or 2*prec-1 limbs is based on the input
+   exponent.  With b=2^GMP_NUMB_BITS is the limb base then we can think of
+   effectively taking out a factor b^(2k), for suitable k, to get to an
+   integer input of the desired size ready for mpn_sqrtrem.  It must be an
+   even power taken out, ie. an even number of limbs, so the square root
+   gives b^k, keeping the radix point on a limb boundary.  So if EXP(r) is
+   even we'll get an even number of input limbs (ie. 2*prec), or if EXP(r)
+   is odd we get an odd number (ie. 2*prec-1).
+
+   Further limbs below the 2*prec or 2*prec-1 used don't affect the result
+   and are simply truncated.  This can be seen by considering an integer x,
+   with s=floor(sqrt(x)).  s is the unique integer satisfying s^2 <= x <
+   (s+1)^2.  Notice that adding a fraction part to x (ie. some further bits)
+   doesn't change the inequality, s remains the unique solution.  Working
+   suitable factors of 2 into this argument lets it apply to an intended
+   precision at any position for any x, not just the integer binary point.
+
+   If the input is smaller than 2*prec or 2*prec-1, then we just pad with
+   zeros, that of course being our usual interpretation of short inputs.
+   The effect is to extend the root beyond the size of the input (for
+   instance into fractional limbs if u is an integer).  */
+
 void
 mpf_sqrt (mpf_ptr r, mpf_srcptr u)
 {
@@ -46,7 +73,7 @@ mpf_sqrt (mpf_ptr r, mpf_srcptr u)
 
   prec = r->_mp_prec;
   rexp = (u->_mp_exp + 1) >> 1;	/* round towards -inf */
-  tsize = 2 * prec + (u->_mp_exp & 1);
+  tsize = 2 * prec - (u->_mp_exp & 1);
 
   up = u->_mp_d;
   tp = (mp_ptr) TMP_ALLOC (tsize * BYTES_PER_MP_LIMB);