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 28long double 29ldexpl (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 39long double 40ldexpl (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 81int 82main (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