Better support of signalling NaNs.
[pspp] / lib / ldexpl.c
1 /* Emulation for ldexpl.
2    Contributed by Paolo Bonzini
3
4    Copyright 2002, 2003, 2007 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 2, or (at your option)
11    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 along
19    with this program; if not, write to the Free Software Foundation,
20    Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.  */
21
22 #include <config.h>
23
24 /* Specification.  */
25 #include <math.h>
26
27 #include <float.h>
28 #include "isnanl.h"
29
30 long double
31 ldexpl(long double x, int exp)
32 {
33   long double factor;
34   int bit;
35
36   /* Check for zero, nan and infinity. */
37   if (isnanl (x) || x + x == x)
38     return x;
39
40   if (exp < 0)
41     {
42       exp = -exp;
43       factor = 0.5L;
44     }
45   else
46     factor = 2.0L;
47
48   for (bit = 1; bit <= exp; bit <<= 1, factor *= factor)
49     if (exp & bit)
50       x *= factor;
51
52   return x;
53 }
54
55 #if 0
56 int
57 main (void)
58 {
59   long double x;
60   int y;
61   for (y = 0; y < 29; y++)
62     printf ("%5d %.16Lg %.16Lg\n", y, ldexpl(0.8L, y), ldexpl(0.8L, -y) * ldexpl(0.8L, y));
63 }
64 #endif