Formerly ieee754/ldexp.c.~3~
[kopensolaris-gnu/glibc.git] / sysdeps / ieee754 / ldexp.c
1 /* Copyright (C) 1992 Free Software Foundation, Inc.
2 This file is part of the GNU C Library.
3
4 The GNU C Library is free software; you can redistribute it and/or
5 modify it under the terms of the GNU Library General Public License as
6 published by the Free Software Foundation; either version 2 of the
7 License, or (at your option) any later version.
8
9 The GNU C Library is distributed in the hope that it will be useful,
10 but WITHOUT ANY WARRANTY; without even the implied warranty of
11 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
12 Library General Public License for more details.
13
14 You should have received a copy of the GNU Library General Public
15 License along with the GNU C Library; see the file COPYING.LIB.  If
16 not, write to the Free Software Foundation, Inc., 675 Mass Ave,
17 Cambridge, MA 02139, USA.  */
18
19 /*
20  * Copyright (c) 1985 Regents of the University of California.
21  * All rights reserved.
22  *
23  * Redistribution and use in source and binary forms are permitted provided
24  * that: (1) source distributions retain this entire copyright notice and
25  * comment, and (2) distributions including binaries display the following
26  * acknowledgement:  ``This product includes software developed by the
27  * University of California, Berkeley and its contributors'' in the
28  * documentation or other materials provided with the distribution and in
29  * all advertising materials mentioning features or use of this software.
30  * Neither the name of the University nor the names of its contributors may
31  * be used to endorse or promote products derived from this software without
32  * specific prior written permission.
33  * THIS SOFTWARE IS PROVIDED ``AS IS'' AND WITHOUT ANY EXPRESS OR IMPLIED
34  * WARRANTIES, INCLUDING, WITHOUT LIMITATION, THE IMPLIED WARRANTIES OF
35  * MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
36  */
37
38 #include <ansidecl.h>
39 #include <math.h>
40 #include <float.h>
41 #include "ieee754.h"
42
43 double
44 DEFUN(ldexp, (x, exp),
45       double x AND int exp)
46 {
47   union ieee754_double u;
48   unsigned int exponent;
49
50   u.d = x;
51 #define x u.d
52
53   exponent = u.ieee.exponent;
54
55   /* The order of the tests is carefully chosen to handle
56      the usual case first, with no branches taken.  */
57
58   if (exponent != 0)
59     {
60       /* X is nonzero and not denormalized.  */
61
62       if (exponent <= DBL_MAX_EXP - DBL_MIN_EXP + 1)
63         {
64           /* X is finite.  When EXP < 0, overflow is actually underflow.  */
65
66           exponent += exp;
67
68           if (exponent != 0)
69             {
70               if (exponent <= DBL_MAX_EXP - DBL_MIN_EXP + 1)
71                 {
72                   /* In range.  */
73                   u.ieee.exponent = exponent;
74                   return x;
75                 }
76
77               if (exp >= 0)
78               overflow:
79                 {
80                   CONST int negative = u.ieee.negative;
81                   u.d = HUGE_VAL;
82                   u.ieee.negative = negative;
83                   errno = ERANGE;
84                   return u.d;
85                 }
86
87               if (exponent <= - (unsigned int) (DBL_MANT_DIG + 1))
88                 {
89                   /* Underflow.  */
90                   CONST int negative = u.ieee.negative;
91                   u.d = 0.0;
92                   u.ieee.negative = negative;
93                   errno = ERANGE;
94                   return u.d;
95                 }
96             }
97
98           /* Gradual underflow.  */
99           u.ieee.exponent = 1;
100           u.d *= ldexp (1.0, (int) exponent - 1);
101           if (u.ieee.mantissa0 == 0 && u.ieee.mantissa1 == 0)
102             /* Underflow.  */
103             errno = ERANGE;
104           return u.d;
105         }
106
107       /* X is +-infinity or NaN.  */
108       if (u.ieee.mantissa0 == 0 && u.ieee.mantissa1 == 0)
109         {
110           /* X is +-infinity.  */
111           if (exp >= 0)
112             goto overflow;
113           else
114             {
115               /* (infinity * number < 1).  With infinite precision,
116                  (infinity / finite) would be infinity, but otherwise it's
117                  safest to regard (infinity / 2) as indeterminate.  The
118                  infinity might be (2 * finite).  */
119               CONST int negative = u.ieee.negative;
120               u.d = NAN;
121               u.ieee.negative = negative;
122               errno = EDOM;
123               return u.d;
124             }
125         }
126
127       /* X is NaN.  */
128       errno = EDOM;
129       return u.d;
130     }
131
132   /* X is zero or denormalized.  */
133   if (u.ieee.mantissa0 == 0 && u.ieee.mantissa1 == 0)
134     /* X is +-0.0. */
135     return x;
136
137   /* X is denormalized.
138      Multiplying by 2 ** DBL_MANT_DIG normalizes it;
139      we then subtract the DBL_MANT_DIG we added to the exponent.  */
140   return ldexp (x * ldexp (1.0, DBL_MANT_DIG), exp - DBL_MANT_DIG);
141 }