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