overby.me
this repo has no description
1/*
2 * Copyright (c) 2002 Apple Computer, Inc. All rights reserved.
3 *
4 * @APPLE_LICENSE_HEADER_START@
5 *
6 * The contents of this file constitute Original Code as defined in and
7 * are subject to the Apple Public Source License Version 1.1 (the
8 * "License"). You may not use this file except in compliance with the
9 * License. Please obtain a copy of the License at
10 * http://www.apple.com/publicsource and read it before using this file.
11 *
12 * This Original Code and all software distributed under the License are
13 * distributed on an "AS IS" basis, WITHOUT WARRANTY OF ANY KIND, EITHER
14 * EXPRESS OR IMPLIED, AND APPLE HEREBY DISCLAIMS ALL SUCH WARRANTIES,
15 * INCLUDING WITHOUT LIMITATION, ANY WARRANTIES OF MERCHANTABILITY,
16 * FITNESS FOR A PARTICULAR PURPOSE OR NON-INFRINGEMENT. Please see the
17 * License for the specific language governing rights and limitations
18 * under the License.
19 *
20 * @APPLE_LICENSE_HEADER_END@
21 */
22/*******************************************************************************
23* *
24* File logb.c, *
25* Functions logb. *
26* Implementation of IEEE-754 logb for the PowerPC platforms. *
27* *
28* Copyright � 1991-2001 Apple Computer, Inc. All rights reserved. *
29* *
30* Written by A. Sazegari, started on June 1991. *
31* Modified and ported by Robert A. Murley (ram) for Mac OS X. *
32* *
33* A MathLib v4 file. *
34* *
35* Change History (most recent last): *
36* *
37* August 26 1991: removed CFront Version 1.1d17 warnings. *
38* August 27 1991: no errors reported by the test suite. *
39* November 11 1991: changed CLASSEXTENDED to the macro CLASSIFY and *
40* + or - infinity to constants. *
41* November 18 1991: changed the macro CLASSIFY to CLASSEXTENDEDint to *
42* improve performance. *
43* February 07 1992: changed bit operations to macros ( object size is *
44* unchanged ). *
45* September 24 1992: took the "#include support.h" out. *
46* December 03 1992: first rs/6000 port. *
47* August 30 1992: set the divide by zero for the zero argument case. *
48* October 05 1993: corrected the environment. *
49* October 17 1994: replaced all environmental functions with __setflm. *
50* May 28 1997: made speed improvements. *
51* April 30 2001: first mac os x port using gcc. *
52* July 16 2001: replaced __setflm with FEGETENVD/FESETENVD. *
53* August 28 2001: added description of logb function. *
54* September 06 2001: added #if __ppc__. *
55* September 09 2001: added more comments. *
56* September 10 2001: added macros to detect PowerPC and correct compiler. *
57* October 08 2001: removed <CoreServices/CoreServices.h>. *
58* changed compiler errors to warnings. *
59* November 06 2001: commented out warning about Intel architectures. *
60* *
61* W A R N I N G: *
62* These routines require a 64-bit double precision IEEE-754 model. *
63* They are written for PowerPC only and are expecting the compiler *
64* to generate the correct sequence of multiply-add fused instructions. *
65* *
66* These routines are not intended for 32-bit Intel architectures. *
67* *
68* A version of gcc higher than 932 is required. *
69* *
70* GCC compiler options: *
71* optimization level 3 (-O3) *
72* -fschedule-insns -finline-functions -funroll-all-loops *
73* *
74********************************************************************************
75* The C math library offers a similar function called "frexp". It is *
76* different in details from logb, but similar in spirit. This current *
77* implementation of logb follows the recommendation in IEEE Standard 854 *
78* which is different in its handling of denormalized numbers from the IEEE *
79* Standard 754. *
80*******************************************************************************/
81
82#include "math.h"
83#include "fenv.h"
84#include "fp_private.h"
85#include "fenv_private.h"
86
87static const double twoTo52 = 0x1.0p+52; // 4.50359962737049600e15;
88static const double klTod = 4503601774854144.0; // 0x1.000008p52
89static const hexdouble minusInf = HEXDOUBLE(0xfff00000, 0x00000000);
90
91static const float twoTo23 = 0x1.0p+23f; // 8388608.0e0;
92static const hexsingle minusInff = { 0xff800000 };
93
94/*******************************************************************************
95* L O G B *
96********************************************************************************
97* *
98* logb extracts the exponent of its argument, as a signed integral *
99* value. A subnormal argument is treated as though it were first *
100* normalized. Thus *
101* 1 <= x * 2^( - Logb ( x ) ) < 2 *
102*******************************************************************************/
103
104double logb ( double x )
105{
106 hexdouble xInHex;
107 int32_t shiftedExp;
108
109 xInHex.d = x;
110 __NOOP;
111 __NOOP;
112 __NOOP;
113 shiftedExp = ( xInHex.i.hi & 0x7ff00000 ) >> 20;
114
115 if (unlikely( shiftedExp == 2047 ))
116 { // NaN or INF
117 if ( ( ( xInHex.i.hi & 0x80000000 ) == 0 ) || ( x != x ) )
118 return x; // NaN or +INF return x
119 else
120 return -x; // -INF returns +INF
121 }
122
123 if (likely( shiftedExp != 0 )) // normal number
124 shiftedExp -= 1023; // unbias exponent
125 else if ( x == 0.0 )
126 { // zero
127 hexdouble OldEnvironment;
128 FEGETENVD_GRP( OldEnvironment.d ); // raise zero divide for DOMAIN error
129 OldEnvironment.i.lo |= FE_DIVBYZERO;
130 FESETENVD_GRP( OldEnvironment.d );
131 return ( minusInf.d ); // return -infinity
132 }
133 else
134 { // subnormal number
135 xInHex.d *= twoTo52; // scale up
136 __NOOP;
137 __NOOP;
138 __NOOP;
139 shiftedExp = ( xInHex.i.hi & 0x7ff00000 ) >> 20;
140 shiftedExp -= 1075; // unbias exponent
141 }
142
143 if (unlikely( shiftedExp == 0 )) // zero result
144 return ( 0.0 );
145 else
146 { // nonzero result
147 xInHex.d = klTod;
148 __NOOP;
149 __NOOP;
150 __NOOP;
151 xInHex.i.lo += shiftedExp;
152 return ( xInHex.d - klTod );
153 }
154}
155
156int ilogb ( double x )
157{
158 hexdouble xInHex;
159 int32_t shiftedExp;
160
161 xInHex.d = x;
162 __NOOP;
163 __NOOP;
164 __NOOP;
165 shiftedExp = ( xInHex.i.hi & 0x7ff00000 ) >> 20;
166
167 if (unlikely( shiftedExp == 2047 ))
168 { // NaN or INF
169 if (x != x)
170 return FP_ILOGBNAN;
171 else
172 {
173 feraiseexcept( FE_INVALID ); // raise invalid for DOMAIN error
174 return INT32_MAX;
175 }
176 }
177
178 if (likely( shiftedExp != 0 )) // normal number
179 shiftedExp -= 1023; // unbias exponent
180 else if ( x == 0.0 )
181 { // zero
182 feraiseexcept( FE_INVALID ); // raise invalid for DOMAIN error
183 return FP_ILOGB0; // return -infinity
184 }
185 else
186 { // subnormal number
187 xInHex.d *= twoTo52; // scale up
188 __NOOP;
189 __NOOP;
190 __NOOP;
191 shiftedExp = ( xInHex.i.hi & 0x7ff00000 ) >> 20;
192 shiftedExp -= 1075; // unbias exponent
193 }
194
195 return shiftedExp;
196}
197
198float logbf ( float x )
199{
200 hexsingle xInHex;
201 int32_t shiftedExp;
202
203 xInHex.fval = x;
204 __NOOP;
205 __NOOP;
206 __NOOP;
207 shiftedExp = ( xInHex.lval & 0x7f800000 ) >> 23;
208
209 if (unlikely( shiftedExp == 255 ))
210 { // NaN or INF
211 if ( ( ( xInHex.lval & 0x80000000 ) == 0 ) || ( x != x ) )
212 return x; // NaN or +INF return x
213 else
214 return -x; // -INF returns +INF
215 }
216
217 if (likely( shiftedExp != 0 )) // normal number
218 shiftedExp -= 127; // unbias exponent
219 else if ( x == 0.0 )
220 { // zero
221 hexdouble OldEnvironment;
222 FEGETENVD_GRP( OldEnvironment.d ); // raise zero divide for DOMAIN error
223 OldEnvironment.i.lo |= FE_DIVBYZERO;
224 FESETENVD_GRP( OldEnvironment.d );
225 return ( minusInff.fval ); // return -infinity
226 }
227 else
228 { // subnormal number
229 xInHex.fval *= twoTo23; // scale up
230 __NOOP;
231 __NOOP;
232 __NOOP;
233 shiftedExp = ( xInHex.lval & 0x7f800000 ) >> 23;
234 shiftedExp -= 150; // unbias exponent
235 }
236
237 return (float)shiftedExp;
238}
239
240int ilogbf ( float x )
241{
242 hexsingle xInHex;
243 int32_t shiftedExp;
244
245 xInHex.fval = x;
246 __NOOP;
247 __NOOP;
248 __NOOP;
249 shiftedExp = ( xInHex.lval & 0x7f800000 ) >> 23;
250
251 if (unlikely( shiftedExp == 255 ))
252 { // NaN or INF
253 if (x != x)
254 return FP_ILOGBNAN;
255 else
256 {
257 feraiseexcept( FE_INVALID ); // raise invalid for DOMAIN error
258 return INT32_MAX;
259 }
260 }
261
262 if (likely( shiftedExp != 0 )) // normal number
263 shiftedExp -= 127; // unbias exponent
264 else if ( x == 0.0 )
265 { // zero
266 feraiseexcept( FE_INVALID ); // raise invalid for DOMAIN error
267 return FP_ILOGB0; // return -infinity
268 }
269 else
270 { // subnormal number
271 xInHex.fval *= twoTo23; // scale up
272 __NOOP;
273 __NOOP;
274 __NOOP;
275 shiftedExp = ( xInHex.lval & 0x7f800000 ) >> 23;
276 shiftedExp -= 150; // unbias exponent
277 }
278
279 return shiftedExp;
280}