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