1/* mpf_sqrt_ui -- Compute the square root of an unsigned integer. 2 3Copyright 1993, 1994, 1996, 2000, 2001, 2004, 2005, 2015 Free Software 4Foundation, Inc. 5 6This file is part of the GNU MP Library. 7 8The GNU MP Library is free software; you can redistribute it and/or modify 9it 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 15or 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 21or both in parallel, as here. 22 23The GNU MP Library is distributed in the hope that it will be useful, but 24WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY 25or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License 26for more details. 27 28You should have received copies of the GNU General Public License and the 29GNU Lesser General Public License along with the GNU MP Library. If not, 30see https://www.gnu.org/licenses/. */ 31 32#include <stdio.h> /* for NULL */ 33#include "gmp-impl.h" 34 35 36/* As usual the aim is to produce PREC(r) limbs of result with the high limb 37 non-zero. That high limb will end up floor(sqrt(u)), and limbs below are 38 produced by padding the input with zeros, two for each desired result 39 limb, being 2*(prec-1) for a total 2*prec-1 limbs passed to mpn_sqrtrem. 40 The way mpn_sqrtrem calculates floor(sqrt(x)) ensures the root is correct 41 to the intended accuracy, ie. truncated to prec limbs. 42 43 With nails, u might be two limbs, in which case a total 2*prec limbs is 44 passed to mpn_sqrtrem (still giving a prec limb result). If uhigh is 45 zero we adjust back to 2*prec-1, since mpn_sqrtrem requires the high 46 non-zero. 2*prec limbs are always allocated, even when uhigh is zero, so 47 the store of uhigh can be done without a conditional. 48 49 u==0 is a special case so the rest of the code can assume the result is 50 non-zero (ie. will have a non-zero high limb on the result). 51 52 Not done: 53 54 No attempt is made to identify perfect squares. It's considered this can 55 be left to an application if it might occur with any frequency. As it 56 stands, mpn_sqrtrem does its normal amount of work on a perfect square 57 followed by zero limbs, though of course only an mpn_sqrtrem1 would be 58 actually needed. We also end up leaving our mpf result with lots of low 59 trailing zeros, slowing down subsequent operations. 60 61 We're not aware of any optimizations that can be made using the fact the 62 input has lots of trailing zeros (apart from the perfect square 63 case). */ 64 65 66/* 1 if we (might) need two limbs for u */ 67#define U2 (GMP_NUMB_BITS < BITS_PER_ULONG) 68 69void 70mpf_sqrt_ui (mpf_ptr r, unsigned long int u) 71{ 72 mp_size_t rsize, zeros; 73 mp_ptr tp; 74 mp_size_t prec; 75 TMP_DECL; 76 77 if (UNLIKELY (u <= 1)) 78 { 79 SIZ (r) = EXP (r) = u; 80 *PTR (r) = u; 81 return; 82 } 83 84 TMP_MARK; 85 86 prec = PREC (r); 87 zeros = 2 * prec - 2; 88 rsize = zeros + 1 + U2; 89 90 tp = TMP_ALLOC_LIMBS (rsize); 91 92 MPN_ZERO (tp, zeros); 93 tp[zeros] = u & GMP_NUMB_MASK; 94 95#if U2 96 { 97 mp_limb_t uhigh = u >> GMP_NUMB_BITS; 98 tp[zeros + 1] = uhigh; 99 rsize -= (uhigh == 0); 100 } 101#endif 102 103 mpn_sqrtrem (PTR (r), NULL, tp, rsize); 104 105 SIZ (r) = prec; 106 EXP (r) = 1; 107 TMP_FREE; 108} 109