logl.c 6.5 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292
  1. /* logl.c
  2. *
  3. * Natural logarithm, long double precision
  4. *
  5. *
  6. *
  7. * SYNOPSIS:
  8. *
  9. * long double x, y, logl();
  10. *
  11. * y = logl( x );
  12. *
  13. *
  14. *
  15. * DESCRIPTION:
  16. *
  17. * Returns the base e (2.718...) logarithm 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 domain # trials peak rms
  35. * IEEE 0.5, 2.0 150000 8.71e-20 2.75e-20
  36. * IEEE exp(+-10000) 100000 5.39e-20 2.34e-20
  37. *
  38. * In the tests over the interval exp(+-10000), the logarithms
  39. * of the random arguments were uniformly distributed over
  40. * [-10000, +10000].
  41. *
  42. * ERROR MESSAGES:
  43. *
  44. * log singularity: x = 0; returns -INFINITYL
  45. * log domain: x < 0; returns NANL
  46. */
  47. /*
  48. Cephes Math Library Release 2.7: May, 1998
  49. Copyright 1984, 1990, 1998 by Stephen L. Moshier
  50. */
  51. #include <math.h>
  52. /* Coefficients for log(1+x) = x - x**2/2 + x**3 P(x)/Q(x)
  53. * 1/sqrt(2) <= x < sqrt(2)
  54. * Theoretical peak relative error = 2.32e-20
  55. */
  56. #ifdef UNK
  57. static long double P[] = {
  58. 4.5270000862445199635215E-5L,
  59. 4.9854102823193375972212E-1L,
  60. 6.5787325942061044846969E0L,
  61. 2.9911919328553073277375E1L,
  62. 6.0949667980987787057556E1L,
  63. 5.7112963590585538103336E1L,
  64. 2.0039553499201281259648E1L,
  65. };
  66. static long double Q[] = {
  67. /* 1.0000000000000000000000E0,*/
  68. 1.5062909083469192043167E1L,
  69. 8.3047565967967209469434E1L,
  70. 2.2176239823732856465394E2L,
  71. 3.0909872225312059774938E2L,
  72. 2.1642788614495947685003E2L,
  73. 6.0118660497603843919306E1L,
  74. };
  75. #endif
  76. #ifdef IBMPC
  77. static short P[] = {
  78. 0x51b9,0x9cae,0x4b15,0xbde0,0x3ff0, XPD
  79. 0x19cf,0xf0d4,0xc507,0xff40,0x3ffd, XPD
  80. 0x9942,0xa7d2,0xfa37,0xd284,0x4001, XPD
  81. 0x4add,0x65ce,0x9c5c,0xef4b,0x4003, XPD
  82. 0x8445,0x619a,0x75c3,0xf3cc,0x4004, XPD
  83. 0x81ab,0x3cd0,0xacba,0xe473,0x4004, XPD
  84. 0x4cbf,0xcc18,0x016c,0xa051,0x4003, XPD
  85. };
  86. static short Q[] = {
  87. /*0x0000,0x0000,0x0000,0x8000,0x3fff,*/
  88. 0xb8b7,0x81f1,0xacf4,0xf101,0x4002, XPD
  89. 0xbc31,0x09a4,0x5a91,0xa618,0x4005, XPD
  90. 0xaeec,0xe7da,0x2c87,0xddc3,0x4006, XPD
  91. 0x2bde,0x4845,0xa2ee,0x9a8c,0x4007, XPD
  92. 0x3120,0x4703,0x89f2,0xd86d,0x4006, XPD
  93. 0x7347,0x3224,0x8223,0xf079,0x4004, XPD
  94. };
  95. #endif
  96. #ifdef MIEEE
  97. static long P[] = {
  98. 0x3ff00000,0xbde04b15,0x9cae51b9,
  99. 0x3ffd0000,0xff40c507,0xf0d419cf,
  100. 0x40010000,0xd284fa37,0xa7d29942,
  101. 0x40030000,0xef4b9c5c,0x65ce4add,
  102. 0x40040000,0xf3cc75c3,0x619a8445,
  103. 0x40040000,0xe473acba,0x3cd081ab,
  104. 0x40030000,0xa051016c,0xcc184cbf,
  105. };
  106. static long Q[] = {
  107. /*0x3fff0000,0x80000000,0x00000000,*/
  108. 0x40020000,0xf101acf4,0x81f1b8b7,
  109. 0x40050000,0xa6185a91,0x09a4bc31,
  110. 0x40060000,0xddc32c87,0xe7daaeec,
  111. 0x40070000,0x9a8ca2ee,0x48452bde,
  112. 0x40060000,0xd86d89f2,0x47033120,
  113. 0x40040000,0xf0798223,0x32247347,
  114. };
  115. #endif
  116. /* Coefficients for log(x) = z + z^3 P(z^2)/Q(z^2),
  117. * where z = 2(x-1)/(x+1)
  118. * 1/sqrt(2) <= x < sqrt(2)
  119. * Theoretical peak relative error = 6.16e-22
  120. */
  121. #ifdef UNK
  122. static long double R[4] = {
  123. 1.9757429581415468984296E-3L,
  124. -7.1990767473014147232598E-1L,
  125. 1.0777257190312272158094E1L,
  126. -3.5717684488096787370998E1L,
  127. };
  128. static long double S[4] = {
  129. /* 1.00000000000000000000E0L,*/
  130. -2.6201045551331104417768E1L,
  131. 1.9361891836232102174846E2L,
  132. -4.2861221385716144629696E2L,
  133. };
  134. static long double C1 = 6.9314575195312500000000E-1L;
  135. static long double C2 = 1.4286068203094172321215E-6L;
  136. #endif
  137. #ifdef IBMPC
  138. static short R[] = {
  139. 0x6ef4,0xf922,0x7763,0x817b,0x3ff6, XPD
  140. 0x15fd,0x1af9,0xde8f,0xb84b,0xbffe, XPD
  141. 0x8b96,0x4f8d,0xa53c,0xac6f,0x4002, XPD
  142. 0x8932,0xb4e3,0xe8ae,0x8ede,0xc004, XPD
  143. };
  144. static short S[] = {
  145. /*0x0000,0x0000,0x0000,0x8000,0x3fff,*/
  146. 0x7ce4,0x1fc9,0xbdc5,0xd19b,0xc003, XPD
  147. 0x0af3,0x0d10,0x716f,0xc19e,0x4006, XPD
  148. 0x4d7d,0x0f55,0x5d06,0xd64e,0xc007, XPD
  149. };
  150. static short sc1[] = {0x0000,0x0000,0x0000,0xb172,0x3ffe, XPD};
  151. #define C1 (*(long double *)sc1)
  152. static short sc2[] = {0x4f1e,0xcd5e,0x8e7b,0xbfbe,0x3feb, XPD};
  153. #define C2 (*(long double *)sc2)
  154. #endif
  155. #ifdef MIEEE
  156. static long R[12] = {
  157. 0x3ff60000,0x817b7763,0xf9226ef4,
  158. 0xbffe0000,0xb84bde8f,0x1af915fd,
  159. 0x40020000,0xac6fa53c,0x4f8d8b96,
  160. 0xc0040000,0x8edee8ae,0xb4e38932,
  161. };
  162. static long S[9] = {
  163. /*0x3fff0000,0x80000000,0x00000000,*/
  164. 0xc0030000,0xd19bbdc5,0x1fc97ce4,
  165. 0x40060000,0xc19e716f,0x0d100af3,
  166. 0xc0070000,0xd64e5d06,0x0f554d7d,
  167. };
  168. static long sc1[] = {0x3ffe0000,0xb1720000,0x00000000};
  169. #define C1 (*(long double *)sc1)
  170. static long sc2[] = {0x3feb0000,0xbfbe8e7b,0xcd5e4f1e};
  171. #define C2 (*(long double *)sc2)
  172. #endif
  173. #define SQRTH 0.70710678118654752440L
  174. extern long double MINLOGL;
  175. #ifdef ANSIPROT
  176. extern long double frexpl ( long double, int * );
  177. extern long double ldexpl ( long double, int );
  178. extern long double polevll ( long double, void *, int );
  179. extern long double p1evll ( long double, void *, int );
  180. extern int isnanl ( long double );
  181. #else
  182. long double frexpl(), ldexpl(), polevll(), p1evll(), isnanl();
  183. #endif
  184. #ifdef INFINITIES
  185. extern long double INFINITYL;
  186. #endif
  187. #ifdef NANS
  188. extern long double NANL;
  189. #endif
  190. long double logl(x)
  191. long double x;
  192. {
  193. long double y, z;
  194. int e;
  195. #ifdef NANS
  196. if( isnanl(x) )
  197. return(x);
  198. #endif
  199. #ifdef INFINITIES
  200. if( x == INFINITYL )
  201. return(x);
  202. #endif
  203. /* Test for domain */
  204. if( x <= 0.0L )
  205. {
  206. if( x == 0.0L )
  207. {
  208. #ifdef INFINITIES
  209. return( -INFINITYL );
  210. #else
  211. mtherr( "logl", SING );
  212. return( MINLOGL );
  213. #endif
  214. }
  215. else
  216. {
  217. #ifdef NANS
  218. return( NANL );
  219. #else
  220. mtherr( "logl", DOMAIN );
  221. return( MINLOGL );
  222. #endif
  223. }
  224. }
  225. /* separate mantissa from exponent */
  226. /* Note, frexp is used so that denormal numbers
  227. * will be handled properly.
  228. */
  229. x = frexpl( x, &e );
  230. /* logarithm using log(x) = z + z**3 P(z)/Q(z),
  231. * where z = 2(x-1)/x+1)
  232. */
  233. if( (e > 2) || (e < -2) )
  234. {
  235. if( x < SQRTH )
  236. { /* 2( 2x-1 )/( 2x+1 ) */
  237. e -= 1;
  238. z = x - 0.5L;
  239. y = 0.5L * z + 0.5L;
  240. }
  241. else
  242. { /* 2 (x-1)/(x+1) */
  243. z = x - 0.5L;
  244. z -= 0.5L;
  245. y = 0.5L * x + 0.5L;
  246. }
  247. x = z / y;
  248. z = x*x;
  249. z = x * ( z * polevll( z, R, 3 ) / p1evll( z, S, 3 ) );
  250. z = z + e * C2;
  251. z = z + x;
  252. z = z + e * C1;
  253. return( z );
  254. }
  255. /* logarithm using log(1+x) = x - .5x**2 + x**3 P(x)/Q(x) */
  256. if( x < SQRTH )
  257. {
  258. e -= 1;
  259. x = ldexpl( x, 1 ) - 1.0L; /* 2x - 1 */
  260. }
  261. else
  262. {
  263. x = x - 1.0L;
  264. }
  265. z = x*x;
  266. y = x * ( z * polevll( x, P, 6 ) / p1evll( x, Q, 6 ) );
  267. y = y + e * C2;
  268. z = y - ldexpl( z, -1 ); /* y - 0.5 * z */
  269. /* Note, the sum of above terms does not exceed x/4,
  270. * so it contributes at most about 1/4 lsb to the error.
  271. */
  272. z = z + x;
  273. z = z + e * C1; /* This sum has an error of 1/2 lsb. */
  274. return( z );
  275. }