1/* __gmpfr_isqrt && __gmpfr_cuberoot -- Integer square root and cube root 2 3Copyright 2004, 2005, 2006, 2007, 2008, 2009, 2010, 2011, 2012, 2013 Free Software Foundation, Inc. 4Contributed by the AriC and Caramel projects, INRIA. 5 6This file is part of the GNU MPFR Library. 7 8The GNU MPFR Library is free software; you can redistribute it and/or modify 9it under the terms of the GNU Lesser General Public License as published by 10the Free Software Foundation; either version 3 of the License, or (at your 11option) any later version. 12 13The GNU MPFR Library is distributed in the hope that it will be useful, but 14WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY 15or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public 16License for more details. 17 18You should have received a copy of the GNU Lesser General Public License 19along with the GNU MPFR Library; see the file COPYING.LESSER. If not, see 20http://www.gnu.org/licenses/ or write to the Free Software Foundation, Inc., 2151 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA. */ 22 23#include "mpfr-impl.h" 24 25/* returns floor(sqrt(n)) */ 26unsigned long 27__gmpfr_isqrt (unsigned long n) 28{ 29 unsigned long i, s; 30 31 /* First find an approximation to floor(sqrt(n)) of the form 2^k. */ 32 i = n; 33 s = 1; 34 while (i >= 2) 35 { 36 i >>= 2; 37 s <<= 1; 38 } 39 40 do 41 { 42 s = (s + n / s) / 2; 43 } 44 while (!(s*s <= n && (s*s > s*(s+2) || n <= s*(s+2)))); 45 /* Short explanation: As mathematically s*(s+2) < 2*ULONG_MAX, 46 the condition s*s > s*(s+2) is evaluated as true when s*(s+2) 47 "overflows" but not s*s. This implies that mathematically, one 48 has s*s <= n <= s*(s+2). If s*s "overflows", this means that n 49 is "large" and the inequality n <= s*(s+2) cannot be satisfied. */ 50 return s; 51} 52 53/* returns floor(n^(1/3)) */ 54unsigned long 55__gmpfr_cuberoot (unsigned long n) 56{ 57 unsigned long i, s; 58 59 /* First find an approximation to floor(cbrt(n)) of the form 2^k. */ 60 i = n; 61 s = 1; 62 while (i >= 4) 63 { 64 i >>= 3; 65 s <<= 1; 66 } 67 68 /* Improve the approximation (this is necessary if n is large, so that 69 mathematically (s+1)*(s+1)*(s+1) isn't much larger than ULONG_MAX). */ 70 if (n >= 256) 71 { 72 s = (2 * s + n / (s * s)) / 3; 73 s = (2 * s + n / (s * s)) / 3; 74 s = (2 * s + n / (s * s)) / 3; 75 } 76 77 do 78 { 79 s = (2 * s + n / (s * s)) / 3; 80 } 81 while (!(s*s*s <= n && (s*s*s > (s+1)*(s+1)*(s+1) || 82 n < (s+1)*(s+1)*(s+1)))); 83 return s; 84} 85