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