home *** CD-ROM | disk | FTP | other *** search
/ Simtel MSDOS - Coast to Coast / simteldosarchivecoasttocoast2.iso / astrnomy / de118i.zip / FLOORL.C < prev    next >
C/C++ Source or Header  |  1993-01-31  |  5KB  |  356 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.  
  65. /* Turn UNK into something else */
  66. #if 0
  67. #if UNK
  68. #undef IBMPC
  69. #define IBMPC 1
  70. #undef UNK
  71. #endif
  72. #endif
  73.  
  74.  
  75. #define DENORMAL 1
  76.  
  77. #ifdef UNK
  78. char *unkmsg = "ceill(), floorl(), frexpl(), ldexpl() must be rewritten!\n";
  79. #endif
  80.  
  81. #ifdef IBMPC
  82. #define EXPMSK 0x800f
  83. #define MEXP 0x7ff
  84. #define NBITS 64
  85. #endif
  86.  
  87. #ifdef MIEEE
  88. #define EXPMSK 0x800f
  89. #define MEXP 0x7ff
  90. #define NBITS 64
  91. #endif
  92.  
  93. extern double MAXNUML;
  94.  
  95.  
  96. long double fabsl(x)
  97. long double x;
  98. {
  99. if( x < 0 )
  100.     return( -x );
  101. else
  102.     return( x );
  103. }
  104.  
  105.  
  106.  
  107. long double ceill(x)
  108. long double x;
  109. {
  110. long double y;
  111. long double floorl();
  112.  
  113. #ifdef UNK
  114. mtherr( "ceill", DOMAIN );
  115. return(0.0L);
  116. #endif
  117.  
  118. y = floorl(x);
  119. if( y < x )
  120.     y += 1.0L;
  121. return(y);
  122. }
  123.  
  124.  
  125.  
  126.  
  127. /* Bit clearing masks: */
  128.  
  129. static unsigned short bmask[] = {
  130. 0xffff,
  131. 0xfffe,
  132. 0xfffc,
  133. 0xfff8,
  134. 0xfff0,
  135. 0xffe0,
  136. 0xffc0,
  137. 0xff80,
  138. 0xff00,
  139. 0xfe00,
  140. 0xfc00,
  141. 0xf800,
  142. 0xf000,
  143. 0xe000,
  144. 0xc000,
  145. 0x8000,
  146. 0x0000,
  147. };
  148.  
  149.  
  150.  
  151.  
  152. long double floorl(x)
  153. long double x;
  154. {
  155. unsigned short *p;
  156. long double y;
  157. int e, i;
  158.  
  159. #ifdef UNK
  160. mtherr( "floor", DOMAIN );
  161. return(0.0L);
  162. #endif
  163.  
  164. y = x;
  165. /* find the exponent (power of 2) */
  166. #ifdef IBMPC
  167. p = (unsigned short *)&y + 4;
  168. e = (*p & 0x7fff) - 0x3fff;
  169. p -= 4;
  170. #endif
  171.  
  172. #ifdef MIEEE
  173. p = (unsigned short *)&y;
  174. e = (*p & 0x7fff) - 0x3fff;
  175. p += 5;
  176. #endif
  177.  
  178. if( e < 0 )
  179.     {
  180.     if( y < 0.0L )
  181.         return( -1.0L );
  182.     else
  183.         return( 0.0L );
  184.     }
  185.  
  186. e = (NBITS -1) - e;
  187. /* clean out 16 bits at a time */
  188. while( e >= 16 )
  189.     {
  190. #ifdef IBMPC
  191.     *p++ = 0;
  192. #endif
  193.  
  194. #ifdef MIEEE
  195.     *p-- = 0;
  196. #endif
  197.     e -= 16;
  198.     }
  199.  
  200. /* clear the remaining bits */
  201. if( e > 0 )
  202.     *p &= bmask[e];
  203.  
  204. if( (x < 0) && (y != x) )
  205.     y -= 1.0L;
  206.  
  207. return(y);
  208. }
  209.  
  210.  
  211.  
  212. long double frexpl( x, pw2 )
  213. long double x;
  214. int *pw2;
  215. {
  216. long double y;
  217. int i, k;
  218. short *q;
  219.  
  220. y = x;
  221.  
  222. #ifdef UNK
  223. mtherr( "frexp", DOMAIN );
  224. return(0.0L);
  225. #endif
  226.  
  227. /* find the exponent (power of 2) */
  228. #ifdef IBMPC
  229. q = (short *)&y + 4;
  230. i  = *q & 0x7fff;
  231. #endif
  232.  
  233. #ifdef MIEEE
  234. q = (short *)&y;
  235. i  = *q & 0x7fff;
  236. #endif
  237.  
  238.  
  239. if( i == 0 )
  240.     {
  241.     if( y == 0.0L )
  242.         {
  243.         *pw2 = 0;
  244.         return(0.0L);
  245.         }
  246. /* Number is denormal or zero */
  247. #if DENORMAL
  248. /* Handle denormal number. */
  249. do
  250.     {
  251.     y *= 2.0L;
  252.     i -= 1;
  253.     k  = *q & 0x7fff;
  254.     }
  255. while( (k == 0) && (i > -66) );
  256. i = i + k;
  257. #else
  258.     *pw2 = 0;
  259.     return(0.0L);
  260. #endif /* DENORMAL */
  261.     }
  262.  
  263. *pw2 = i - 0x3ffe;
  264. *q = 0x3ffe;
  265. return( y );
  266. }
  267.  
  268.  
  269.  
  270.  
  271.  
  272.  
  273. long double ldexpl( x, pw2 )
  274. long double x;
  275. int pw2;
  276. {
  277. long double y;
  278. unsigned short *q;
  279. long e;
  280.  
  281. #ifdef UNK
  282. mtherr( "ldexp", DOMAIN );
  283. return(0.0L);
  284. #endif
  285.  
  286. y = x;
  287. #ifdef IBMPC
  288. q = (unsigned short *)&y + 4;
  289. #endif
  290. #ifdef MIEEE
  291. q = (unsigned short *)&y;
  292. #endif
  293. while( (e = (*q & 0x7fffL)) == 0 )
  294.     {
  295. #if DENORMAL
  296.     if( y == 0.0L )
  297.         {
  298.         return( 0.0L );
  299.         }
  300. /* Input is denormal. */
  301.     if( pw2 > 0 )
  302.         {
  303.         y *= 2.0L;
  304.         pw2 -= 1;
  305.         }
  306.     if( pw2 < 0 )
  307.         {
  308.         if( pw2 < -64 )
  309.             return(0.0L);
  310.         y *= 0.5L;
  311.         pw2 += 1;
  312.         }
  313.     if( pw2 == 0 )
  314.         return(y);
  315. #else
  316.     return( 0.0L );
  317. #endif
  318.     }
  319.  
  320. e = e + pw2;
  321.  
  322. /* Handle overflow */
  323. if( e > 0x7fffL )
  324.     {
  325.     return( MAXNUML );
  326.     }
  327. *q &= 0x8000;
  328. /* Handle denormalized results */
  329. if( e < 1 )
  330.     {
  331. #if DENORMAL
  332.     if( e < -64 )
  333.         return(0.0L);
  334.  
  335. #ifdef IBMPC
  336.     *(q-1) |= 0x8000;
  337. #endif
  338. #ifdef MIEEE
  339.     *(q+2) |= 0x8000;
  340. #endif
  341.  
  342.     while( e < 1 )
  343.         {
  344.         y *= 0.5L;
  345.         e += 1;
  346.         }
  347.     e = 0;
  348. #else
  349.     return(0.0L);
  350. #endif
  351.     }
  352.  
  353. *q |= e & 0x7fff;
  354. return(y);
  355. }
  356.