contrib/gmp/mpf/sqrt.c

   1 /* mpf_sqrt -- Compute the square root of a float.
   2
   3 Copyright 1993, 1994, 1996, 2000, 2001, 2004, 2005 Free Software Foundation,
   4 Inc.
   5
   6 This file is part of the GNU MP Library.
   7
   8 The GNU MP Library is free software; you can redistribute it and/or modify
   9 it under the terms of the GNU Lesser General Public License as published by
  10 the Free Software Foundation; either version 3 of the License, or (at your
  11 option) any later version.
  12
  13 The GNU MP Library is distributed in the hope that it will be useful, but
  14 WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
  15 or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
  16 License for more details.
  17
  18 You should have received a copy of the GNU Lesser General Public License
  19 along with the GNU MP Library.  If not, see http://www.gnu.org/licenses/.  */
  20
  21 #include <stdio.h> /* for NULL */
  22 #include "gmp.h"
  23 #include "gmp-impl.h"
  24
  25
  26 /* As usual, the aim is to produce PREC(r) limbs of result, with the high
  27    limb non-zero.  This is accomplished by applying mpn_sqrtrem to either
  28    2*prec or 2*prec-1 limbs, both such sizes resulting in prec limbs.
  29
  30    The choice between 2*prec or 2*prec-1 limbs is based on the input
  31    exponent.  With b=2^GMP_NUMB_BITS the limb base then we can think of
  32    effectively taking out a factor b^(2k), for suitable k, to get to an
  33    integer input of the desired size ready for mpn_sqrtrem.  It must be an
  34    even power taken out, ie. an even number of limbs, so the square root
  35    gives factor b^k and the radix point is still on a limb boundary.  So if
  36    EXP(r) is even we'll get an even number of input limbs 2*prec, or if
  37    EXP(r) is odd we get an odd number 2*prec-1.
  38
  39    Further limbs below the 2*prec or 2*prec-1 used don't affect the result
  40    and are simply truncated.  This can be seen by considering an integer x,
  41    with s=floor(sqrt(x)).  s is the unique integer satisfying s^2 <= x <
  42    (s+1)^2.  Notice that adding a fraction part to x (ie. some further bits)
  43    doesn't change the inequality, s remains the unique solution.  Working
  44    suitable factors of 2 into this argument lets it apply to an intended
  45    precision at any position for any x, not just the integer binary point.
  46
  47    If the input is smaller than 2*prec or 2*prec-1, then we just pad with
  48    zeros, that of course being our usual interpretation of short inputs.
  49    The effect is to extend the root beyond the size of the input (for
  50    instance into fractional limbs if u is an integer).  */
  51
  52 void
  53 mpf_sqrt (mpf_ptr r, mpf_srcptr u)
  54 {
  55   mp_size_t usize;
  56   mp_ptr up, tp;
  57   mp_size_t prec, tsize;
  58   mp_exp_t uexp, expodd;
  59   TMP_DECL;
  60
  61   usize = u->_mp_size;
  62   if (usize <= 0)
  63     {
  64       if (usize < 0)
  65         SQRT_OF_NEGATIVE;
  66       r->_mp_size = 0;
  67       r->_mp_exp = 0;
  68       return;
  69     }
  70
  71   TMP_MARK;
  72
  73   uexp = u->_mp_exp;
  74   prec = r->_mp_prec;
  75   up = u->_mp_d;
  76
  77   expodd = (uexp & 1);
  78   tsize = 2 * prec - expodd;
  79   r->_mp_size = prec;
  80   r->_mp_exp = (uexp + expodd) / 2;    /* ceil(uexp/2) */
  81
  82   /* root size is ceil(tsize/2), this will be our desired "prec" limbs */
  83   ASSERT ((tsize + 1) / 2 == prec);
  84
  85   tp = (mp_ptr) TMP_ALLOC (tsize * BYTES_PER_MP_LIMB);
  86
  87   if (usize > tsize)
  88     {
  89       up += usize - tsize;
  90       usize = tsize;
  91       MPN_COPY (tp, up, tsize);
  92     }
  93   else
  94     {
  95       MPN_ZERO (tp, tsize - usize);
  96       MPN_COPY (tp + (tsize - usize), up, usize);
  97     }
  98
  99   mpn_sqrtrem (r->_mp_d, NULL, tp, tsize);
 100
 101   TMP_FREE;
 102 }