1/* mpfr_scale2 -- multiply a double float by 2^exp 2 3Copyright 1999, 2000, 2001, 2002, 2003, 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 <float.h> /* for DBL_EPSILON */ 24#include "mpfr-impl.h" 25 26/* Note: we could use the ldexp function, but since we want not to depend on 27 math.h, we write our own implementation. */ 28 29/* multiplies 1/2 <= d <= 1 by 2^exp */ 30double 31mpfr_scale2 (double d, int exp) 32{ 33#if _GMP_IEEE_FLOATS 34 { 35 union ieee_double_extract x; 36 37 if (MPFR_UNLIKELY (d == 1.0)) 38 { 39 d = 0.5; 40 exp ++; 41 } 42 43 /* now 1/2 <= d < 1 */ 44 45 /* infinities and zeroes have already been checked */ 46 MPFR_ASSERTD (-1073 <= exp && exp <= 1025); 47 48 x.d = d; 49 if (MPFR_UNLIKELY (exp < -1021)) /* subnormal case */ 50 { 51 x.s.exp += exp + 52; 52 x.d *= DBL_EPSILON; 53 } 54 else /* normalized case */ 55 { 56 x.s.exp += exp; 57 } 58 return x.d; 59 } 60#else /* _GMP_IEEE_FLOATS */ 61 { 62 double factor; 63 64 /* An overflow may occurs (example: 0.5*2^1024) */ 65 if (d < 1.0) 66 { 67 d += d; 68 exp--; 69 } 70 /* Now 1.0 <= d < 2.0 */ 71 72 if (exp < 0) 73 { 74 factor = 0.5; 75 exp = -exp; 76 } 77 else 78 { 79 factor = 2.0; 80 } 81 while (exp != 0) 82 { 83 if ((exp & 1) != 0) 84 d *= factor; 85 exp >>= 1; 86 factor *= factor; 87 } 88 return d; 89 } 90#endif 91} 92