1 /* Emulation for ldexpl. 2 Contributed by Paolo Bonzini 3 4 Copyright 2002-2003, 2007-2012 Free Software Foundation, Inc. 5 6 This file is part of gnulib. 7 8 This program is free software: you can redistribute it and/or modify 9 it under the terms of the GNU General Public License as published by 10 the Free Software Foundation; either version 3 of the License, or 11 (at your option) any later version. 12 13 This program is distributed in the hope that it will be useful, 14 but WITHOUT ANY WARRANTY; without even the implied warranty of 15 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 16 GNU General Public License for more details. 17 18 You should have received a copy of the GNU General Public License 19 along with this program. If not, see <http://www.gnu.org/licenses/>. */ 20 21 #include <config.h> 22 23 /* Specification. */ 24 #include <math.h> 25 26 #if HAVE_SAME_LONG_DOUBLE_AS_DOUBLE 27 28 long double 29 ldexpl (long double x, int exp) 30 { 31 return ldexp (x, exp); 32 } 33 34 #else 35 36 # include <float.h> 37 # include "fpucw.h" 38 39 long double 40 ldexpl (long double x, int exp) 41 { 42 long double factor; 43 int bit; 44 DECL_LONG_DOUBLE_ROUNDING 45 46 BEGIN_LONG_DOUBLE_ROUNDING (); 47 48 /* Check for zero, nan and infinity. */ 49 if (!(isnanl (x) || x + x == x)) 50 { 51 if (exp < 0) 52 { 53 exp = -exp; 54 factor = 0.5L; 55 } 56 else 57 factor = 2.0L; 58 59 if (exp > 0) 60 for (bit = 1;;) 61 { 62 /* Invariant: Here bit = 2^i, factor = 2^-2^i or = 2^2^i, 63 and bit <= exp. */ 64 if (exp & bit) 65 x *= factor; 66 bit <<= 1; 67 if (bit > exp) 68 break; 69 factor = factor * factor; 70 } 71 } 72 73 END_LONG_DOUBLE_ROUNDING (); 74 75 return x; 76 } 77 78 #endif 79 80 #if 0 81 int 82 main (void) 83 { 84 long double x; 85 int y; 86 for (y = 0; y < 29; y++) 87 printf ("%5d %.16Lg %.16Lg\n", y, ldexpl (0.8L, y), ldexpl (0.8L, -y) * ldexpl (0.8L, y)); 88 } 89 #endif 90
