blob: 790dd7ba032519135ae308c9d6a1fea2656d7fc2 (
plain)
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
|
/*******************************************************************************
* *
* File ceilfloor.c, *
* Function ceil(x) and floor(x), *
* Implementation of ceil and floor for the PowerPC. *
* *
* Copyright © 1991 Apple Computer, Inc. All rights reserved. *
* *
* Written by Ali Sazegari, started on November 1991, *
* *
* based on math.h, library code for Macintoshes with a 68881/68882 *
* by Jim Thomas. *
* *
* W A R N I N G: This routine expects a 64 bit double model. *
* *
* December 03 1992: first rs6000 port. *
* July 14 1993: comment changes and addition of #pragma fenv_access. *
* May 06 1997: port of the ibm/taligent ceil and floor routines. *
* April 11 2001: first port to os x using gcc. *
* June 13 2001: replaced __setflm with in-line assembly *
* *
*******************************************************************************/
static const double twoTo52 = 4503599627370496.0;
static const unsigned long signMask = 0x80000000ul;
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;
/*******************************************************************************
* Functions needed for the computation. *
*******************************************************************************/
/*******************************************************************************
* Ceil(x) returns the smallest integer not less than x. *
*******************************************************************************/
double ceil ( double x )
{
DblInHex xInHex,OldEnvironment;
register double y;
register unsigned long int xhi;
register int target;
xInHex.dbl = x;
xhi = xInHex.words.hi & 0x7fffffffUL; // xhi is the high half of |x|
target = ( xInHex.words.hi < signMask );
if ( xhi < 0x43300000ul )
/*******************************************************************************
* Is |x| < 2.0^52? *
*******************************************************************************/
{
if ( xhi < 0x3ff00000ul )
/*******************************************************************************
* Is |x| < 1.0? *
*******************************************************************************/
{
if ( ( xhi | xInHex.words.lo ) == 0ul ) // zero x is exact case
return ( x );
else
{ // inexact case
asm ("mffs %0" : "=f" (OldEnvironment.dbl));
OldEnvironment.words.lo |= 0x02000000ul;
asm ("mtfsf 255,%0" : /*NULLOUT*/ : /*IN*/ "f" ( OldEnvironment.dbl ));
if ( target )
return ( 1.0 );
else
return ( -0.0 );
}
}
/*******************************************************************************
* Is 1.0 < |x| < 2.0^52? *
*******************************************************************************/
if ( target )
{
y = ( x + twoTo52 ) - twoTo52; // round at binary pt.
if ( y < x )
return ( y + 1.0 );
else
return ( y );
}
else
{
y = ( x - twoTo52 ) + twoTo52; // round at binary pt.
if ( y < x )
return ( y + 1.0 );
else
return ( y );
}
}
/*******************************************************************************
* |x| >= 2.0^52 or x is a NaN. *
*******************************************************************************/
return ( x );
}
|
