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