floor.c 7.1 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531
  1. /* ceil()
  2. * floor()
  3. * frexp()
  4. * ldexp()
  5. * signbit()
  6. * isnan()
  7. * isfinite()
  8. *
  9. * Floating point numeric utilities
  10. *
  11. *
  12. *
  13. * SYNOPSIS:
  14. *
  15. * double ceil(), floor(), frexp(), ldexp();
  16. * int signbit(), isnan(), isfinite();
  17. * double x, y;
  18. * int expnt, n;
  19. *
  20. * y = floor(x);
  21. * y = ceil(x);
  22. * y = frexp( x, &expnt );
  23. * y = ldexp( x, n );
  24. * n = signbit(x);
  25. * n = isnan(x);
  26. * n = isfinite(x);
  27. *
  28. *
  29. *
  30. * DESCRIPTION:
  31. *
  32. * All four routines return a double precision floating point
  33. * result.
  34. *
  35. * floor() returns the largest integer less than or equal to x.
  36. * It truncates toward minus infinity.
  37. *
  38. * ceil() returns the smallest integer greater than or equal
  39. * to x. It truncates toward plus infinity.
  40. *
  41. * frexp() extracts the exponent from x. It returns an integer
  42. * power of two to expnt and the significand between 0.5 and 1
  43. * to y. Thus x = y * 2**expn.
  44. *
  45. * ldexp() multiplies x by 2**n.
  46. *
  47. * signbit(x) returns 1 if the sign bit of x is 1, else 0.
  48. *
  49. * These functions are part of the standard C run time library
  50. * for many but not all C compilers. The ones supplied are
  51. * written in C for either DEC or IEEE arithmetic. They should
  52. * be used only if your compiler library does not already have
  53. * them.
  54. *
  55. * The IEEE versions assume that denormal numbers are implemented
  56. * in the arithmetic. Some modifications will be required if
  57. * the arithmetic has abrupt rather than gradual underflow.
  58. */
  59. /*
  60. Cephes Math Library Release 2.8: June, 2000
  61. Copyright 1984, 1995, 2000 by Stephen L. Moshier
  62. */
  63. #include <math.h>
  64. #ifdef UNK
  65. /* ceil(), floor(), frexp(), ldexp() may need to be rewritten. */
  66. #undef UNK
  67. #if BIGENDIAN
  68. #define MIEEE 1
  69. #else
  70. #define IBMPC 1
  71. #endif
  72. #endif
  73. #ifdef DEC
  74. #define EXPMSK 0x807f
  75. #define MEXP 255
  76. #define NBITS 56
  77. #endif
  78. #ifdef IBMPC
  79. #define EXPMSK 0x800f
  80. #define MEXP 0x7ff
  81. #define NBITS 53
  82. #endif
  83. #ifdef MIEEE
  84. #define EXPMSK 0x800f
  85. #define MEXP 0x7ff
  86. #define NBITS 53
  87. #endif
  88. extern double MAXNUM, NEGZERO;
  89. #ifdef ANSIPROT
  90. double floor ( double );
  91. int isnan ( double );
  92. int isfinite ( double );
  93. double ldexp ( double, int );
  94. #else
  95. double floor();
  96. int isnan(), isfinite();
  97. double ldexp();
  98. #endif
  99. double ceil(x)
  100. double x;
  101. {
  102. double y;
  103. #ifdef UNK
  104. mtherr( "ceil", DOMAIN );
  105. return(0.0);
  106. #endif
  107. #ifdef NANS
  108. if( isnan(x) )
  109. return( x );
  110. #endif
  111. #ifdef INFINITIES
  112. if(!isfinite(x))
  113. return(x);
  114. #endif
  115. y = floor(x);
  116. if( y < x )
  117. y += 1.0;
  118. #ifdef MINUSZERO
  119. if( y == 0.0 && x < 0.0 )
  120. return( NEGZERO );
  121. #endif
  122. return(y);
  123. }
  124. /* Bit clearing masks: */
  125. static unsigned short bmask[] = {
  126. 0xffff,
  127. 0xfffe,
  128. 0xfffc,
  129. 0xfff8,
  130. 0xfff0,
  131. 0xffe0,
  132. 0xffc0,
  133. 0xff80,
  134. 0xff00,
  135. 0xfe00,
  136. 0xfc00,
  137. 0xf800,
  138. 0xf000,
  139. 0xe000,
  140. 0xc000,
  141. 0x8000,
  142. 0x0000,
  143. };
  144. double floor(x)
  145. double x;
  146. {
  147. union
  148. {
  149. double y;
  150. unsigned short sh[4];
  151. } u;
  152. unsigned short *p;
  153. int e;
  154. #ifdef UNK
  155. mtherr( "floor", DOMAIN );
  156. return(0.0);
  157. #endif
  158. #ifdef NANS
  159. if( isnan(x) )
  160. return( x );
  161. #endif
  162. #ifdef INFINITIES
  163. if(!isfinite(x))
  164. return(x);
  165. #endif
  166. #ifdef MINUSZERO
  167. if(x == 0.0L)
  168. return(x);
  169. #endif
  170. u.y = x;
  171. /* find the exponent (power of 2) */
  172. #ifdef DEC
  173. p = (unsigned short *)&u.sh[0];
  174. e = (( *p >> 7) & 0377) - 0201;
  175. p += 3;
  176. #endif
  177. #ifdef IBMPC
  178. p = (unsigned short *)&u.sh[3];
  179. e = (( *p >> 4) & 0x7ff) - 0x3ff;
  180. p -= 3;
  181. #endif
  182. #ifdef MIEEE
  183. p = (unsigned short *)&u.sh[0];
  184. e = (( *p >> 4) & 0x7ff) - 0x3ff;
  185. p += 3;
  186. #endif
  187. if( e < 0 )
  188. {
  189. if( u.y < 0.0 )
  190. return( -1.0 );
  191. else
  192. return( 0.0 );
  193. }
  194. e = (NBITS -1) - e;
  195. /* clean out 16 bits at a time */
  196. while( e >= 16 )
  197. {
  198. #ifdef IBMPC
  199. *p++ = 0;
  200. #endif
  201. #ifdef DEC
  202. *p-- = 0;
  203. #endif
  204. #ifdef MIEEE
  205. *p-- = 0;
  206. #endif
  207. e -= 16;
  208. }
  209. /* clear the remaining bits */
  210. if( e > 0 )
  211. *p &= bmask[e];
  212. if( (x < 0) && (u.y != x) )
  213. u.y -= 1.0;
  214. return(u.y);
  215. }
  216. double frexp( x, pw2 )
  217. double x;
  218. int *pw2;
  219. {
  220. union
  221. {
  222. double y;
  223. unsigned short sh[4];
  224. } u;
  225. int i;
  226. #ifdef DENORMAL
  227. int k;
  228. #endif
  229. short *q;
  230. u.y = x;
  231. #ifdef UNK
  232. mtherr( "frexp", DOMAIN );
  233. return(0.0);
  234. #endif
  235. #ifdef IBMPC
  236. q = (short *)&u.sh[3];
  237. #endif
  238. #ifdef DEC
  239. q = (short *)&u.sh[0];
  240. #endif
  241. #ifdef MIEEE
  242. q = (short *)&u.sh[0];
  243. #endif
  244. /* find the exponent (power of 2) */
  245. #ifdef DEC
  246. i = ( *q >> 7) & 0377;
  247. if( i == 0 )
  248. {
  249. *pw2 = 0;
  250. return(0.0);
  251. }
  252. i -= 0200;
  253. *pw2 = i;
  254. *q &= 0x807f; /* strip all exponent bits */
  255. *q |= 040000; /* mantissa between 0.5 and 1 */
  256. return(u.y);
  257. #endif
  258. #ifdef IBMPC
  259. i = ( *q >> 4) & 0x7ff;
  260. if( i != 0 )
  261. goto ieeedon;
  262. #endif
  263. #ifdef MIEEE
  264. i = *q >> 4;
  265. i &= 0x7ff;
  266. if( i != 0 )
  267. goto ieeedon;
  268. #ifdef DENORMAL
  269. #else
  270. *pw2 = 0;
  271. return(0.0);
  272. #endif
  273. #endif
  274. #ifndef DEC
  275. /* Number is denormal or zero */
  276. #ifdef DENORMAL
  277. if( u.y == 0.0 )
  278. {
  279. *pw2 = 0;
  280. return( 0.0 );
  281. }
  282. /* Handle denormal number. */
  283. do
  284. {
  285. u.y *= 2.0;
  286. i -= 1;
  287. k = ( *q >> 4) & 0x7ff;
  288. }
  289. while( k == 0 );
  290. i = i + k;
  291. #endif /* DENORMAL */
  292. ieeedon:
  293. i -= 0x3fe;
  294. *pw2 = i;
  295. *q &= 0x800f;
  296. *q |= 0x3fe0;
  297. return( u.y );
  298. #endif
  299. }
  300. double ldexp( x, pw2 )
  301. double x;
  302. int pw2;
  303. {
  304. union
  305. {
  306. double y;
  307. unsigned short sh[4];
  308. } u;
  309. short *q;
  310. int e;
  311. #ifdef UNK
  312. mtherr( "ldexp", DOMAIN );
  313. return(0.0);
  314. #endif
  315. u.y = x;
  316. #ifdef DEC
  317. q = (short *)&u.sh[0];
  318. e = ( *q >> 7) & 0377;
  319. if( e == 0 )
  320. return(0.0);
  321. #else
  322. #ifdef IBMPC
  323. q = (short *)&u.sh[3];
  324. #endif
  325. #ifdef MIEEE
  326. q = (short *)&u.sh[0];
  327. #endif
  328. while( (e = (*q & 0x7ff0) >> 4) == 0 )
  329. {
  330. if( u.y == 0.0 )
  331. {
  332. return( 0.0 );
  333. }
  334. /* Input is denormal. */
  335. if( pw2 > 0 )
  336. {
  337. u.y *= 2.0;
  338. pw2 -= 1;
  339. }
  340. if( pw2 < 0 )
  341. {
  342. if( pw2 < -53 )
  343. return(0.0);
  344. u.y /= 2.0;
  345. pw2 += 1;
  346. }
  347. if( pw2 == 0 )
  348. return(u.y);
  349. }
  350. #endif /* not DEC */
  351. e += pw2;
  352. /* Handle overflow */
  353. #ifdef DEC
  354. if( e > MEXP )
  355. return( MAXNUM );
  356. #else
  357. if( e >= MEXP )
  358. return( 2.0*MAXNUM );
  359. #endif
  360. /* Handle denormalized results */
  361. if( e < 1 )
  362. {
  363. #ifdef DENORMAL
  364. if( e < -53 )
  365. return(0.0);
  366. *q &= 0x800f;
  367. *q |= 0x10;
  368. /* For denormals, significant bits may be lost even
  369. when dividing by 2. Construct 2^-(1-e) so the result
  370. is obtained with only one multiplication. */
  371. u.y *= ldexp(1.0, e-1);
  372. return(u.y);
  373. #else
  374. return(0.0);
  375. #endif
  376. }
  377. else
  378. {
  379. #ifdef DEC
  380. *q &= 0x807f; /* strip all exponent bits */
  381. *q |= (e & 0xff) << 7;
  382. #else
  383. *q &= 0x800f;
  384. *q |= (e & 0x7ff) << 4;
  385. #endif
  386. return(u.y);
  387. }
  388. }
  389. /**********************************************************************/
  390. /*
  391. * trunc is just a slightly modified version of floor above.
  392. */
  393. double trunc(double x)
  394. {
  395. union {
  396. double y;
  397. unsigned short sh[4];
  398. } u;
  399. unsigned short *p;
  400. int e;
  401. #ifdef UNK
  402. mtherr( "trunc", DOMAIN );
  403. return(0.0);
  404. #endif
  405. #ifdef NANS
  406. if( isnan(x) )
  407. return( x );
  408. #endif
  409. #ifdef INFINITIES
  410. if(!isfinite(x))
  411. return(x);
  412. #endif
  413. #ifdef MINUSZERO
  414. if(x == 0.0L)
  415. return(x);
  416. #endif
  417. u.y = x;
  418. /* find the exponent (power of 2) */
  419. #ifdef DEC
  420. p = (unsigned short *)&u.sh[0];
  421. e = (( *p >> 7) & 0377) - 0201;
  422. p += 3;
  423. #endif
  424. #ifdef IBMPC
  425. p = (unsigned short *)&u.sh[3];
  426. e = (( *p >> 4) & 0x7ff) - 0x3ff;
  427. p -= 3;
  428. #endif
  429. #ifdef MIEEE
  430. p = (unsigned short *)&u.sh[0];
  431. e = (( *p >> 4) & 0x7ff) - 0x3ff;
  432. p += 3;
  433. #endif
  434. if( e < 0 )
  435. return( 0.0 );
  436. e = (NBITS -1) - e;
  437. /* clean out 16 bits at a time */
  438. while( e >= 16 )
  439. {
  440. #ifdef IBMPC
  441. *p++ = 0;
  442. #endif
  443. #ifdef DEC
  444. *p-- = 0;
  445. #endif
  446. #ifdef MIEEE
  447. *p-- = 0;
  448. #endif
  449. e -= 16;
  450. }
  451. /* clear the remaining bits */
  452. if( e > 0 )
  453. *p &= bmask[e];
  454. return(u.y);
  455. }