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
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
|
/*******************************************************************************
** File: rndint.c
**
** Contains: C source code for implementations of floating-point
** functions which round to integral value or format, as
** defined in header <fp.h>. In particular, this file
** contains implementations of functions rinttol, roundtol,
** modf and modfl. This file targets PowrPC or Power platforms.
**
** Written by: A. Sazegari, Apple AltiVec Group
** Created originally by Jon Okada, Apple Numerics Group
**
** Copyright: © 1992-2001 by Apple Computer, Inc., all rights reserved
**
** Change History (most recent first):
**
** 13 Jul 01 ram replaced --setflm calls with inline assembly
** 03 Mar 01 ali first port to os x using gcc, added the crucial __setflm
** definition.
** 1. removed double_t, put in double for now.
** 2. removed iclass from nearbyint.
** 3. removed wrong comments intrunc.
** 4.
** 13 May 97 ali made performance improvements in rint, rinttol, roundtol
** and trunc by folding some of the taligent ideas into this
** implementation. nearbyint is faster than the one in taligent,
** rint is more elegant, but slower by %30 than the taligent one.
** 09 Apr 97 ali deleted modfl and deferred to AuxiliaryDD.c
** 15 Sep 94 ali Major overhaul and performance improvements of all functions.
** 20 Jul 94 PAF New faster version
** 16 Jul 93 ali Added the modfl function.
** 18 Feb 93 ali Changed the return value of fenv functions
** feclearexcept and feraiseexcept to their new
** NCEG X3J11.1/93-001 definitions.
** 16 Dec 92 JPO Removed __itrunc implementation to a
** separate file.
** 15 Dec 92 JPO Added __itrunc implementation and modified
** rinttol to include conversion from double
** to long int format. Modified roundtol to
** call __itrunc.
** 10 Dec 92 JPO Added modf (double) implementation.
** 04 Dec 92 JPO First created.
**
*******************************************************************************/
#include <limits.h>
#include <math.h>
#include <endian.h>
#define SET_INVALID 0x01000000UL
typedef union
{
struct {
#if (__BYTE_ORDER == __BIG_ENDIAN)
unsigned long int hi;
unsigned long int lo;
#else
unsigned long int lo;
unsigned long int hi;
#endif
} words;
double dbl;
} DblInHex;
static const unsigned long int signMask = 0x80000000ul;
static const double twoTo52 = 4503599627370496.0;
static const double doubleToLong = 4503603922337792.0; // 2^52
static const DblInHex TOWARDZERO = {{ 0x00000000, 0x00000001 }};
/*******************************************************************************
* *
* The function rinttol converts its double argument to integral value *
* according to the current rounding direction and returns the result in *
* long int format. This conversion signals invalid if the argument is a *
* NaN or the rounded intermediate result is out of range of the *
* destination long int format, and it delivers an unspecified result in *
* this case. This function signals inexact if the rounded result is *
* within range of the long int format but unequal to the operand. *
* *
*******************************************************************************/
long int rinttol ( double x )
{
register double y;
DblInHex argument, OldEnvironment;
unsigned long int xHead;
register long int target;
argument.dbl = x;
target = ( argument.words.hi < signMask ); // flag positive sign
xHead = argument.words.hi & 0x7ffffffful; // high 32 bits of x
if ( target )
/*******************************************************************************
* Sign of x is positive. *
*******************************************************************************/
{
if ( xHead < 0x41dffffful )
{ // x is safely in long range
y = ( x + twoTo52 ) - twoTo52; // round at binary point
argument.dbl = y + doubleToLong; // force result into argument.words.lo
return ( ( long ) argument.words.lo );
}
asm ("mffs %0" : "=f" (OldEnvironment.dbl)); // get environment
if ( xHead > 0x41dffffful )
{ // x is safely out of long range
OldEnvironment.words.lo |= SET_INVALID;
asm ("mtfsf 255,%0" : /*NULLOUT*/ : /*IN*/ "f" ( OldEnvironment.dbl ));
return ( LONG_MAX );
}
/*******************************************************************************
* x > 0.0 and may or may not be out of range of long. *
*******************************************************************************/
y = ( x + twoTo52 ) - twoTo52; // do rounding
if ( y > ( double ) LONG_MAX )
{ // out of range of long
OldEnvironment.words.lo |= SET_INVALID;
asm ("mtfsf 255,%0" : /*NULLOUT*/ : /*IN*/ "f" ( OldEnvironment.dbl ));
return ( LONG_MAX );
}
argument.dbl = y + doubleToLong; // in range
return ( ( long ) argument.words.lo ); // return result & flags
}
/*******************************************************************************
* Sign of x is negative. *
*******************************************************************************/
if ( xHead < 0x41e00000ul )
{ // x is safely in long range
y = ( x - twoTo52 ) + twoTo52;
argument.dbl = y + doubleToLong;
return ( ( long ) argument.words.lo );
}
asm ("mffs %0" : "=f" (OldEnvironment.dbl)); // get environment
if ( xHead > 0x41e00000ul )
{ // x is safely out of long range
OldEnvironment.words.lo |= SET_INVALID;
asm ("mtfsf 255,%0" : /*NULLOUT*/ : /*IN*/ "f" ( OldEnvironment.dbl ));
return ( LONG_MIN );
}
/*******************************************************************************
* x < 0.0 and may or may not be out of range of long. *
*******************************************************************************/
y = ( x - twoTo52 ) + twoTo52; // do rounding
if ( y < ( double ) LONG_MIN )
{ // out of range of long
OldEnvironment.words.lo |= SET_INVALID;
asm ("mtfsf 255,%0" : /*NULLOUT*/ : /*IN*/ "f" ( OldEnvironment.dbl ));
return ( LONG_MIN );
}
argument.dbl = y + doubleToLong; // in range
return ( ( long ) argument.words.lo ); // return result & flags
}
/*******************************************************************************
* *
* The function roundtol converts its double argument to integral format *
* according to the "add half to the magnitude and chop" rounding mode of *
* Pascal's Round function and FORTRAN's NINT function. This conversion *
* signals invalid if the argument is a NaN or the rounded intermediate *
* result is out of range of the destination long int format, and it *
* delivers an unspecified result in this case. This function signals *
* inexact if the rounded result is within range of the long int format but *
* unequal to the operand. *
* *
*******************************************************************************/
long int roundtol ( double x )
{
register double y, z;
DblInHex argument, OldEnvironment;
register unsigned long int xhi;
register long int target;
const DblInHex kTZ = {{ 0x0, 0x1 }};
const DblInHex kUP = {{ 0x0, 0x2 }};
argument.dbl = x;
xhi = argument.words.hi & 0x7ffffffful; // high 32 bits of x
target = ( argument.words.hi < signMask ); // flag positive sign
if ( xhi > 0x41e00000ul )
/*******************************************************************************
* Is x is out of long range or NaN? *
*******************************************************************************/
{
asm ("mffs %0" : "=f" (OldEnvironment.dbl)); // get environment
OldEnvironment.words.lo |= SET_INVALID;
asm ("mtfsf 255,%0" : /*NULLOUT*/ : /*IN*/ "f" ( OldEnvironment.dbl ));
if ( target ) // pin result
return ( LONG_MAX );
else
return ( LONG_MIN );
}
if ( target )
/*******************************************************************************
* Is sign of x is "+"? *
*******************************************************************************/
{
if ( x < 2147483647.5 )
/*******************************************************************************
* x is in the range of a long. *
*******************************************************************************/
{
y = ( x + doubleToLong ) - doubleToLong; // round at binary point
if ( y != x )
{ // inexact case
asm ("mffs %0" : "=f" (OldEnvironment.dbl)); // save environment
asm ("mtfsf 255,%0" : /*NULLOUT*/ : /*IN*/ "f" ( kTZ.dbl )); // truncate rounding
z = x + 0.5; // truncate x + 0.5
argument.dbl = z + doubleToLong;
asm ("mtfsf 255,%0" : /*NULLOUT*/ : /*IN*/ "f" ( OldEnvironment.dbl ));
return ( ( long ) argument.words.lo );
}
argument.dbl = y + doubleToLong; // force result into argument.words.lo
return ( ( long ) argument.words.lo ); // return long result
}
/*******************************************************************************
* Rounded positive x is out of the range of a long. *
*******************************************************************************/
asm ("mffs %0" : "=f" (OldEnvironment.dbl));
OldEnvironment.words.lo |= SET_INVALID;
asm ("mtfsf 255,%0" : /*NULLOUT*/ : /*IN*/ "f" ( OldEnvironment.dbl ));
return ( LONG_MAX ); // return pinned result
}
/*******************************************************************************
* x < 0.0 and may or may not be out of the range of a long. *
*******************************************************************************/
if ( x > -2147483648.5 )
/*******************************************************************************
* x is in the range of a long. *
*******************************************************************************/
{
y = ( x + doubleToLong ) - doubleToLong; // round at binary point
if ( y != x )
{ // inexact case
asm ("mffs %0" : "=f" (OldEnvironment.dbl)); // save environment
asm ("mtfsf 255,%0" : /*NULLOUT*/ : /*IN*/ "f" ( kUP.dbl )); // round up
z = x - 0.5; // truncate x - 0.5
argument.dbl = z + doubleToLong;
asm ("mtfsf 255,%0" : /*NULLOUT*/ : /*IN*/ "f" ( OldEnvironment.dbl ));
return ( ( long ) argument.words.lo );
}
argument.dbl = y + doubleToLong;
return ( ( long ) argument.words.lo ); // return long result
}
/*******************************************************************************
* Rounded negative x is out of the range of a long. *
*******************************************************************************/
asm ("mffs %0" : "=f" (OldEnvironment.dbl));
OldEnvironment.words.lo |= SET_INVALID;
asm ("mtfsf 255,%0" : /*NULLOUT*/ : /*IN*/ "f" ( OldEnvironment.dbl ));
return ( LONG_MIN ); // return pinned result
}
/*******************************************************************************
* The modf family of functions separate a floating-point number into its *
* fractional and integral parts, returning the fractional part and writing *
* the integral part in floating-point format to the object pointed to by a *
* pointer argument. If the input argument is integral or infinite in *
* value, the return value is a zero with the sign of the input argument. *
* The modf family of functions raises no floating-point exceptions. older *
* implemenation set the INVALID flag due to signaling NaN input. *
* *
*******************************************************************************/
/*******************************************************************************
* modf is the double implementation. *
*******************************************************************************/
libm_hidden_proto(modf)
double modf ( double x, double *iptr )
{
register double OldEnvironment, xtrunc;
register unsigned long int xHead, signBit;
DblInHex argument;
argument.dbl = x;
xHead = argument.words.hi & 0x7ffffffful; // |x| high bit pattern
signBit = ( argument.words.hi & 0x80000000ul ); // isolate sign bit
if (xHead == 0x7ff81fe0)
signBit = signBit | 0;
if ( xHead < 0x43300000ul )
/*******************************************************************************
* Is |x| < 2.0^53? *
*******************************************************************************/
{
if ( xHead < 0x3ff00000ul )
/*******************************************************************************
* Is |x| < 1.0? *
*******************************************************************************/
{
argument.words.hi = signBit; // truncate to zero
argument.words.lo = 0ul;
*iptr = argument.dbl;
return ( x );
}
/*******************************************************************************
* Is 1.0 < |x| < 2.0^52? *
*******************************************************************************/
asm ("mffs %0" : "=f" (OldEnvironment)); // save environment
// round toward zero
asm ("mtfsf 255,%0" : /*NULLOUT*/ : /*IN*/ "f" ( TOWARDZERO.dbl ));
if ( signBit == 0ul ) // truncate to integer
xtrunc = ( x + twoTo52 ) - twoTo52;
else
xtrunc = ( x - twoTo52 ) + twoTo52;
// restore caller's env
asm ("mtfsf 255,%0" : /*NULLOUT*/ : /*IN*/ "f" ( OldEnvironment ));
*iptr = xtrunc; // store integral part
if ( x != xtrunc ) // nonzero fraction
return ( x - xtrunc );
else
{ // zero with x's sign
argument.words.hi = signBit;
argument.words.lo = 0ul;
return ( argument.dbl );
}
}
*iptr = x; // x is integral or NaN
if ( x != x ) // NaN is returned
return x;
else
{ // zero with x's sign
argument.words.hi = signBit;
argument.words.lo = 0ul;
return ( argument.dbl );
}
}
libm_hidden_def(modf)
|
