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 ceilfloor.c, *
25* Function ceil(x) and floor(x), *
26* Implementation of ceil and floor for the PowerPC. *
27* *
28* Copyright � 1991-2001 Apple Computer, Inc. All rights reserved. *
29* *
30* Written by Ali Sazegari, started on November 1991, *
31* Modified and ported by Robert A. Murley (ram) for Mac OS X. *
32* *
33* A MathLib v4 file. *
34* *
35* December 03 1992: first rs6000 port. *
36* July 14 1993: comment changes and addition of #pragma fenv_access. *
37* May 06 1997: port of the ibm/taligent ceil and floor routines. *
38* April 11 2001: first port to os x using gcc. *
39* June 13 2001: replaced __setflm with FEGETENVD/FESETENVD; *
40* replaced DblInHex typedef with hexdouble. *
41* used standard exception symbols from fenv.h. *
42* Sept 06 2001: added #ifdef __ppc__. *
43* Sept 09 2001: added more comments. *
44* Sept 10 2001: added macros to detect PowerPC and correct compiler. *
45* Sept 17 2001: deleted "fp.h" and changed "fenv.h" to <fenv.h>. *
46* Sept 18 2001: added <CoreServices/CoreServices.h> to get <fenv.h>. *
47* October 08 2001: removed <CoreServices/CoreServices.h>. *
48* changed compiler errors to warnings. *
49* November 06 2001: commented out warning about Intel architectures. *
50* *
51* W A R N I N G: *
52* These routines require a 64-bit double precision IEEE-754 model. *
53* They are written for PowerPC only and are expecting the compiler *
54* to generate the correct sequence of multiply-add fused instructions. *
55* *
56* These routines are not intended for 32-bit Intel architectures. *
57* *
58* A version of gcc higher than 932 is required. *
59* *
60* GCC compiler options: *
61* optimization level 3 (-O3) *
62* -fschedule-insns -finline-functions -funroll-all-loops *
63* *
64*******************************************************************************/
65
66#include "math.h"
67#include "fp_private.h"
68#include "fenv_private.h"
69
70static const double twoTo52 = 0x1.0p+52; // 4503599627370496.0;
71static const float twoTo23 = 0x1.0p+23F; // 8388608.0;
72
73/*******************************************************************************
74* Ceil(x) returns the smallest integer not less than x. *
75*******************************************************************************/
76
77double ceil ( double x )
78{
79 register double y;
80 register int target;
81
82 register double FPR_absx, FPR_Two52, FPR_one, FPR_zero, FPR_Mzero, FPR_pi4;
83
84 FPR_absx = __FABS( x ); FPR_zero = 0.0;
85 FPR_Two52 = twoTo52; FPR_one = 1.0;
86 __ENSURE( FPR_zero, FPR_Two52, FPR_one );
87 FPR_pi4 = M_PI_4; FPR_Mzero = -0.0;
88 target = ( x > FPR_zero ); __ENSURE( FPR_Mzero, FPR_Two52, FPR_pi4 );
89
90 if (likely( FPR_absx < FPR_Two52 ))
91/*******************************************************************************
92* Is |x| < 2.0^52? *
93*******************************************************************************/
94 {
95 if ( FPR_absx < FPR_one )
96/*******************************************************************************
97* Is |x| < 1.0? *
98*******************************************************************************/
99 {
100 if ( x == FPR_zero ) // zero x is exact case
101 return ( x );
102 else
103 { // inexact case
104 __PROG_INEXACT( FPR_pi4 );
105 if ( target )
106 return ( FPR_one );
107 else
108 return ( FPR_Mzero );
109 }
110 }
111/*******************************************************************************
112* Is 1.0 < |x| < 2.0^52? *
113*******************************************************************************/
114 if ( target )
115 {
116 y = ( x + FPR_Two52 ) - FPR_Two52; // round at binary pt.
117 if ( y < x )
118 return ( y + FPR_one );
119 else
120 return ( y );
121 }
122
123 else
124 {
125 y = ( x - FPR_Two52 ) + FPR_Two52; // round at binary pt.
126 if ( y < x )
127 return ( y + FPR_one );
128 else
129 return ( y );
130 }
131 }
132/*******************************************************************************
133* |x| >= 2.0^52 or x is a NaN. *
134*******************************************************************************/
135 return ( x );
136}
137
138float ceilf ( float x )
139{
140 register float y;
141 register int target;
142
143 register float FPR_absx, FPR_Two23, FPR_one, FPR_zero, FPR_Mzero;
144 register double FPR_pi4;
145
146 FPR_absx = __FABSF( x ); FPR_zero = 0.0f;
147 FPR_Two23 = twoTo23; FPR_one = 1.0f;
148 __ENSURE( FPR_zero, FPR_Two23, FPR_one );
149 FPR_pi4 = M_PI_4; FPR_Mzero = -0.0f;
150 target = ( x > FPR_zero ); __ENSURE( FPR_Mzero, FPR_Two23, FPR_pi4 );
151
152 if (likely( FPR_absx < FPR_Two23 ))
153/*******************************************************************************
154* Is |x| < 2.0^23? *
155*******************************************************************************/
156 {
157 if ( FPR_absx < FPR_one )
158/*******************************************************************************
159* Is |x| < 1.0? *
160*******************************************************************************/
161 {
162 if ( x == FPR_zero ) // zero x is exact case
163 return ( x );
164 else
165 { // inexact case
166 __PROG_INEXACT( FPR_pi4 );
167 if ( target )
168 return ( FPR_one );
169 else
170 return ( FPR_Mzero );
171 }
172 }
173/*******************************************************************************
174* Is 1.0 < |x| < 2.0^23? *
175*******************************************************************************/
176 if ( target )
177 {
178 y = ( x + FPR_Two23 ) - FPR_Two23; // round at binary pt.
179 if ( y < x )
180 return ( y + FPR_one );
181 else
182 return ( y );
183 }
184
185 else
186 {
187 y = ( x - FPR_Two23 ) + FPR_Two23; // round at binary pt.
188 if ( y < x )
189 return ( y + FPR_one );
190 else
191 return ( y );
192 }
193 }
194/*******************************************************************************
195* |x| >= 2.0^23 or x is a NaN. *
196*******************************************************************************/
197 return ( x );
198}
199
200/*******************************************************************************
201* Floor(x) returns the largest integer not greater than x. *
202*******************************************************************************/
203
204double floor ( double x )
205{
206 register double y;
207 register int target;
208
209 register double FPR_absx, FPR_Two52, FPR_one, FPR_zero, FPR_Mone, FPR_pi4;
210
211 FPR_absx = __FABS( x ); FPR_zero = 0.0;
212 FPR_Two52 = twoTo52; FPR_one = 1.0;
213 __ENSURE( FPR_zero, FPR_Two52, FPR_one );
214 FPR_pi4 = M_PI_4; FPR_Mone = -1.0;
215 target = ( x > FPR_zero ); __ENSURE( FPR_Mone, FPR_zero, FPR_pi4 );
216
217 if (likely( FPR_absx < FPR_Two52 ))
218/*******************************************************************************
219* Is |x| < 2.0^52? *
220*******************************************************************************/
221 {
222 if ( FPR_absx < FPR_one )
223/*******************************************************************************
224* Is |x| < 1.0? *
225*******************************************************************************/
226 {
227 if ( x == FPR_zero ) // zero x is exact case
228 return ( x );
229 else
230 { // inexact case
231 __PROG_INEXACT( FPR_pi4 );
232 if ( target )
233 return ( FPR_zero );
234 else
235 return ( FPR_Mone );
236 }
237 }
238/*******************************************************************************
239* Is 1.0 < |x| < 2.0^52? *
240*******************************************************************************/
241 if ( target )
242 {
243 y = ( x + FPR_Two52 ) - FPR_Two52; // round at binary pt.
244 if ( y > x )
245 return ( y - FPR_one );
246 else
247 return ( y );
248 }
249
250 else
251 {
252 y = ( x - FPR_Two52 ) + FPR_Two52; // round at binary pt.
253 if ( y > x )
254 return ( y - FPR_one );
255 else
256 return ( y );
257 }
258 }
259/*******************************************************************************
260* |x| >= 2.0^52 or x is a NaN. *
261*******************************************************************************/
262 return ( x );
263}
264
265float floorf ( float x )
266{
267 register float y;
268 register int target;
269
270 register float FPR_absx, FPR_Two23, FPR_one, FPR_zero, FPR_Mone;
271 register double FPR_pi4;
272
273 FPR_absx = __FABSF( x ); FPR_zero = 0.0f;
274 FPR_Two23 = twoTo23; FPR_one = 1.0f;
275 __ENSURE( FPR_zero, FPR_Two23, FPR_one );
276 FPR_pi4 = M_PI_4; FPR_Mone = -1.0f;
277 target = ( x > FPR_zero ); __ENSURE( FPR_Mone, FPR_zero, FPR_pi4 );
278
279 if (likely( FPR_absx < FPR_Two23 ))
280/*******************************************************************************
281* Is |x| < 2.0^23? *
282*******************************************************************************/
283 {
284 if ( FPR_absx < FPR_one )
285/*******************************************************************************
286* Is |x| < 1.0? *
287*******************************************************************************/
288 {
289 if ( x == FPR_zero ) // zero x is exact case
290 return ( x );
291 else
292 { // inexact case
293 __PROG_INEXACT( FPR_pi4 );
294 if ( target )
295 return ( FPR_zero );
296 else
297 return ( FPR_Mone );
298 }
299 }
300/*******************************************************************************
301* Is 1.0 < |x| < 2.0^23? *
302*******************************************************************************/
303 if ( target )
304 {
305 y = ( x + FPR_Two23 ) - FPR_Two23; // round at binary pt.
306 if ( y > x )
307 return ( y - FPR_one );
308 else
309 return ( y );
310 }
311
312 else
313 {
314 y = ( x - FPR_Two23 ) + FPR_Two23; // round at binary pt.
315 if ( y > x )
316 return ( y - FPR_one );
317 else
318 return ( y );
319 }
320 }
321/*******************************************************************************
322* |x| >= 2.0^23 or x is a NaN. *
323*******************************************************************************/
324 return ( x );
325}