base/math/ceilfloor.c

Go to the documentation of this file.
00001 #if defined(__ppc__) 00002 /******************************************************************************* 00003 * * 00004 * File ceilfloor.c, * 00005 * Function ceil(x) and floor(x), * 00006 * Implementation of ceil and floor for the PowerPC. * 00007 * * 00008 * Copyright © 1991 Apple Computer, Inc. All rights reserved. * 00009 * * 00010 * Written by Ali Sazegari, started on November 1991, * 00011 * * 00012 * based on math.h, library code for Macintoshes with a 68881/68882 * 00013 * by Jim Thomas. * 00014 * * 00015 * W A R N I N G: This routine expects a 64 bit double model. * 00016 * * 00017 * December 03 1992: first rs6000 port. * 00018 * July 14 1993: comment changes and addition of #pragma fenv_access. * 00019 * May 06 1997: port of the ibm/taligent ceil and floor routines. * 00020 * April 11 2001: first port to os x using gcc. * 00021 * June 13 2001: replaced __setflm with in-line assembly * 00022 * * 00023 *******************************************************************************/ 00024 00025 #if !defined(__ppc__) 00026 #define asm(x) 00027 #endif 00028 00029 static const double twoTo52 = 4503599627370496.0; 00030 static const unsigned long signMask = 0x80000000ul; 00031 00032 typedef union 00033 { 00034 struct { 00035 #if defined(__BIG_ENDIAN__) 00036 unsigned long int hi; 00037 unsigned long int lo; 00038 #else 00039 unsigned long int lo; 00040 unsigned long int hi; 00041 #endif 00042 } words; 00043 double dbl; 00044 } DblInHex; 00045 00046 /******************************************************************************* 00047 * Functions needed for the computation. * 00048 *******************************************************************************/ 00049 00050 /******************************************************************************* 00051 * Ceil(x) returns the smallest integer not less than x. * 00052 *******************************************************************************/ 00053 00054 double ceil ( double x ) 00055 { 00056 DblInHex xInHex,OldEnvironment; 00057 register double y; 00058 register unsigned long int xhi; 00059 register int target; 00060 00061 xInHex.dbl = x; 00062 xhi = xInHex.words.hi & 0x7fffffffUL; // xhi is the high half of |x| 00063 target = ( xInHex.words.hi < signMask ); 00064 00065 if ( xhi < 0x43300000ul ) 00066 /******************************************************************************* 00067 * Is |x| < 2.0^52? * 00068 *******************************************************************************/ 00069 { 00070 if ( xhi < 0x3ff00000ul ) 00071 /******************************************************************************* 00072 * Is |x| < 1.0? * 00073 *******************************************************************************/ 00074 { 00075 if ( ( xhi | xInHex.words.lo ) == 0ul ) // zero x is exact case 00076 return ( x ); 00077 else 00078 { // inexact case 00079 asm ("mffs %0" : "=f" (OldEnvironment.dbl)); 00080 OldEnvironment.words.lo |= 0x02000000ul; 00081 asm ("mtfsf 255,%0" : /*NULLOUT*/ : /*IN*/ "f" ( OldEnvironment.dbl )); 00082 if ( target ) 00083 return ( 1.0 ); 00084 else 00085 return ( -0.0 ); 00086 } 00087 } 00088 /******************************************************************************* 00089 * Is 1.0 < |x| < 2.0^52? * 00090 *******************************************************************************/ 00091 if ( target ) 00092 { 00093 y = ( x + twoTo52 ) - twoTo52; // round at binary pt. 00094 if ( y < x ) 00095 return ( y + 1.0 ); 00096 else 00097 return ( y ); 00098 } 00099 00100 else 00101 { 00102 y = ( x - twoTo52 ) + twoTo52; // round at binary pt. 00103 if ( y < x ) 00104 return ( y + 1.0 ); 00105 else 00106 return ( y ); 00107 } 00108 } 00109 /******************************************************************************* 00110 * |x| >= 2.0^52 or x is a NaN. * 00111 *******************************************************************************/ 00112 return ( x ); 00113 } 00114 00115 /******************************************************************************* 00116 * Floor(x) returns the largest integer not greater than x. * 00117 *******************************************************************************/ 00118 00119 double floor ( double x ) 00120 { 00121 DblInHex xInHex,OldEnvironment; 00122 register double y; 00123 register unsigned long int xhi; 00124 register long int target; 00125 00126 xInHex.dbl = x; 00127 xhi = xInHex.words.hi & 0x7fffffffUL; // xhi is the high half of |x| 00128 target = ( xInHex.words.hi < signMask ); 00129 00130 if ( xhi < 0x43300000ul ) 00131 /******************************************************************************* 00132 * Is |x| < 2.0^52? * 00133 *******************************************************************************/ 00134 { 00135 if ( xhi < 0x3ff00000ul ) 00136 /******************************************************************************* 00137 * Is |x| < 1.0? * 00138 *******************************************************************************/ 00139 { 00140 if ( ( xhi | xInHex.words.lo ) == 0ul ) // zero x is exact case 00141 return ( x ); 00142 else 00143 { // inexact case 00144 asm ("mffs %0" : "=f" (OldEnvironment.dbl)); 00145 OldEnvironment.words.lo |= 0x02000000ul; 00146 asm ("mtfsf 255,%0" : /*NULLOUT*/ : /*IN*/ "f" ( OldEnvironment.dbl )); 00147 if ( target ) 00148 return ( 0.0 ); 00149 else 00150 return ( -1.0 ); 00151 } 00152 } 00153 /******************************************************************************* 00154 * Is 1.0 < |x| < 2.0^52? * 00155 *******************************************************************************/ 00156 if ( target ) 00157 { 00158 y = ( x + twoTo52 ) - twoTo52; // round at binary pt. 00159 if ( y > x ) 00160 return ( y - 1.0 ); 00161 else 00162 return ( y ); 00163 } 00164 00165 else 00166 { 00167 y = ( x - twoTo52 ) + twoTo52; // round at binary pt. 00168 if ( y > x ) 00169 return ( y - 1.0 ); 00170 else 00171 return ( y ); 00172 } 00173 } 00174 /******************************************************************************* 00175 * |x| >= 2.0^52 or x is a NaN. * 00176 *******************************************************************************/ 00177 return ( x ); 00178 } 00179 #endif /* __ppc__ */

Generated on Thu Nov 20 11:49:51 2008 for RTAI API by doxygen 1.3.8