home *** CD-ROM | disk | FTP | other *** search
/ Simtel MSDOS 1992 June / SIMTEL_0692.cdr / msdos / c / cephes.arc / LOG.C < prev    next >
Encoding:
C/C++ Source or Header  |  1987-06-20  |  7.0 KB  |  318 lines
  1. /*                            log.c
  2.  *
  3.  *    Natural logarithm
  4.  *
  5.  *
  6.  *
  7.  * SYNOPSIS:
  8.  *
  9.  * double x, y, log();
  10.  *
  11.  * y = log( x );
  12.  *
  13.  *
  14.  *
  15.  * DESCRIPTION:
  16.  *
  17.  * Returns logarithm to the base e (2.718...) of x.
  18.  *
  19.  * The argument is separated into its exponent and fractional
  20.  * parts.  If the exponent is between -1 and +1, the logarithm
  21.  * of the fraction is approximated by
  22.  *
  23.  *     log(1+x) = x - 0.5 x**2 + x**3 P(x)/Q(x).
  24.  *
  25.  * Otherwise, setting  z = 2(x-1)/x+1),
  26.  * 
  27.  *     log(x) = z + z**3 P(z)/Q(z).
  28.  *
  29.  *
  30.  *
  31.  * ACCURACY:
  32.  *
  33.  *                      Relative error:
  34.  * arithmetic   range      # trials      peak         rms
  35.  *    IEEE      0.5, 2.0    35000       1.8e-16     5.5e-17
  36.  *    IEEE      1, MAXNUM   10000       1.4e-16     4.8e-17
  37.  *    DEC       0.5, 2.0    20000       2.0e-17     7.0e-18
  38.  *    DEC       1, MAXNUM   17700       1.9e-17     6.1e-18
  39.  *
  40.  * In the tests over the interval [1, MAXNUM], the logarithms
  41.  * of the random arguments were uniformly distributed over
  42.  * [0, MAXLOG].
  43.  *
  44.  * ERROR MESSAGES:
  45.  *
  46.  * log singularity:  x = 0; returns MINLOG
  47.  * log domain:       x < 0; returns MINLOG
  48.  */
  49.  
  50. /*
  51. Cephes Math Library Release 2.0:  April, 1987
  52. Copyright 1984, 1987 by Stephen L. Moshier
  53. Direct inquiries to 30 Frost Street, Cambridge, MA 02140
  54. */
  55.  
  56. #include "mconf.h"
  57. static char fname[] = {"log"};
  58.  
  59. /* Coefficients for log(1+x) = x - x**2/2 + x**3 P(x)/Q(x)
  60.  * 1/sqrt(2) <= x < sqrt(2)
  61.  */
  62. #ifdef UNK
  63. static double P[] = {
  64.   4.58482948458143443514E-5,
  65.   4.98531067254050724270E-1,
  66.   6.56312093769992875930E0,
  67.   2.97877425097986925891E1,
  68.   6.06127134467767258030E1,
  69.   5.67349287391754285487E1,
  70.   1.98892446572874072159E1
  71. };
  72. static double Q[] = {
  73. /* 1.00000000000000000000E0, */
  74.   1.50314182634250003249E1,
  75.   8.27410449222435217021E1,
  76.   2.20664384982121929218E2,
  77.   3.07254189979530058263E2,
  78.   2.14955586696422947765E2,
  79.   5.96677339718622216300E1
  80. };
  81. #endif
  82.  
  83. #ifdef DEC
  84. static short P[] = {
  85. 0034500,0046473,0051374,0135174,
  86. 0037777,0037566,0145712,0150321,
  87. 0040722,0002426,0031543,0123107,
  88. 0041356,0046513,0170752,0004346,
  89. 0041562,0071553,0023536,0163343,
  90. 0041542,0170221,0024316,0114216,
  91. 0041237,0016454,0046611,0104602
  92. };
  93. static short Q[] = {
  94. /*0040200,0000000,0000000,0000000,*/
  95. 0041160,0100260,0067736,0102424,
  96. 0041645,0075552,0036563,0147072,
  97. 0042134,0125025,0021132,0025320,
  98. 0042231,0120211,0046030,0103271,
  99. 0042126,0172241,0052151,0120426,
  100. 0041556,0125702,0072116,0047103
  101. };
  102. #endif
  103.  
  104. #ifdef IBMPC
  105. static short P[] = {
  106. 0x974f,0x6a5f,0x09a7,0x3f08,
  107. 0x5a1a,0xd979,0xe7ee,0x3fdf,
  108. 0x74c9,0xc66c,0x40a2,0x401a,
  109. 0x411d,0x7e3d,0xc9a9,0x403d,
  110. 0xdcdc,0x64eb,0x4e6d,0x404e,
  111. 0xd312,0x2519,0x5e12,0x404c,
  112. 0x3130,0x89b1,0xe3a5,0x4033
  113. };
  114. static short Q[] = {
  115. /*0x0000,0x0000,0x0000,0x3ff0,*/
  116. 0xd0a2,0x0dfb,0x1016,0x402e,
  117. 0x79c7,0x47ae,0xaf6d,0x4054,
  118. 0x455a,0xa44b,0x9542,0x406b,
  119. 0x10d7,0x2983,0x3411,0x4073,
  120. 0x3423,0x2a8d,0xde94,0x406a,
  121. 0xc9c8,0x4e89,0xd578,0x404d
  122. };
  123. #endif
  124.  
  125. #ifdef MIEEE
  126. static short P[] = {
  127. 0x3f08,0x09a7,0x6a5f,0x974f,
  128. 0x3fdf,0xe7ee,0xd979,0x5a1a,
  129. 0x401a,0x40a2,0xc66c,0x74c9,
  130. 0x403d,0xc9a9,0x7e3d,0x411d,
  131. 0x404e,0x4e6d,0x64eb,0xdcdc,
  132. 0x404c,0x5e12,0x2519,0xd312,
  133. 0x4033,0xe3a5,0x89b1,0x3130
  134. };
  135. static short Q[] = {
  136. 0x402e,0x1016,0x0dfb,0xd0a2,
  137. 0x4054,0xaf6d,0x47ae,0x79c7,
  138. 0x406b,0x9542,0xa44b,0x455a,
  139. 0x4073,0x3411,0x2983,0x10d7,
  140. 0x406a,0xde94,0x2a8d,0x3423,
  141. 0x404d,0xd578,0x4e89,0xc9c8
  142. };
  143. #endif
  144.  
  145. /* Coefficients for log(x) = z + z**3 P(z)/Q(z),
  146.  * where z = 2(x-1)/(x+1)
  147.  * 1/sqrt(2) <= x < sqrt(2)
  148.  */
  149.  
  150. #ifdef UNK
  151. static double R[] = {
  152. -7.90003033195003570794E-1,
  153.  1.64470394105188802674E1,
  154. -6.44977287217605054794E1,
  155. };
  156. static double S[] = {
  157. /* 1.00000000000000000000E0, */
  158. -3.57676812896550927595E1,
  159.  3.13460384625398562220E2,
  160. -7.73972744661126066616E2,
  161. };
  162. #endif
  163.  
  164. #ifdef DEC
  165. static short R[] = {
  166. 0140112,0036643,0103520,0033473,
  167. 0041203,0111611,0063001,0116432,
  168. 0141600,0177326,0046214,0075606,
  169. };
  170. static short S[] = {
  171.  
  172. 0141417,0011033,0005503,0043546,
  173. 0042234,0135355,0161046,0152602,
  174. 0142501,0077101,0071322,0134511,
  175. };
  176. #endif
  177.  
  178. #ifdef IBMPC
  179. static short R[] = {
  180. 0x06e7,0x70ea,0x47b4,0xbfe9,
  181. 0x33a3,0x2cc0,0x7271,0x4030,
  182. 0x8f71,0xc991,0x1fda,0xc050
  183. };
  184. static short S[] = {
  185.  
  186. 0x68ed,0x6168,0xe243,0xc041,
  187. 0xdab0,0xbc44,0x975d,0x4073,
  188. 0x5729,0x2e5a,0x2fc8,0xc088
  189. };
  190. #endif
  191.  
  192. #ifdef MIEEE
  193. static short R[] = {
  194. 0xbfe9,0x47b4,0x70ea,0x06e7,
  195. 0x4030,0x7271,0x2cc0,0x33a3,
  196. 0xc050,0x1fda,0xc991,0x8f71
  197. };
  198. static short S[] = {
  199. 0xc041,0xe243,0x6168,0x68ed,
  200. 0x4073,0x975d,0xbc44,0xdab0,
  201. 0xc088,0x2fc8,0x2e5a,0x5729
  202. };
  203. #endif
  204.  
  205.  
  206. #define SQRTH 0.70710678118654752440
  207. extern double LOGE2, SQRT2, MAXLOG, MINLOG;
  208.  
  209. double log(x)
  210. double x;
  211. {
  212. short e, *q;
  213. double y, z;
  214. double frexp(), ldexp(), polevl(), p1evl();
  215.  
  216. /* Test for domain */
  217. if( x <= 0.0 )
  218.     {
  219.     if( x == 0.0 )
  220.         mtherr( fname, SING );
  221.     else
  222.         mtherr( fname, DOMAIN );
  223.     return( MINLOG );
  224.     }
  225.  
  226. /* separate mantissa from exponent */
  227.  
  228. #ifdef DEC
  229. q = (short *)&x;
  230. e = *q;            /* short containing exponent */
  231. e = ((e >> 7) & 0377) - 0200;    /* the exponent */
  232. *q &= 0177;    /* strip exponent from x */
  233. *q |= 040000;    /* x now between 0.5 and 1 */
  234. #endif
  235.  
  236. #ifdef IBMPC
  237. q = (short *)&x;
  238. q += 3;
  239. e = *q;
  240. e = ((e >> 4) & 0x0fff) - 0x3fe;
  241. *q &= 0x0f;
  242. *q |= 0x3fe0;
  243. #endif
  244.  
  245. /* Equivalent C language standard library function: */
  246. #ifdef UNK
  247. x = frexp( x, &e );
  248. #endif
  249.  
  250. #ifdef MIEEE
  251. x = frexp( x, &e );
  252. #endif
  253.  
  254.  
  255.  
  256. /* logarithm using log(x) = z + z**3 P(z)/Q(z),
  257.  * where z = 2(x-1)/x+1)
  258.  */
  259.  
  260. if( (e > 2) || (e < -2) )
  261. {
  262. if( x < SQRTH )
  263.     { /* 2( 2x-1 )/( 2x+1 ) */
  264.     e -= 1;
  265.     z = x - 0.5;
  266.     y = 0.5 * z + 0.5;
  267.     }    
  268. else
  269.     { /*  2 (x-1)/(x+1)   */
  270.     z = x - 0.5;
  271.     z -= 0.5;
  272.     y = 0.5 * x  + 0.5;
  273.     }
  274.  
  275. x = z / y;
  276.  
  277.  
  278. /* rational form */
  279. z = x*x;
  280. z = x + x * ( z * polevl( z, R, 2 ) / p1evl( z, S, 3 ) );
  281.  
  282. goto ldone;
  283. }
  284.  
  285.  
  286.  
  287. /* logarithm using log(1+x) = x - .5x**2 + x**3 P(x)/Q(x) */
  288.  
  289. if( x < SQRTH )
  290.     {
  291.     e -= 1;
  292.     x = ldexp( x, 1 ) - 1.0; /*  2x - 1  */
  293.     }    
  294. else
  295.     {
  296.     x = x - 1.0;
  297.     }
  298.  
  299.  
  300. /* rational form */
  301. z = x*x;
  302. y = x * ( z * polevl( x, P, 6 ) / p1evl( x, Q, 6 ) );
  303. y = y - ldexp( z, -1 );   /*  y - 0.5 * z  */
  304. z = x + y;
  305.  
  306. ldone:
  307.  
  308. /* recombine with exponent term */
  309. if( e != 0 )
  310.     {
  311.     y = e;
  312.     z = z - y * 2.121944400546905827679e-4;
  313.     z = z + y * 0.693359375;
  314.     }
  315.  
  316. return( z );
  317. }
  318.