home *** CD-ROM | disk | FTP | other *** search
/ Simtel MSDOS 1992 June / SIMTEL_0692.cdr / msdos / c / cephes.arc / LOG10.C < prev    next >
Encoding:
C/C++ Source or Header  |  1987-06-20  |  4.8 KB  |  224 lines
  1. /*                            log10.c
  2.  *
  3.  *    Common logarithm
  4.  *
  5.  *
  6.  *
  7.  * SYNOPSIS:
  8.  *
  9.  * double x, y, log10();
  10.  *
  11.  * y = log10( x );
  12.  *
  13.  *
  14.  *
  15.  * DESCRIPTION:
  16.  *
  17.  * Returns logarithm to the base 10 of x.
  18.  *
  19.  * The argument is separated into its exponent and fractional
  20.  * parts.  The logarithm of the fraction is approximated by
  21.  *
  22.  *     log(1+x) = x - 0.5 x**2 + x**3 P(x)/Q(x).
  23.  *
  24.  *
  25.  *
  26.  * ACCURACY:
  27.  *
  28.  *                      Relative error:
  29.  * arithmetic   range      # trials      peak         rms
  30.  *    IEEE      0.5, 2.0    
  31.  *    IEEE      0, MAXNUM    10000      1.1e-16     4.7e-17
  32.  *    DEC       0.5, 2.0    
  33.  *    DEC       1, MAXNUM    14000      2.5e-17     6.0e-18
  34.  *
  35.  * In the tests over the interval [1, MAXNUM], the logarithms
  36.  * of the random arguments were uniformly distributed over
  37.  * [0, MAXLOG].
  38.  *
  39.  * ERROR MESSAGES:
  40.  *
  41.  * log10 singularity:  x = 0; returns MINLOG
  42.  * log10 domain:       x < 0; returns MINLOG
  43.  */
  44.  
  45. /*
  46. Cephes Math Library Release 2.0:  April, 1987
  47. Copyright 1984, 1987 by Stephen L. Moshier
  48. Direct inquiries to 30 Frost Street, Cambridge, MA 02140
  49. */
  50.  
  51. #include "mconf.h"
  52. static char fname[] = {"log10"};
  53.  
  54. /* Coefficients for log(1+x) = x - x**2/2 + x**3 P(x)/Q(x)
  55.  * 1/sqrt(2) <= x < sqrt(2)
  56.  */
  57. #ifdef UNK
  58. static double P[] = {
  59.   4.58482948458143443514E-5,
  60.   4.98531067254050724270E-1,
  61.   6.56312093769992875930E0,
  62.   2.97877425097986925891E1,
  63.   6.06127134467767258030E1,
  64.   5.67349287391754285487E1,
  65.   1.98892446572874072159E1
  66. };
  67. static double Q[] = {
  68. /* 1.00000000000000000000E0, */
  69.   1.50314182634250003249E1,
  70.   8.27410449222435217021E1,
  71.   2.20664384982121929218E2,
  72.   3.07254189979530058263E2,
  73.   2.14955586696422947765E2,
  74.   5.96677339718622216300E1
  75. };
  76. #endif
  77.  
  78. #ifdef DEC
  79. static short P[] = {
  80. 0034500,0046473,0051374,0135174,
  81. 0037777,0037566,0145712,0150321,
  82. 0040722,0002426,0031543,0123107,
  83. 0041356,0046513,0170752,0004346,
  84. 0041562,0071553,0023536,0163343,
  85. 0041542,0170221,0024316,0114216,
  86. 0041237,0016454,0046611,0104602
  87. };
  88. static short Q[] = {
  89. /*0040200,0000000,0000000,0000000,*/
  90. 0041160,0100260,0067736,0102424,
  91. 0041645,0075552,0036563,0147072,
  92. 0042134,0125025,0021132,0025320,
  93. 0042231,0120211,0046030,0103271,
  94. 0042126,0172241,0052151,0120426,
  95. 0041556,0125702,0072116,0047103
  96. };
  97. #endif
  98.  
  99. #ifdef IBMPC
  100. static short P[] = {
  101. 0x974f,0x6a5f,0x09a7,0x3f08,
  102. 0x5a1a,0xd979,0xe7ee,0x3fdf,
  103. 0x74c9,0xc66c,0x40a2,0x401a,
  104. 0x411d,0x7e3d,0xc9a9,0x403d,
  105. 0xdcdc,0x64eb,0x4e6d,0x404e,
  106. 0xd312,0x2519,0x5e12,0x404c,
  107. 0x3130,0x89b1,0xe3a5,0x4033
  108. };
  109. static short Q[] = {
  110. /*0x0000,0x0000,0x0000,0x3ff0,*/
  111. 0xd0a2,0x0dfb,0x1016,0x402e,
  112. 0x79c7,0x47ae,0xaf6d,0x4054,
  113. 0x455a,0xa44b,0x9542,0x406b,
  114. 0x10d7,0x2983,0x3411,0x4073,
  115. 0x3423,0x2a8d,0xde94,0x406a,
  116. 0xc9c8,0x4e89,0xd578,0x404d
  117. };
  118. #endif
  119.  
  120. #ifdef MIEEE
  121. static short P[] = {
  122. 0x3f08,0x09a7,0x6a5f,0x974f,
  123. 0x3fdf,0xe7ee,0xd979,0x5a1a,
  124. 0x401a,0x40a2,0xc66c,0x74c9,
  125. 0x403d,0xc9a9,0x7e3d,0x411d,
  126. 0x404e,0x4e6d,0x64eb,0xdcdc,
  127. 0x404c,0x5e12,0x2519,0xd312,
  128. 0x4033,0xe3a5,0x89b1,0x3130
  129. };
  130. static short Q[] = {
  131. 0x402e,0x1016,0x0dfb,0xd0a2,
  132. 0x4054,0xaf6d,0x47ae,0x79c7,
  133. 0x406b,0x9542,0xa44b,0x455a,
  134. 0x4073,0x3411,0x2983,0x10d7,
  135. 0x406a,0xde94,0x2a8d,0x3423,
  136. 0x404d,0xd578,0x4e89,0xc9c8
  137. };
  138. #endif
  139.  
  140. #define SQRTH 0.70710678118654752440
  141. #define L102A 3.0078125E-1
  142. #define L102B 2.48745663981195213739E-4
  143. #define L10EA 4.3359375E-1
  144. #define L10EB 7.00731903251827651129E-4
  145.  
  146. extern double LOGE2, SQRT2, MAXLOG, MINLOG;
  147.  
  148. double log10(x)
  149. double x;
  150. {
  151. short e, *q;
  152. double w, y, z;
  153. double frexp(), ldexp(), polevl(), p1evl();
  154.  
  155. /* Test for domain */
  156. if( x <= 0.0 )
  157.     {
  158.     if( x == 0.0 )
  159.         mtherr( fname, SING );
  160.     else
  161.         mtherr( fname, DOMAIN );
  162.     return( MINLOG );
  163.     }
  164.  
  165. /* separate mantissa from exponent */
  166.  
  167. #ifdef DEC
  168. q = (short *)&x;
  169. e = *q;            /* short containing exponent */
  170. e = ((e >> 7) & 0377) - 0200;    /* the exponent */
  171. *q &= 0177;    /* strip exponent from x */
  172. *q |= 040000;    /* x now between 0.5 and 1 */
  173. #endif
  174.  
  175. #ifdef IBMPC
  176. q = (short *)&x;
  177. q += 3;
  178. e = *q;
  179. e = ((e >> 4) & 0x0fff) - 0x3fe;
  180. *q &= 0x0f;
  181. *q |= 0x3fe0;
  182. #endif
  183.  
  184. /* Equivalent C language standard library function: */
  185. #ifdef UNK
  186. x = frexp( x, &e );
  187. #endif
  188.  
  189. #ifdef MIEEE
  190. x = frexp( x, &e );
  191. #endif
  192.  
  193. /* logarithm using log(1+x) = x - .5x**2 + x**3 P(x)/Q(x) */
  194.  
  195. if( x < SQRTH )
  196.     {
  197.     e -= 1;
  198.     x = ldexp( x, 1 ) - 1.0; /*  2x - 1  */
  199.     }    
  200. else
  201.     {
  202.     x = x - 1.0;
  203.     }
  204.  
  205.  
  206. /* rational form */
  207. z = x*x;
  208. y = x * ( z * polevl( x, P, 6 ) / p1evl( x, Q, 6 ) );
  209. y = y - ldexp( z, -1 );   /*  y - 0.5 * x**2  */
  210.  
  211. /* multiply log of fraction by log10(e)
  212.  * and base 2 exponent by log10(2)
  213.  */
  214. z = (x + y) * L10EB;  /* accumulate terms in order of size */
  215. z += y * L10EA;
  216. z += x * L10EA;
  217. w = e;
  218. z += w * L102B;
  219. z += w * L102A;
  220.  
  221.  
  222. return( z );
  223. }
  224.