Updated from GMP 1.906.7
[kopensolaris-gnu/glibc.git] / sysdeps / ieee754 / ldexp.c
1 /* Copyright (C) 1992, 1995 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 <errno.h>
42 #include "ieee754.h"
43
44 double
45 DEFUN(ldexp, (x, exp),
46       double x AND int exp)
47 {
48   union ieee754_double u;
49   unsigned int exponent;
50
51   u.d = x;
52 #define x u.d
53
54   exponent = u.ieee.exponent;
55
56   /* The order of the tests is carefully chosen to handle
57      the usual case first, with no branches taken.  */
58
59   if (exponent != 0)
60     {
61       /* X is nonzero and not denormalized.  */
62
63       if (exponent <= DBL_MAX_EXP - DBL_MIN_EXP + 1)
64         {
65           /* X is finite.  When EXP < 0, overflow is actually underflow.  */
66
67           exponent += exp;
68
69           if (exponent != 0)
70             {
71               if (exponent <= DBL_MAX_EXP - DBL_MIN_EXP + 1)
72                 {
73                   /* In range.  */
74                   u.ieee.exponent = exponent;
75                   return x;
76                 }
77
78               if (exp >= 0)
79               overflow:
80                 {
81                   CONST int negative = u.ieee.negative;
82                   u.d = HUGE_VAL;
83                   u.ieee.negative = negative;
84                   errno = ERANGE;
85                   return u.d;
86                 }
87
88               if (exponent <= - (unsigned int) (DBL_MANT_DIG + 1))
89                 {
90                   /* Underflow.  */
91                   CONST int negative = u.ieee.negative;
92                   u.d = 0.0;
93                   u.ieee.negative = negative;
94                   errno = ERANGE;
95                   return u.d;
96                 }
97             }
98
99           /* Gradual underflow.  */
100           u.ieee.exponent = 1;
101           u.d *= ldexp (1.0, (int) exponent - 1);
102           if (u.ieee.mantissa0 == 0 && u.ieee.mantissa1 == 0)
103             /* Underflow.  */
104             errno = ERANGE;
105           return u.d;
106         }
107
108       /* X is +-infinity or NaN.  */
109       if (u.ieee.mantissa0 == 0 && u.ieee.mantissa1 == 0)
110         {
111           /* X is +-infinity.  */
112           if (exp >= 0)
113             goto overflow;
114           else
115             {
116               /* (infinity * number < 1).  With infinite precision,
117                  (infinity / finite) would be infinity, but otherwise it's
118                  safest to regard (infinity / 2) as indeterminate.  The
119                  infinity might be (2 * finite).  */
120               CONST int negative = u.ieee.negative;
121               u.d = NAN;
122               u.ieee.negative = negative;
123               errno = EDOM;
124               return u.d;
125             }
126         }
127
128       /* X is NaN.  */
129       errno = EDOM;
130       return u.d;
131     }
132
133   /* X is zero or denormalized.  */
134   if (u.ieee.mantissa0 == 0 && u.ieee.mantissa1 == 0)
135     /* X is +-0.0. */
136     return x;
137
138   /* X is denormalized.
139      Multiplying by 2 ** DBL_MANT_DIG normalizes it;
140      we then subtract the DBL_MANT_DIG we added to the exponent.  */
141   return ldexp (x * ldexp (1.0, DBL_MANT_DIG), exp - DBL_MANT_DIG);
142 }
143
144 /* Compatibility names for the same function.  */
145 weak_alias (ldexp, __scalb)
146 weak_alias (ldexp, scalb)