s_modf.c 16 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344
  1. /*******************************************************************************
  2. ** File: rndint.c
  3. **
  4. ** Contains: C source code for implementations of floating-point
  5. ** functions which round to integral value or format, as
  6. ** defined in header <fp.h>. In particular, this file
  7. ** contains implementations of functions rinttol, roundtol,
  8. ** modf and modfl. This file targets PowrPC or Power platforms.
  9. **
  10. ** Written by: A. Sazegari, Apple AltiVec Group
  11. ** Created originally by Jon Okada, Apple Numerics Group
  12. **
  13. ** Copyright: © 1992-2001 by Apple Computer, Inc., all rights reserved
  14. **
  15. ** Change History (most recent first):
  16. **
  17. ** 13 Jul 01 ram replaced --setflm calls with inline assembly
  18. ** 03 Mar 01 ali first port to os x using gcc, added the crucial __setflm
  19. ** definition.
  20. ** 1. removed double_t, put in double for now.
  21. ** 2. removed iclass from nearbyint.
  22. ** 3. removed wrong comments intrunc.
  23. ** 4.
  24. ** 13 May 97 ali made performance improvements in rint, rinttol, roundtol
  25. ** and trunc by folding some of the taligent ideas into this
  26. ** implementation. nearbyint is faster than the one in taligent,
  27. ** rint is more elegant, but slower by %30 than the taligent one.
  28. ** 09 Apr 97 ali deleted modfl and deferred to AuxiliaryDD.c
  29. ** 15 Sep 94 ali Major overhaul and performance improvements of all functions.
  30. ** 20 Jul 94 PAF New faster version
  31. ** 16 Jul 93 ali Added the modfl function.
  32. ** 18 Feb 93 ali Changed the return value of fenv functions
  33. ** feclearexcept and feraiseexcept to their new
  34. ** NCEG X3J11.1/93-001 definitions.
  35. ** 16 Dec 92 JPO Removed __itrunc implementation to a
  36. ** separate file.
  37. ** 15 Dec 92 JPO Added __itrunc implementation and modified
  38. ** rinttol to include conversion from double
  39. ** to long int format. Modified roundtol to
  40. ** call __itrunc.
  41. ** 10 Dec 92 JPO Added modf (double) implementation.
  42. ** 04 Dec 92 JPO First created.
  43. **
  44. *******************************************************************************/
  45. #include <limits.h>
  46. #include <math.h>
  47. #include <endian.h>
  48. #define SET_INVALID 0x01000000UL
  49. typedef union
  50. {
  51. struct {
  52. #if (__BYTE_ORDER == __BIG_ENDIAN)
  53. unsigned long int hi;
  54. unsigned long int lo;
  55. #else
  56. unsigned long int lo;
  57. unsigned long int hi;
  58. #endif
  59. } words;
  60. double dbl;
  61. } DblInHex;
  62. static const unsigned long int signMask = 0x80000000ul;
  63. static const double twoTo52 = 4503599627370496.0;
  64. static const double doubleToLong = 4503603922337792.0; // 2^52
  65. static const DblInHex TOWARDZERO = {{ 0x00000000, 0x00000001 }};
  66. /*******************************************************************************
  67. * *
  68. * The function rinttol converts its double argument to integral value *
  69. * according to the current rounding direction and returns the result in *
  70. * long int format. This conversion signals invalid if the argument is a *
  71. * NaN or the rounded intermediate result is out of range of the *
  72. * destination long int format, and it delivers an unspecified result in *
  73. * this case. This function signals inexact if the rounded result is *
  74. * within range of the long int format but unequal to the operand. *
  75. * *
  76. *******************************************************************************/
  77. long int rinttol ( double x )
  78. {
  79. register double y;
  80. DblInHex argument, OldEnvironment;
  81. unsigned long int xHead;
  82. register long int target;
  83. argument.dbl = x;
  84. target = ( argument.words.hi < signMask ); // flag positive sign
  85. xHead = argument.words.hi & 0x7ffffffful; // high 32 bits of x
  86. if ( target )
  87. /*******************************************************************************
  88. * Sign of x is positive. *
  89. *******************************************************************************/
  90. {
  91. if ( xHead < 0x41dffffful )
  92. { // x is safely in long range
  93. y = ( x + twoTo52 ) - twoTo52; // round at binary point
  94. argument.dbl = y + doubleToLong; // force result into argument.words.lo
  95. return ( ( long ) argument.words.lo );
  96. }
  97. asm ("mffs %0" : "=f" (OldEnvironment.dbl)); // get environment
  98. if ( xHead > 0x41dffffful )
  99. { // x is safely out of long range
  100. OldEnvironment.words.lo |= SET_INVALID;
  101. asm ("mtfsf 255,%0" : /*NULLOUT*/ : /*IN*/ "f" ( OldEnvironment.dbl ));
  102. return ( LONG_MAX );
  103. }
  104. /*******************************************************************************
  105. * x > 0.0 and may or may not be out of range of long. *
  106. *******************************************************************************/
  107. y = ( x + twoTo52 ) - twoTo52; // do rounding
  108. if ( y > ( double ) LONG_MAX )
  109. { // out of range of long
  110. OldEnvironment.words.lo |= SET_INVALID;
  111. asm ("mtfsf 255,%0" : /*NULLOUT*/ : /*IN*/ "f" ( OldEnvironment.dbl ));
  112. return ( LONG_MAX );
  113. }
  114. argument.dbl = y + doubleToLong; // in range
  115. return ( ( long ) argument.words.lo ); // return result & flags
  116. }
  117. /*******************************************************************************
  118. * Sign of x is negative. *
  119. *******************************************************************************/
  120. if ( xHead < 0x41e00000ul )
  121. { // x is safely in long range
  122. y = ( x - twoTo52 ) + twoTo52;
  123. argument.dbl = y + doubleToLong;
  124. return ( ( long ) argument.words.lo );
  125. }
  126. asm ("mffs %0" : "=f" (OldEnvironment.dbl)); // get environment
  127. if ( xHead > 0x41e00000ul )
  128. { // x is safely out of long range
  129. OldEnvironment.words.lo |= SET_INVALID;
  130. asm ("mtfsf 255,%0" : /*NULLOUT*/ : /*IN*/ "f" ( OldEnvironment.dbl ));
  131. return ( LONG_MIN );
  132. }
  133. /*******************************************************************************
  134. * x < 0.0 and may or may not be out of range of long. *
  135. *******************************************************************************/
  136. y = ( x - twoTo52 ) + twoTo52; // do rounding
  137. if ( y < ( double ) LONG_MIN )
  138. { // out of range of long
  139. OldEnvironment.words.lo |= SET_INVALID;
  140. asm ("mtfsf 255,%0" : /*NULLOUT*/ : /*IN*/ "f" ( OldEnvironment.dbl ));
  141. return ( LONG_MIN );
  142. }
  143. argument.dbl = y + doubleToLong; // in range
  144. return ( ( long ) argument.words.lo ); // return result & flags
  145. }
  146. /*******************************************************************************
  147. * *
  148. * The function roundtol converts its double argument to integral format *
  149. * according to the "add half to the magnitude and chop" rounding mode of *
  150. * Pascal's Round function and FORTRAN's NINT function. This conversion *
  151. * signals invalid if the argument is a NaN or the rounded intermediate *
  152. * result is out of range of the destination long int format, and it *
  153. * delivers an unspecified result in this case. This function signals *
  154. * inexact if the rounded result is within range of the long int format but *
  155. * unequal to the operand. *
  156. * *
  157. *******************************************************************************/
  158. long int roundtol ( double x )
  159. {
  160. register double y, z;
  161. DblInHex argument, OldEnvironment;
  162. register unsigned long int xhi;
  163. register long int target;
  164. const DblInHex kTZ = {{ 0x0, 0x1 }};
  165. const DblInHex kUP = {{ 0x0, 0x2 }};
  166. argument.dbl = x;
  167. xhi = argument.words.hi & 0x7ffffffful; // high 32 bits of x
  168. target = ( argument.words.hi < signMask ); // flag positive sign
  169. if ( xhi > 0x41e00000ul )
  170. /*******************************************************************************
  171. * Is x is out of long range or NaN? *
  172. *******************************************************************************/
  173. {
  174. asm ("mffs %0" : "=f" (OldEnvironment.dbl)); // get environment
  175. OldEnvironment.words.lo |= SET_INVALID;
  176. asm ("mtfsf 255,%0" : /*NULLOUT*/ : /*IN*/ "f" ( OldEnvironment.dbl ));
  177. if ( target ) // pin result
  178. return ( LONG_MAX );
  179. else
  180. return ( LONG_MIN );
  181. }
  182. if ( target )
  183. /*******************************************************************************
  184. * Is sign of x is "+"? *
  185. *******************************************************************************/
  186. {
  187. if ( x < 2147483647.5 )
  188. /*******************************************************************************
  189. * x is in the range of a long. *
  190. *******************************************************************************/
  191. {
  192. y = ( x + doubleToLong ) - doubleToLong; // round at binary point
  193. if ( y != x )
  194. { // inexact case
  195. asm ("mffs %0" : "=f" (OldEnvironment.dbl)); // save environment
  196. asm ("mtfsf 255,%0" : /*NULLOUT*/ : /*IN*/ "f" ( kTZ.dbl )); // truncate rounding
  197. z = x + 0.5; // truncate x + 0.5
  198. argument.dbl = z + doubleToLong;
  199. asm ("mtfsf 255,%0" : /*NULLOUT*/ : /*IN*/ "f" ( OldEnvironment.dbl ));
  200. return ( ( long ) argument.words.lo );
  201. }
  202. argument.dbl = y + doubleToLong; // force result into argument.words.lo
  203. return ( ( long ) argument.words.lo ); // return long result
  204. }
  205. /*******************************************************************************
  206. * Rounded positive x is out of the range of a long. *
  207. *******************************************************************************/
  208. asm ("mffs %0" : "=f" (OldEnvironment.dbl));
  209. OldEnvironment.words.lo |= SET_INVALID;
  210. asm ("mtfsf 255,%0" : /*NULLOUT*/ : /*IN*/ "f" ( OldEnvironment.dbl ));
  211. return ( LONG_MAX ); // return pinned result
  212. }
  213. /*******************************************************************************
  214. * x < 0.0 and may or may not be out of the range of a long. *
  215. *******************************************************************************/
  216. if ( x > -2147483648.5 )
  217. /*******************************************************************************
  218. * x is in the range of a long. *
  219. *******************************************************************************/
  220. {
  221. y = ( x + doubleToLong ) - doubleToLong; // round at binary point
  222. if ( y != x )
  223. { // inexact case
  224. asm ("mffs %0" : "=f" (OldEnvironment.dbl)); // save environment
  225. asm ("mtfsf 255,%0" : /*NULLOUT*/ : /*IN*/ "f" ( kUP.dbl )); // round up
  226. z = x - 0.5; // truncate x - 0.5
  227. argument.dbl = z + doubleToLong;
  228. asm ("mtfsf 255,%0" : /*NULLOUT*/ : /*IN*/ "f" ( OldEnvironment.dbl ));
  229. return ( ( long ) argument.words.lo );
  230. }
  231. argument.dbl = y + doubleToLong;
  232. return ( ( long ) argument.words.lo ); // return long result
  233. }
  234. /*******************************************************************************
  235. * Rounded negative x is out of the range of a long. *
  236. *******************************************************************************/
  237. asm ("mffs %0" : "=f" (OldEnvironment.dbl));
  238. OldEnvironment.words.lo |= SET_INVALID;
  239. asm ("mtfsf 255,%0" : /*NULLOUT*/ : /*IN*/ "f" ( OldEnvironment.dbl ));
  240. return ( LONG_MIN ); // return pinned result
  241. }
  242. /*******************************************************************************
  243. * The modf family of functions separate a floating-point number into its *
  244. * fractional and integral parts, returning the fractional part and writing *
  245. * the integral part in floating-point format to the object pointed to by a *
  246. * pointer argument. If the input argument is integral or infinite in *
  247. * value, the return value is a zero with the sign of the input argument. *
  248. * The modf family of functions raises no floating-point exceptions. older *
  249. * implemenation set the INVALID flag due to signaling NaN input. *
  250. * *
  251. *******************************************************************************/
  252. /*******************************************************************************
  253. * modf is the double implementation. *
  254. *******************************************************************************/
  255. libm_hidden_proto(modf)
  256. double modf ( double x, double *iptr )
  257. {
  258. register double OldEnvironment, xtrunc;
  259. register unsigned long int xHead, signBit;
  260. DblInHex argument;
  261. argument.dbl = x;
  262. xHead = argument.words.hi & 0x7ffffffful; // |x| high bit pattern
  263. signBit = ( argument.words.hi & 0x80000000ul ); // isolate sign bit
  264. if (xHead == 0x7ff81fe0)
  265. signBit = signBit | 0;
  266. if ( xHead < 0x43300000ul )
  267. /*******************************************************************************
  268. * Is |x| < 2.0^53? *
  269. *******************************************************************************/
  270. {
  271. if ( xHead < 0x3ff00000ul )
  272. /*******************************************************************************
  273. * Is |x| < 1.0? *
  274. *******************************************************************************/
  275. {
  276. argument.words.hi = signBit; // truncate to zero
  277. argument.words.lo = 0ul;
  278. *iptr = argument.dbl;
  279. return ( x );
  280. }
  281. /*******************************************************************************
  282. * Is 1.0 < |x| < 2.0^52? *
  283. *******************************************************************************/
  284. asm ("mffs %0" : "=f" (OldEnvironment)); // save environment
  285. // round toward zero
  286. asm ("mtfsf 255,%0" : /*NULLOUT*/ : /*IN*/ "f" ( TOWARDZERO.dbl ));
  287. if ( signBit == 0ul ) // truncate to integer
  288. xtrunc = ( x + twoTo52 ) - twoTo52;
  289. else
  290. xtrunc = ( x - twoTo52 ) + twoTo52;
  291. // restore caller's env
  292. asm ("mtfsf 255,%0" : /*NULLOUT*/ : /*IN*/ "f" ( OldEnvironment ));
  293. *iptr = xtrunc; // store integral part
  294. if ( x != xtrunc ) // nonzero fraction
  295. return ( x - xtrunc );
  296. else
  297. { // zero with x's sign
  298. argument.words.hi = signBit;
  299. argument.words.lo = 0ul;
  300. return ( argument.dbl );
  301. }
  302. }
  303. *iptr = x; // x is integral or NaN
  304. if ( x != x ) // NaN is returned
  305. return x;
  306. else
  307. { // zero with x's sign
  308. argument.words.hi = signBit;
  309. argument.words.lo = 0ul;
  310. return ( argument.dbl );
  311. }
  312. }
  313. libm_hidden_def(modf)