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: scalb.c *
25* *
26* Contains: C source code for implementation of the IEEE-754 scalb *
27* function for double format on PowerPC platforms. *
28* *
29* Copyright � 1992-2001 Apple Computer, Inc. All rights reserved. *
30* *
31* Written by Jon Okada, started on December 1992. *
32* Modified by A. Sazegari (ali) for MathLib v3. *
33* Modified and ported by Robert A. Murley (ram) for Mac OS X. *
34* *
35* A MathLib v4 file. *
36* *
37* Change History ( most recent first ): *
38* *
39* 06 Nov 01 ram commented out warning about Intel architectures. *
40* 10 Oct 01 ram changed compiler errors to warnings. *
41* 10 Sep 01 ali added macros to detect PowerPC and correct compiler. *
42* 09 Sep 01 ali added more comments. *
43* 28 Aug 01 ram added #ifdef __ppc__. *
44* 16 Jul 01 ram replaced DblInHex typedef with hexdouble. *
45* 28 May 97 ali made a speed improvement for large n, *
46* removed scalbl. *
47* 12 Dec 92 JPO First created. *
48* *
49* W A R N I N G: *
50* These routines require a 64-bit double precision IEEE-754 model. *
51* They are written for PowerPC only and are expecting the compiler *
52* to generate the correct sequence of multiply-add fused instructions. *
53* *
54* These routines are not intended for 32-bit Intel architectures. *
55* *
56* A version of gcc higher than 932 is required. *
57* *
58* GCC compiler options: *
59* optimization level 3 (-O3) *
60* -fschedule-insns -finline-functions -funroll-all-loops *
61* *
62*******************************************************************************/
63
64// Put definition of __DARWIN_ALIAS() in sys/cdefs.h in scope
65#define __DARWIN_UNIX03 1
66#include "sys/cdefs.h"
67#include "math.h"
68
69#include "fp_private.h"
70
71static const double twoTo1023 = 0x1.0p+1023;
72static const double twoToM1022 = 0x1.0p-1022;
73static const double twoTo127 = 0x1.0p+127;
74static const double twoToM126 = 0x1.0p-126;
75
76/***********************************************************************
77 Function scalbn
78 Returns its argument x scaled by the factor 2^m. NaNs, signed
79 zeros, and infinities are propagated by this function regardless
80 of the value of n.
81
82 Exceptions: OVERFLOW/INEXACT or UNDERFLOW/INEXACT may occur;
83 INVALID for signaling NaN inputs (quiet NaN returned).
84***********************************************************************/
85
86double scalbn ( double x, int n )
87{
88 hexdouble xInHex;
89
90 xInHex.i.lo = 0u; // init. low half of xInHex
91
92 if ( n > 1023 )
93 { // large positive scaling
94 if ( n > 2097 ) // huge scaling
95 {
96 register volatile double s, t, u;
97 s = x * twoTo1023;
98 t = s * twoTo1023;
99 u = t * twoTo1023;
100 return u;
101 }
102 while ( n > 1023 )
103 { // scale reduction loop
104 x *= twoTo1023; // scale x by 2^1023
105 n -= 1023; // reduce n by 1023
106 }
107 }
108
109 else if ( n < -1022 )
110 { // large negative scaling
111 if ( n < -2098 ) // huge negative scaling
112 {
113 register volatile double s, t, u;
114 s = x * twoToM1022;
115 t = s * twoToM1022;
116 u = t * twoToM1022;
117 return u;
118 }
119 while ( n < -1022 )
120 { // scale reduction loop
121 x *= twoToM1022; // scale x by 2^( -1022 )
122 n += 1022; // incr n by 1022
123 }
124 }
125
126 if ( -127 < n && n < 128 )
127 {
128/*******************************************************************************
129* -126 <= n <= -127; convert n to single scale factor. *
130* Allows a store-forward to execute successfully *
131*******************************************************************************/
132 hexsingle XInHex;
133
134 XInHex.lval = ( ( uint32_t ) ( n + 127 ) ) << 23;
135
136 __NOOP;
137 __NOOP;
138 __NOOP;
139 return ( x * XInHex.fval );
140 }
141
142/*******************************************************************************
143* -1022 <= n <= 1023; convert n to double scale factor. *
144*******************************************************************************/
145
146 xInHex.i.hi = ( ( uint32_t ) ( n + 1023 ) ) << 20;
147 return ( x * xInHex.d );
148}
149
150double scalbln ( double x, long int n )
151{
152 int m;
153
154 // Clip n
155 if (unlikely(n > 2097))
156 m = 2098;
157 else if (unlikely(n < -2098))
158 m = -2099;
159 else
160 m = (int) n;
161
162 return scalbn(x, m);
163}
164
165float scalbnf ( float x, int n )
166{
167 hexsingle xInHex;
168
169 if ( n > 127 )
170 { // large positive scaling
171 if ( n > 276 ) // huge scaling
172 {
173 register volatile float s, t, u;
174 s = x * (float) twoTo127;
175 t = s * (float) twoTo127;
176 u = t * (float) twoTo127;
177 return u;
178 }
179 while ( n > 127 )
180 { // scale reduction loop
181 x *= (float) twoTo127; // scale x by 2^127
182 n -= 127; // reduce n by 127
183 }
184 }
185
186 else if ( n < -126 )
187 { // large negative scaling
188 if ( n < -277 ) // huge negative scaling
189 {
190 register volatile float s, t, u;
191 s = x * (float) twoToM126;
192 t = s * (float) twoToM126;
193 u = t * (float) twoToM126;
194 return u;
195 }
196 while ( n < -126 )
197 { // scale reduction loop
198 x *= (float) twoToM126; // scale x by 2^( -126 )
199 n += 126; // incr n by 126
200 }
201 }
202
203/*******************************************************************************
204* -126 <= n <= 127; convert n to float scale factor. *
205*******************************************************************************/
206
207 xInHex.lval = ( ( uint32_t ) ( n + 127 ) ) << 23;
208
209 // Force the fetch for xInHex.fval to the next cycle to avoid Store/Load hazard.
210 __NOOP;
211 __NOOP;
212 __NOOP;
213
214 return ( x * xInHex.fval );
215}
216
217float scalblnf ( float x, long int n )
218{
219 int m;
220
221 // Clip n
222 if (unlikely(n > 276))
223 m = 277;
224 else if (unlikely(n < -277))
225 m = -278;
226 else
227 m = (int) n;
228
229 return scalbnf(x, m);
230}
231
232// POSIX mandated signature for "scalb"
233extern double scalb ( double, double ) __DARWIN_ALIAS(scalb);
234double scalb ( double x, double n )
235{
236 int m;
237
238 if ( n > 2098.0 )
239 m = 2098;
240 else if ( n < -2099.0 )
241 m = -2099;
242 else m = (int) n;
243
244 return scalbn( x, m );
245}