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
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
|
/*******************************************************************************
** 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 rint, nearbyint,
** rinttol, round, roundtol, trunc, modf and modfl. This file
** targets PowerPC 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>
#if !defined(__ppc__)
#define asm(x)
#endif
#define SET_INVALID 0x01000000UL
typedef union
{
struct {
#if defined(__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 Huge = {{ 0x7FF00000, 0x00000000 }};
static const DblInHex TOWARDZERO = {{ 0x00000000, 0x00000001 }};
/*******************************************************************************
* *
* The function rint rounds its double argument to integral value *
* according to the current rounding direction and returns the result in *
* double format. This function signals inexact if an ordered return *
* value is not equal to the operand. *
* *
********************************************************************************
* *
* This function calls: fabs. *
* *
*******************************************************************************/
/*******************************************************************************
* First, an elegant implementation. *
********************************************************************************
*
*double rint ( double x )
* {
* double y;
*
* y = twoTo52.fval;
*
* if ( fabs ( x ) >= y ) // huge case is exact
* return x;
* if ( x < 0 ) y = -y; // negative case
* y = ( x + y ) - y; // force rounding
* if ( y == 0.0 ) // zero results mirror sign of x
* y = copysign ( y, x );
* return ( y );
* }
********************************************************************************
* Now a bit twidling version that is about %30 faster. *
*******************************************************************************/
#if defined(__ppc__)
double rint ( double x )
{
DblInHex argument;
register double y;
unsigned long int xHead;
register long int target;
argument.dbl = x;
xHead = argument.words.hi & 0x7fffffffUL; // xHead <- high half of |x|
target = ( argument.words.hi < signMask ); // flags positive sign
if ( xHead < 0x43300000ul )
/*******************************************************************************
* Is |x| < 2.0^52? *
*******************************************************************************/
{
if ( xHead < 0x3ff00000ul )
/*******************************************************************************
* Is |x| < 1.0? *
*******************************************************************************/
{
if ( target )
y = ( x + twoTo52 ) - twoTo52; // round at binary point
else
y = ( x - twoTo52 ) + twoTo52; // round at binary point
if ( y == 0.0 )
{ // fix sign of zero result
if ( target )
return ( 0.0 );
else
return ( -0.0 );
}
return y;
}
/*******************************************************************************
* Is 1.0 < |x| < 2.0^52? *
*******************************************************************************/
if ( target )
return ( ( x + twoTo52 ) - twoTo52 ); // round at binary pt.
else
return ( ( x - twoTo52 ) + twoTo52 );
}
/*******************************************************************************
* |x| >= 2.0^52 or x is a NaN. *
*******************************************************************************/
return ( x );
}
#endif /* __ppc__ */
/*******************************************************************************
* *
* The function nearbyint rounds its double argument to integral value *
* according to the current rounding direction and returns the result in *
* double format. This function does not signal inexact. *
* *
********************************************************************************
* *
* This function calls fabs and copysign. *
* *
*******************************************************************************/
double nearbyint ( double x )
{
double y;
#if defined(__ppc__)
double OldEnvironment;
#endif /* __ppc__ */
y = twoTo52;
asm ("mffs %0" : "=f" (OldEnvironment)); /* get the environement */
if ( fabs ( x ) >= y ) /* huge case is exact */
return x;
if ( x < 0 ) y = -y; /* negative case */
y = ( x + y ) - y; /* force rounding */
if ( y == 0.0 ) /* zero results mirror sign of x */
y = copysign ( y, x );
// restore old flags
asm ("mtfsf 255,%0" : /*NULLOUT*/ : /*IN*/ "f" ( OldEnvironment ));
return ( y );
}
/*******************************************************************************
* *
* 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 round rounds its double argument to integral value *
* according to the "add half to the magnitude and truncate" rounding of *
* Pascal's Round function and FORTRAN's ANINT function and returns the *
* result in double format. This function signals inexact if an ordered *
* return value is not equal to the operand. *
* *
*******************************************************************************/
double round ( double x )
{
DblInHex argument, OldEnvironment;
register double y, z;
register unsigned long int xHead;
register long int target;
argument.dbl = x;
xHead = argument.words.hi & 0x7fffffffUL; // xHead <- high half of |x|
target = ( argument.words.hi < signMask ); // flag positive sign
if ( xHead < 0x43300000ul )
/*******************************************************************************
* Is |x| < 2.0^52? *
*******************************************************************************/
{
if ( xHead < 0x3ff00000ul )
/*******************************************************************************
* Is |x| < 1.0? *
*******************************************************************************/
{
asm ("mffs %0" : "=f" (OldEnvironment.dbl)); // get environment
if ( xHead < 0x3fe00000ul )
/*******************************************************************************
* Is |x| < 0.5? *
*******************************************************************************/
{
if ( ( xHead | argument.words.lo ) != 0ul )
OldEnvironment.words.lo |= 0x02000000ul;
asm ("mtfsf 255,%0" : /*NULLOUT*/ : /*IN*/ "f" ( OldEnvironment.dbl ));
if ( target )
return ( 0.0 );
else
return ( -0.0 );
}
/*******************************************************************************
* Is 0.5 |x| < 1.0? *
*******************************************************************************/
OldEnvironment.words.lo |= 0x02000000ul;
asm ("mtfsf 255,%0" : /*NULLOUT*/ : /*IN*/ "f" ( OldEnvironment.dbl ));
if ( target )
return ( 1.0 );
else
return ( -1.0 );
}
/*******************************************************************************
* Is 1.0 < |x| < 2.0^52? *
*******************************************************************************/
if ( target )
{ // positive x
y = ( x + twoTo52 ) - twoTo52; // round at binary point
if ( y == x ) // exact case
return ( x );
z = x + 0.5; // inexact case
y = ( z + twoTo52 ) - twoTo52; // round at binary point
if ( y > z )
return ( y - 1.0 );
else
return ( y );
}
/*******************************************************************************
* Is x < 0? *
*******************************************************************************/
else
{
y = ( x - twoTo52 ) + twoTo52; // round at binary point
if ( y == x )
return ( x );
z = x - 0.5;
y = ( z - twoTo52 ) + twoTo52; // round at binary point
if ( y < z )
return ( y + 1.0 );
else
return ( y );
}
}
/*******************************************************************************
* |x| >= 2.0^52 or x is a NaN. *
*******************************************************************************/
return ( x );
}
/*******************************************************************************
* *
* 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;
#if defined(__ppc__)
const DblInHex kTZ = {{ 0x0, 0x1 }};
const DblInHex kUP = {{ 0x0, 0x2 }};
#endif /* __ppc__ */
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 function trunc truncates its double argument to integral value *
* and returns the result in double format. This function signals *
* inexact if an ordered return value is not equal to the operand. *
* *
*******************************************************************************/
double trunc ( double x )
{
DblInHex argument,OldEnvironment;
register double y;
register unsigned long int xhi;
register long int target;
argument.dbl = x;
xhi = argument.words.hi & 0x7fffffffUL; // xhi <- high half of |x|
target = ( argument.words.hi < signMask ); // flag positive sign
if ( xhi < 0x43300000ul )
/*******************************************************************************
* Is |x| < 2.0^53? *
*******************************************************************************/
{
if ( xhi < 0x3ff00000ul )
/*******************************************************************************
* Is |x| < 1.0? *
*******************************************************************************/
{
if ( ( xhi | argument.words.lo ) != 0ul )
{ // raise deserved INEXACT
asm ("mffs %0" : "=f" (OldEnvironment.dbl));
OldEnvironment.words.lo |= 0x02000000ul;
asm ("mtfsf 255,%0" : /*NULLOUT*/ : /*IN*/ "f" ( OldEnvironment.dbl ));
}
if ( target ) // return properly signed zero
return ( 0.0 );
else
return ( -0.0 );
}
/*******************************************************************************
* Is 1.0 < |x| < 2.0^52? *
*******************************************************************************/
if ( target )
{
y = ( x + twoTo52 ) - twoTo52; // round at binary point
if ( y > x )
return ( y - 1.0 );
else
return ( y );
}
else
{
y = ( x - twoTo52 ) + twoTo52; // round at binary point.
if ( y < x )
return ( y + 1.0 );
else
return ( y );
}
}
/*******************************************************************************
* Is |x| >= 2.0^52 or x is a NaN. *
*******************************************************************************/
return ( x );
}
/*******************************************************************************
* 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. *
*******************************************************************************/
#if defined(__ppc__)
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 );
}
}
#endif /* __ppc__ */
|
