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: rndint.c *
25* *
26* Contains: C source code for implementations of floating-point *
27* functions which round to integral value or format, as *
28* defined in C99. In particular, this file contains *
29* implementations of the following functions: *
30* rint, nearbyint, rinttol, round, roundtol, trunc and modf. *
31* *
32* Copyright � 1992-2001 by Apple Computer, Inc. All rights reserved. *
33* *
34* Written by Jon Okada, started on December 1992. *
35* Modified by Paul Finlayson (PAF) for MathLib v2. *
36* Modified by A. Sazegari (ali) for MathLib v3. *
37* Modified and ported by Robert A. Murley (ram) for Mac OS X. *
38* *
39* A MathLib v4 file. *
40* *
41* Change History (most recent first): *
42* *
43* 06 Nov 01 ram commented out warning about Intel architectures. *
44* changed i386 stubs to call abort(). *
45* 02 Nov 01 ram added stubs for i386 routines. *
46* 08 Oct 01 ram removed <Limits.h> and <CoreServices/CoreServices.h>. *
47* changed compiler errors to warnings. *
48* 05 Oct 01 ram added defines for INT32_MAX and INT32_MIN *
49* 18 Sep 01 ali <CoreServices/CoreServices.h> replaced "fp.h" & "fenv.h". *
50* 10 Sep 01 ali added more comments. *
51* 09 Sep 01 ali added macros to detect PowerPC and correct compiler. *
52* 28 Aug 01 ram added #ifdef __ppc__. *
53* 13 Jul 01 ram Replaced __setflm with FEGETENVD/FESETENVD. *
54* replaced DblInHex typedef with hexdouble. *
55* 03 Mar 01 ali first port to os x using gcc, added the crucial *
56* __setflm definition. *
57* 1. removed double_t, put in double for now. *
58* 2. removed iclass from nearbyint. *
59* 3. removed wrong comments in trunc. *
60* 13 May 97 ali made performance improvements in rint, rinttol, *
61* roundtol and trunc by folding some of the taligent *
62* ideas into this implementation. nearbyint is faster *
63* than the one in taligent, rint is more elegant, but *
64* slower by %30 than the taligent one. *
65* 09 Apr 97 ali deleted modfl and deferred to AuxiliaryDD.c *
66* 15 Sep 94 ali Major overhaul and performance improvements of all *
67* functions. *
68* 20 Jul 94 PAF New faster version. *
69* 16 Jul 93 ali Added the modfl function. *
70* 18 Feb 93 ali Changed the return value of fenv functions *
71* feclearexcept and feraiseexcept to their new *
72* NCEG X3J11.1/93-001 definitions. *
73* 16 Dec 92 JPO Removed __itrunc implementation to a separate file. *
74* 15 Dec 92 JPO Added __itrunc implementation and modified rinttol *
75* to include conversion from double to long int format. *
76* Modified roundtol to call __itrunc. *
77* 10 Dec 92 JPO Added modf (double) implementation. *
78* 04 Dec 92 JPO First created. *
79* *
80* W A R N I N G: *
81* These routines require a 64-bit double precision IEEE-754 model. *
82* They are written for PowerPC only and are expecting the compiler *
83* to generate the correct sequence of multiply-add fused instructions. *
84* *
85* These routines are not intended for 32-bit Intel architectures. *
86* *
87* A version of gcc higher than 932 is required. *
88* *
89* GCC compiler options: *
90* optimization level 3 (-O3) *
91* -fschedule-insns -finline-functions -funroll-all-loops *
92* *
93*******************************************************************************/
94
95#include "math.h"
96#include "fp_private.h"
97#include "fenv_private.h"
98#include "limits.h"
99
100
101static const hexdouble TOWARDZERO = HEXDOUBLE(0x00000000, 0x00000001);
102
103static const float twoTo23 = 0x1.0p+23f; // 8388608.0;
104
105#if !defined(BUILDING_FOR_CARBONCORE_LEGACY)
106
107static const double twoTo52 = 0x1.0p+52; // 4503599627370496.0;
108
109/*******************************************************************************
110* *
111* The function rint rounds its double argument to integral value *
112* according to the current rounding direction and returns the result in *
113* double format. This function signals inexact if an ordered return *
114* value is not equal to the operand. *
115* *
116*******************************************************************************/
117
118/*******************************************************************************
119* First, an elegant implementation. *
120********************************************************************************
121*
122*double rint ( double x )
123* {
124* double y;
125*
126* y = twoTo52.fval;
127*
128* if ( fabs ( x ) >= y ) // huge case is exact
129* return x;
130* if ( x < 0 ) y = -y; // negative case
131* y = ( x + y ) - y; // force rounding
132* if ( y == 0.0 ) // zero results mirror sign of x
133* y = copysign ( y, x );
134* return ( y );
135* }
136********************************************************************************
137* Now a bit twidling version that is about %30 faster. *
138*******************************************************************************/
139
140double rint ( double x )
141{
142 hexdouble argument;
143 register double y;
144 uint32_t xHead;
145 register int target;
146
147 argument.d = x;
148 __NOOP;
149 __NOOP;
150 __NOOP;
151
152 xHead = argument.i.hi & 0x7fffffff; // xHead <- high half of |x|
153 target = ( argument.i.hi < 0x80000000 ); // flags positive sign
154
155 if (likely( xHead < 0x43300000u ))
156/*******************************************************************************
157* Is |x| < 2.0^52? *
158*******************************************************************************/
159 {
160 if ( xHead < 0x3ff00000 )
161/*******************************************************************************
162* Is |x| < 1.0? *
163*******************************************************************************/
164 {
165 if ( target )
166 y = ( x + twoTo52 ) - twoTo52; // round at binary point
167 else
168 y = ( x - twoTo52 ) + twoTo52; // round at binary point
169 if ( y == 0.0 )
170 { // fix sign of zero result
171 if ( target )
172 return ( 0.0 );
173 else
174 return ( -0.0 );
175 }
176 return y;
177 }
178
179/*******************************************************************************
180* Is 1.0 < |x| < 2.0^52? *
181*******************************************************************************/
182
183 if ( target )
184 return ( ( x + twoTo52 ) - twoTo52 ); // round at binary pt.
185 else
186 return ( ( x - twoTo52 ) + twoTo52 );
187 }
188
189/*******************************************************************************
190* |x| >= 2.0^52 or x is a NaN. *
191*******************************************************************************/
192 return ( x );
193}
194
195float rintf ( float x )
196{
197 hexsingle argument;
198 register float y;
199 uint32_t xHead;
200 register int target;
201
202 argument.fval = x;
203 __NOOP;
204 __NOOP;
205 __NOOP;
206
207 xHead = argument.lval & 0x7fffffff; // xHead <- |x|
208 target = ( (uint32_t)argument.lval < 0x80000000u ); // flags positive sign
209
210 if (likely( xHead < 0x4b000000u ))
211/*******************************************************************************
212* Is |x| < 2.0^23? *
213*******************************************************************************/
214 {
215 if ( xHead < 0x3f800000 )
216/*******************************************************************************
217* Is |x| < 1.0? *
218*******************************************************************************/
219 {
220 if ( target )
221 y = ( x + twoTo23 ) - twoTo23; // round at binary point
222 else
223 y = ( x - twoTo23 ) + twoTo23; // round at binary point
224 if ( y == 0.0 )
225 { // fix sign of zero result
226 if ( target )
227 return ( 0.0f );
228 else
229 {
230#if defined(__GNUC__) && (__GNUC__<3) /* workaround gcc2.x botch of -0 return. */
231 volatile hexsingle zInHex;
232 zInHex.lval = 0x80000000;
233 return zInHex.fval;
234
235#else
236 return ( -0.0f );
237#endif
238 }
239 }
240 return y;
241 }
242
243/*******************************************************************************
244* Is 1.0 < |x| < 2.0^23? *
245*******************************************************************************/
246
247 if ( target )
248 return ( ( x + twoTo23 ) - twoTo23 ); // round at binary pt.
249 else
250 return ( ( x - twoTo23 ) + twoTo23 );
251 }
252
253/*******************************************************************************
254* |x| >= 2.0^23 or x is a NaN. *
255*******************************************************************************/
256 return ( x );
257}
258
259/*******************************************************************************
260* *
261* The function nearbyint rounds its double argument to integral value *
262* according to the current rounding direction and returns the result in *
263* double format. This function does not signal inexact. *
264* *
265* Functions used in this routine: *
266* fabs and copysign. *
267* *
268*******************************************************************************/
269
270double nearbyint ( double x )
271{
272 double y, OldEnvironment;
273
274 if (unlikely(x != x))
275 return x;
276
277 y = twoTo52;
278
279 FEGETENVD( OldEnvironment ); /* save the environement */
280
281 if (unlikely( fabs ( x ) >= y )) /* huge case is exact */
282 {
283 FESETENVD( OldEnvironment ); /* restore old environment */
284 return x;
285 }
286 if ( x < 0 ) y = -y; /* negative case */
287 y = ( x + y ) - y; /* force rounding */
288 if ( y == 0.0 ) /* zero results mirror sign of x */
289 y = copysign ( y, x );
290 FESETENVD( OldEnvironment ); /* restore old environment */
291 return ( y );
292}
293
294float nearbyintf ( float x )
295{
296 double OldEnvironment;
297 float y;
298
299 if (unlikely(x != x))
300 return x;
301
302 y = twoTo23;
303
304 FEGETENVD( OldEnvironment ); /* save the environement */
305
306 if (unlikely( fabsf ( x ) >= y )) /* huge case is exact */
307 {
308 FESETENVD( OldEnvironment ); /* restore old environment */
309 return x;
310 }
311 if ( x < 0 ) y = -y; /* negative case */
312 y = ( x + y ) - y; /* force rounding */
313 if ( y == 0.0 ) /* zero results mirror sign of x */
314 y = copysignf ( y, x );
315 FESETENVD( OldEnvironment ); /* restore old environment */
316 return ( y );
317}
318
319long int lrint ( double x )
320{
321 if (sizeof(long int) == 4) // PPC32 ABI
322 {
323 hexdouble hx;
324
325 asm volatile ("fctiw %0,%1" : "=f"(hx.d) : "f"(x));
326 __NOOP;
327 __NOOP;
328 __NOOP;
329 return hx.i.lo;
330 }
331 else // (sizeof(long int) == 8) LP64 ABI
332 {
333 union { double d; uint64_t i;} hx;
334 asm volatile ("fctid %0,%1" : "=f"(hx.d) : "f"(x));
335 __NOOP;
336 __NOOP;
337 __NOOP;
338 return hx.i;
339 }
340}
341
342long int lrintf ( float x )
343{
344 if (sizeof(long int) == 4) // PPC32 ABI
345 {
346 hexdouble hx;
347
348 asm volatile ("fctiw %0,%1" : "=f"(hx.d) : "f"(x));
349 __NOOP;
350 __NOOP;
351 __NOOP;
352 return hx.i.lo;
353 }
354 else // (sizeof(long int) == 8) LP64 ABI, PPC64 ISA
355 {
356 union { double d; uint64_t i;} hx;
357 asm volatile ("fctid %0,%1" : "=f"(hx.d) : "f"(x));
358 __NOOP;
359 __NOOP;
360 __NOOP;
361 return hx.i;
362 }
363}
364
365long long int llrint ( double x )
366{
367 double t;
368 long long int result;
369 double fenv;
370
371 if (sizeof(long int) == 8) // LP64 ABI, PPC64 ISA
372 {
373 union { double d; uint64_t i;} hx;
374 asm volatile ("fctid %0,%1" : "=f"(hx.d) : "f"(x));
375 __NOOP;
376 __NOOP;
377 __NOOP;
378 return hx.i;
379 }
380
381 if (unlikely(x != x))
382 {
383 feraiseexcept(FE_INVALID);
384 return LONG_LONG_MIN;
385 }
386
387 FEGETENVD(fenv);
388 t = rint ( x );
389 FESETENVD(fenv);
390
391 if ( t < (double)LONG_LONG_MIN )
392 {
393 feraiseexcept(FE_INVALID);
394 result = LONG_LONG_MIN;
395 }
396 else if ( t > (double)LONG_LONG_MAX )
397 {
398 feraiseexcept(FE_INVALID);
399 result = LONG_LONG_MAX;
400 }
401 else if (t != x)
402 {
403 feraiseexcept(FE_INEXACT);
404 result = (long long int) t;
405 }
406 else
407 {
408 result = (long long int) t;
409 }
410
411 return result;
412}
413
414long long int llrintf (float x)
415{
416 float t;
417 long long int result;
418 double fenv;
419
420 if (sizeof(long int) == 8) // LP64 ABI, PPC64 ISA
421 {
422 union { double d; uint64_t i;} hx;
423 asm volatile ("fctid %0,%1" : "=f"(hx.d) : "f"(x));
424 __NOOP;
425 __NOOP;
426 __NOOP;
427 return hx.i;
428 }
429
430 if (unlikely(x != x))
431 {
432 feraiseexcept(FE_INVALID);
433 return LONG_LONG_MIN;
434 }
435
436 FEGETENVD(fenv);
437 t = rintf ( x );
438 FESETENVD(fenv);
439
440 if ( t < (float)LONG_LONG_MIN )
441 {
442 feraiseexcept(FE_INVALID);
443 result = LONG_LONG_MIN;
444 }
445 else if ( t > (float)LONG_LONG_MAX )
446 {
447 feraiseexcept(FE_INVALID);
448 result = LONG_LONG_MAX;
449 }
450 else if (t != x)
451 {
452 feraiseexcept(FE_INEXACT);
453 result = (long long int) t;
454 }
455 else
456 {
457 result = (long long int) t;
458 }
459
460 return result;
461}
462
463/*******************************************************************************
464* *
465* The function round rounds its double argument to integral value *
466* according to the "add half to the magnitude and truncate" rounding of *
467* Pascal's Round function and FORTRAN's ANINT function and returns the *
468* result in double format. This function signals inexact if an ordered *
469* return value is not equal to the operand. *
470* *
471*******************************************************************************/
472
473double round ( double x )
474{
475 hexdouble argument, OldEnvironment;
476 register double y, z;
477 uint32_t xHead;
478 register int target;
479
480 argument.d = x;
481 __NOOP;
482 __NOOP;
483 __NOOP;
484
485 xHead = argument.i.hi & 0x7fffffff; // xHead <- high half of |x|
486 target = ( argument.i.hi < 0x80000000 ); // flag positive sign
487
488 if (likely( xHead < 0x43300000u ))
489/*******************************************************************************
490* Is |x| < 2.0^52? *
491*******************************************************************************/
492 {
493 if ( xHead < 0x3ff00000 )
494/*******************************************************************************
495* Is |x| < 1.0? *
496*******************************************************************************/
497 {
498 FEGETENVD( OldEnvironment.d ); // get environment
499 if ( xHead < 0x3fe00000u )
500/*******************************************************************************
501* Is |x| < 0.5? *
502*******************************************************************************/
503 {
504 if ( ( xHead | argument.i.lo ) != 0u )
505 OldEnvironment.i.lo |= FE_INEXACT;
506 FESETENVD_GRP( OldEnvironment.d );
507 if ( target )
508 return ( 0.0 );
509 else
510 return ( -0.0 );
511 }
512/*******************************************************************************
513* Is 0.5 � |x| < 1.0? *
514*******************************************************************************/
515 OldEnvironment.i.lo |= FE_INEXACT;
516 FESETENVD_GRP ( OldEnvironment.d );
517 if ( target )
518 return ( 1.0 );
519 else
520 return ( -1.0 );
521 }
522/*******************************************************************************
523* Is 1.0 < |x| < 2.0^52? *
524*******************************************************************************/
525 if ( target )
526 { // positive x
527 y = ( x + twoTo52 ) - twoTo52; // round at binary point
528 if ( y == x ) // exact case
529 return ( x );
530 z = x + 0.5; // inexact case
531 y = ( z + twoTo52 ) - twoTo52; // round at binary point
532 if ( y > z )
533 return ( y - 1.0 );
534 else
535 return ( y );
536 }
537
538/*******************************************************************************
539* Is x < 0? *
540*******************************************************************************/
541 else
542 {
543 y = ( x - twoTo52 ) + twoTo52; // round at binary point
544 if ( y == x )
545 return ( x );
546 z = x - 0.5;
547 y = ( z - twoTo52 ) + twoTo52; // round at binary point
548 if ( y < z )
549 return ( y + 1.0 );
550 else
551 return ( y );
552 }
553 }
554/*******************************************************************************
555* |x| >= 2.0^52 or x is a NaN. *
556*******************************************************************************/
557 return ( x );
558}
559
560float roundf ( float x )
561{
562 hexdouble OldEnvironment;
563 hexsingle argument;
564 register float y, z;
565 uint32_t xHead;
566 register int target;
567
568 argument.fval = x;
569 __NOOP;
570 __NOOP;
571 __NOOP;
572
573 xHead = argument.lval & 0x7fffffff; // xHead <- |x|
574 target = ( (uint32_t)argument.lval < 0x80000000u ); // flags positive sign
575
576 if (likely( xHead < 0x4b000000u ))
577/*******************************************************************************
578* Is |x| < 2.0^52? *
579*******************************************************************************/
580 {
581 if ( xHead < 0x3f800000u )
582/*******************************************************************************
583* Is |x| < 1.0? *
584*******************************************************************************/
585 {
586 FEGETENVD( OldEnvironment.d ); // get environment
587 if ( xHead < 0x3f000000u )
588/*******************************************************************************
589* Is |x| < 0.5? *
590*******************************************************************************/
591 {
592 if ( xHead != 0u )
593 OldEnvironment.i.lo |= FE_INEXACT;
594 FESETENVD_GRP( OldEnvironment.d );
595 if ( target )
596 return ( 0.0f );
597 else
598 {
599#if defined(__GNUC__) && (__GNUC__<3) /* workaround gcc2.x botch of -0 return. */
600 volatile hexsingle zInHex;
601 zInHex.lval = 0x80000000;
602 return zInHex.fval;
603
604#else
605 return ( -0.0f );
606#endif
607 }
608 }
609/*******************************************************************************
610* Is 0.5 � |x| < 1.0? *
611*******************************************************************************/
612 OldEnvironment.i.lo |= FE_INEXACT;
613 FESETENVD_GRP ( OldEnvironment.d );
614 if ( target )
615 return ( 1.0f );
616 else
617 return ( -1.0f );
618 }
619/*******************************************************************************
620* Is 1.0 < |x| < 2.0^23? *
621*******************************************************************************/
622 if ( target )
623 { // positive x
624 y = ( x + twoTo23 ) - twoTo23; // round at binary point
625 if ( y == x ) // exact case
626 return ( x );
627 z = x + 0.5f; // inexact case
628 y = ( z + twoTo23 ) - twoTo23; // round at binary point
629 if ( y > z )
630 return ( y - 1.0f );
631 else
632 return ( y );
633 }
634
635/*******************************************************************************
636* Is x < 0? *
637*******************************************************************************/
638 else
639 {
640 y = ( x - twoTo23 ) + twoTo23; // round at binary point
641 if ( y == x )
642 return ( x );
643 z = x - 0.5f;
644 y = ( z - twoTo23 ) + twoTo23; // round at binary point
645 if ( y < z )
646 return ( y + 1.0f );
647 else
648 return ( y );
649 }
650 }
651/*******************************************************************************
652* |x| >= 2.0^23 or x is a NaN. *
653*******************************************************************************/
654 return ( x );
655}
656
657/*******************************************************************************
658* *
659* The function roundtol converts its double argument to integral format *
660* according to the "add half to the magnitude and chop" rounding mode of *
661* Pascal's Round function and FORTRAN's NINT function. This conversion *
662* signals invalid if the argument is a NaN or the rounded intermediate *
663* result is out of range of the destination long int format, and it *
664* delivers an unspecified result in this case. This function signals *
665* inexact if the rounded result is within range of the long int format but *
666* unequal to the operand. *
667* *
668*******************************************************************************/
669
670// These work just as well for the LP64 ABI
671long int lround ( double x )
672{
673 double t;
674 long int result;
675 double fenv;
676
677 if (unlikely(x != x))
678 {
679 feraiseexcept(FE_INVALID);
680 return LONG_MAX;
681 }
682
683 FEGETENVD(fenv);
684 t = round ( x );
685 FESETENVD(fenv);
686
687 if ( t < (double)LONG_MIN )
688 {
689 feraiseexcept(FE_INVALID);
690 result = LONG_MIN;
691 }
692 else if ( t > (double)LONG_MAX )
693 {
694 feraiseexcept(FE_INVALID);
695 result = LONG_MAX;
696 }
697 else if (t != x)
698 {
699 feraiseexcept(FE_INEXACT);
700 result = (long int) t;
701 }
702 else
703 {
704 result = (long int) t;
705 }
706
707 return result;
708}
709
710long int lroundf ( float x )
711{
712 float t;
713 long int result;
714 double fenv;
715
716 if (unlikely(x != x))
717 {
718 feraiseexcept(FE_INVALID);
719 return LONG_MAX;
720 }
721
722 FEGETENVD(fenv);
723 t = roundf ( x );
724 FESETENVD(fenv);
725
726 if ( t < (float)LONG_MIN )
727 {
728 feraiseexcept(FE_INVALID);
729 result = LONG_MIN;
730 }
731 else if ( t > (float)LONG_MAX )
732 {
733 feraiseexcept(FE_INVALID);
734 result = LONG_MAX;
735 }
736 else if (t != x)
737 {
738 feraiseexcept(FE_INEXACT);
739 result = (long int) t;
740 }
741 else
742 {
743 result = (long int) t;
744 }
745
746 return result;
747}
748
749long long int llround ( double x )
750{
751 double t;
752 long long int result;
753 double fenv;
754
755 if (unlikely(x != x))
756 {
757 feraiseexcept(FE_INVALID);
758 return LONG_LONG_MAX;
759 }
760
761 FEGETENVD(fenv);
762 t = round ( x );
763 FESETENVD(fenv);
764
765 if ( t < (double)LONG_LONG_MIN )
766 {
767 feraiseexcept(FE_INVALID);
768 result = LONG_LONG_MIN;
769 }
770 else if ( t > (double)LONG_LONG_MAX )
771 {
772 feraiseexcept(FE_INVALID);
773 result = LONG_LONG_MAX;
774 }
775 else if (t != x)
776 {
777 feraiseexcept(FE_INEXACT);
778 result = (long long int) t;
779 }
780 else
781 {
782 result = (long long int) t;
783 }
784
785 return result;
786}
787
788long long int llroundf ( float x )
789{
790 float t;
791 long long int result;
792 double fenv;
793
794 if (unlikely(x != x))
795 {
796 feraiseexcept(FE_INVALID);
797 return LONG_LONG_MAX;
798 }
799
800 FEGETENVD(fenv);
801 t = roundf ( x );
802 FESETENVD(fenv);
803
804 if ( t < (float)LONG_LONG_MIN )
805 {
806 feraiseexcept(FE_INVALID);
807 result = LONG_LONG_MIN;
808 }
809 else if ( t > (float)LONG_LONG_MAX )
810 {
811 feraiseexcept(FE_INVALID);
812 result = LONG_LONG_MAX;
813 }
814 else if (t != x)
815 {
816 feraiseexcept(FE_INEXACT);
817 result = (long long int) t;
818 }
819 else
820 {
821 result = (long long int) t;
822 }
823
824 return result;
825}
826
827/*******************************************************************************
828* *
829* The function trunc truncates its double argument to integral value *
830* and returns the result in double format. This function signals *
831* inexact if an ordered return value is not equal to the operand. *
832* *
833*******************************************************************************/
834
835double trunc ( double x )
836{
837 hexdouble argument, OldEnvironment;
838 register double y;
839 uint32_t xhi;
840 register int target;
841
842 argument.d = x;
843 __NOOP;
844 __NOOP;
845 __NOOP;
846 xhi = argument.i.hi & 0x7fffffff; // xhi <- high half of |x|
847 target = ( argument.i.hi < 0x80000000 ); // flag positive sign
848
849 if (likely( xhi < 0x43300000u ))
850/*******************************************************************************
851* Is |x| < 2.0^52? *
852*******************************************************************************/
853 {
854 if ( xhi < 0x3ff00000u )
855/*******************************************************************************
856* Is |x| < 1.0? *
857*******************************************************************************/
858 {
859 if ( ( xhi | argument.i.lo ) != 0u )
860 { // raise deserved INEXACT
861 FEGETENVD_GRP( OldEnvironment.d );
862 OldEnvironment.i.lo |= FE_INEXACT;
863 FESETENVD_GRP( OldEnvironment.d );
864 }
865 if ( target ) // return properly signed zero
866 return ( 0.0 );
867 else
868 return ( -0.0 );
869 }
870/*******************************************************************************
871* Is 1.0 < |x| < 2.0^52? *
872*******************************************************************************/
873 if ( target )
874 {
875 y = ( x + twoTo52 ) - twoTo52; // round at binary point
876 if ( y > x )
877 return ( y - 1.0 );
878 else
879 return ( y );
880 }
881 else
882 {
883 y = ( x - twoTo52 ) + twoTo52; // round at binary point.
884 if ( y < x )
885 return ( y + 1.0 );
886 else
887 return ( y );
888 }
889 }
890/*******************************************************************************
891* Is |x| >= 2.0^52 or x is a NaN. *
892*******************************************************************************/
893 return ( x );
894}
895
896float truncf ( float x )
897{
898 hexdouble OldEnvironment;
899 hexsingle argument;
900 register float y;
901 uint32_t xhi;
902 register int target;
903
904 argument.fval = x;
905 __NOOP;
906 __NOOP;
907 __NOOP;
908
909 xhi = argument.lval & 0x7fffffff; // xhi <- |x|
910 target = ( (uint32_t)argument.lval < 0x80000000u ); // flags positive sign
911
912 if (likely( xhi < 0x4b000000u ))
913/*******************************************************************************
914* Is |x| < 2.0^23? *
915*******************************************************************************/
916 {
917 if ( xhi < 0x3f800000 )
918/*******************************************************************************
919* Is |x| < 1.0? *
920*******************************************************************************/
921 {
922 if ( xhi != 0u )
923 { // raise deserved INEXACT
924 FEGETENVD_GRP( OldEnvironment.d );
925 OldEnvironment.i.lo |= FE_INEXACT;
926 FESETENVD_GRP( OldEnvironment.d );
927 }
928 if ( target ) // return properly signed zero
929 return ( 0.0f );
930 else
931 {
932#if defined(__GNUC__) && (__GNUC__<3) /* workaround gcc2.x botch of -0 return. */
933 volatile hexsingle zInHex;
934 zInHex.lval = 0x80000000;
935 return zInHex.fval;
936
937#else
938 return ( -0.0f );
939#endif
940 }
941 }
942/*******************************************************************************
943* Is 1.0 < |x| < 2.0^23? *
944*******************************************************************************/
945 if ( target )
946 {
947 y = ( x + twoTo23 ) - twoTo23; // round at binary point
948 if ( y > x )
949 return ( y - 1.0f );
950 else
951 return ( y );
952 }
953 else
954 {
955 y = ( x - twoTo23 ) + twoTo23; // round at binary point.
956 if ( y < x )
957 return ( y + 1.0f );
958 else
959 return ( y );
960 }
961 }
962/*******************************************************************************
963* Is |x| >= 2.0^23 or x is a NaN. *
964*******************************************************************************/
965 return ( x );
966}
967
968/*******************************************************************************
969* *
970* The modf family of functions separate a floating-point number into its *
971* fractional and integral parts, returning the fractional part and writing *
972* the integral part in floating-point format to the object pointed to by a *
973* pointer argument. If the input argument is integral or infinite in *
974* value, the return value is a zero with the sign of the input argument. *
975* The modf family of functions raises no floating-point exceptions. older *
976* implemenation set the INVALID flag due to signaling NaN input. *
977* *
978* modf is the double implementation. *
979* *
980*******************************************************************************/
981
982double modf ( double x, double *iptr )
983{
984 register double OldEnvironment, xtrunc;
985 hexdouble argument;
986
987 register double FPR_negZero, FPR_zero, FPR_one, FPR_Two52, FPR_TowardZero, FPR_absx;
988
989 FPR_absx = __FABS( x ); FPR_Two52 = twoTo52;
990 FPR_one = 1.0; argument.d = x;
991
992 FPR_TowardZero = TOWARDZERO.d; FPR_zero = 0.0;
993
994 FPR_negZero = -0.0; __NOOP;
995
996 __ENSURE(FPR_zero, FPR_TowardZero, FPR_Two52); __ENSURE(FPR_zero, FPR_Two52, FPR_one);
997
998/*******************************************************************************
999* Is |x| < 2.0^52? *
1000*******************************************************************************/
1001 if (likely( FPR_absx < FPR_Two52 ))
1002 {
1003/*******************************************************************************
1004* Is |x| < 1.0? *
1005*******************************************************************************/
1006 if ( FPR_absx < FPR_one )
1007 {
1008 if ( argument.i.hi & 0x80000000 ) // isolate sign bit
1009 *iptr = FPR_negZero; // truncate to zero
1010 else
1011 *iptr = FPR_zero; // truncate to zero
1012 return ( x );
1013 }
1014/*******************************************************************************
1015* Is 1.0 < |x| < 2.0^52? *
1016*******************************************************************************/
1017 FEGETENVD( OldEnvironment ); // save environment
1018 // round toward zero
1019 FESETENVD( FPR_TowardZero );
1020 if ( x > FPR_zero ) // truncate to integer
1021 xtrunc = ( x + FPR_Two52 ) - FPR_Two52;
1022 else
1023 xtrunc = ( x - FPR_Two52 ) + FPR_Two52;
1024 *iptr = xtrunc; // store integral part
1025 // restore caller's env
1026 FESETENVD( OldEnvironment ); // restore environment
1027 if ( x != xtrunc ) // nonzero fraction
1028 return ( x - xtrunc );
1029 else
1030 { // zero with x's sign
1031 if ( argument.i.hi & 0x80000000 ) // isolate sign bit
1032 return FPR_negZero; // truncate to zero
1033 else
1034 return FPR_zero; // truncate to zero
1035 }
1036 }
1037
1038 *iptr = x; // x is integral or NaN
1039 if ( x != x ) // NaN is returned
1040 return x;
1041 else
1042 { // zero with x's sign
1043 if ( argument.i.hi & 0x80000000 ) // isolate sign bit
1044 return FPR_negZero; // truncate to zero
1045 else
1046 return FPR_zero; // truncate to zero
1047 }
1048}
1049
1050#else /* BUILDING_FOR_CARBONCORE_LEGACY */
1051
1052float modff ( float x, float *iptr )
1053{
1054 register double OldEnvironment;
1055 register float xtrunc;
1056 register uint32_t xHead, signBit;
1057 hexsingle argument;
1058
1059 argument.fval = x;
1060 __NOOP;
1061 __NOOP;
1062 __NOOP;
1063
1064 xHead = argument.lval & 0x7fffffff; // |x| high bit pattern
1065 signBit = ( argument.lval & 0x80000000 ); // isolate sign bit
1066
1067 if (likely( xHead < 0x4b000000u ))
1068/*******************************************************************************
1069* Is |x| < 2.0^23? *
1070*******************************************************************************/
1071 {
1072 if ( xHead < 0x3f800000 )
1073/*******************************************************************************
1074* Is |x| < 1.0? *
1075*******************************************************************************/
1076 {
1077 argument.lval = signBit; // truncate to zero
1078 __NOOP;
1079 __NOOP;
1080 __NOOP;
1081
1082 *iptr = argument.fval;
1083 return ( x );
1084 }
1085/*******************************************************************************
1086* Is 1.0 < |x| < 2.0^23? *
1087*******************************************************************************/
1088 FEGETENVD( OldEnvironment ); // save environment
1089 // round toward zero
1090 FESETENVD( TOWARDZERO.d );
1091 if ( signBit == 0u ) // truncate to integer
1092 xtrunc = ( x + twoTo23 ) - twoTo23;
1093 else
1094 xtrunc = ( x - twoTo23 ) + twoTo23;
1095 // restore caller's env
1096 FESETENVD( OldEnvironment ); // restore environment
1097 *iptr = xtrunc; // store integral part
1098 if ( x != xtrunc ) // nonzero fraction
1099 return ( x - xtrunc );
1100 else
1101 { // zero with x's sign
1102 argument.lval = signBit;
1103 __NOOP;
1104 __NOOP;
1105 __NOOP;
1106
1107 return ( argument.fval );
1108 }
1109 }
1110
1111 *iptr = x; // x is integral or NaN
1112 if ( x != x ) // NaN is returned
1113 return x;
1114 else
1115 { // zero with x's sign
1116 argument.lval = signBit;
1117 __NOOP;
1118 __NOOP;
1119 __NOOP;
1120
1121 return ( argument.fval );
1122 }
1123}
1124
1125#endif /* BUILDING_FOR_CARBONCORE_LEGACY */