home *** CD-ROM | disk | FTP | other *** search
/ Power-Programmierung / CD2.mdf / c / library / dos / math / cephes / cmath / exp10.c < prev    next >
Encoding:
C/C++ Source or Header  |  1992-11-17  |  4.4 KB  |  219 lines
  1. /*                            exp10.c
  2.  *
  3.  *    Base 10 exponential function
  4.  *      (Common antilogarithm)
  5.  *
  6.  *
  7.  *
  8.  * SYNOPSIS:
  9.  *
  10.  * double x, y, exp10();
  11.  *
  12.  * y = exp10( x );
  13.  *
  14.  *
  15.  *
  16.  * DESCRIPTION:
  17.  *
  18.  * Returns 10 raised to the x power.
  19.  *
  20.  * Range reduction is accomplished by expressing the argument
  21.  * as 10**x = 2**n 10**f, with |f| < 0.5 log10(2).
  22.  * The Pade' form
  23.  *
  24.  *       2x P(x**2)/( Q(x**2) - P(x**2) )
  25.  *
  26.  * is used to approximate 10**f - 1.
  27.  *
  28.  *
  29.  *
  30.  * ACCURACY:
  31.  *
  32.  *                      Relative error:
  33.  * arithmetic   domain     # trials      peak         rms
  34.  *    DEC       -38,+38     47000       2.7e-17     6.2e-18
  35.  *    IEEE     -308,+308    30000       4.4e-16     6.2e-17
  36.  *
  37.  * ERROR MESSAGES:
  38.  *
  39.  *   message         condition      value returned
  40.  * exp10 underflow    x < -MAXL10        0.0
  41.  * exp10 overflow     x > MAXL10       MAXNUM
  42.  *
  43.  * DEC arithmetic: MAXL10 = 38.230809449325611792.
  44.  * IEEE arithmetic: MAXL10 = 308.2547155599167.
  45.  *
  46.  */
  47.  
  48. /*
  49. Cephes Math Library Release 2.1:  December, 1988
  50. Copyright 1984, 1987, 1988 by Stephen L. Moshier
  51. Direct inquiries to 30 Frost Street, Cambridge, MA 02140
  52. */
  53.  
  54.  
  55. #include "mconf.h"
  56.  
  57. #ifdef UNK
  58. static double P[] = {
  59.  7.67176747836011812180E-2,
  60.  6.08131766159694132207E0,
  61.  4.12967835799494649943E1
  62. };
  63. static double Q[] = {
  64. /* 1.00000000000000000000E0,*/
  65.  2.11303917828964245249E1,
  66.  3.58699304582497404950E1
  67. };
  68. static double LOG102 = 3.01029995663981195214e-1;
  69. static double LOG210 = 3.32192809488736234787e0;
  70. static double L102A = 3.00781250000000000000E-1;
  71. static double L102B = 2.48745663981195213739E-4
  72. static double MAXL10 = 38.230809449325611792;
  73. #endif
  74.  
  75. #ifdef DEC
  76. static short P[] = {
  77. 0037235,0017050,0000704,0007227,
  78. 0040702,0115047,0077444,0126202,
  79. 0041445,0027750,0004347,0076663
  80. };
  81. static short Q[] = {
  82. /*0040200,0000000,0000000,0000000,*/
  83. 0041251,0005412,0154331,0124200,
  84. 0041417,0075317,0006317,0164140
  85. };
  86. static short L102[] = {
  87. 0037632,0020232,0102373,0147770
  88. };
  89. #define LOG102 *(double *)L102
  90. static short L210[] = {
  91. 0040524,0115170,0045715,0015613
  92. };
  93. #define LOG210 *(double *)L210
  94. static short L102A[] = {
  95. 0037632,0000000,0000000,0000000
  96. };
  97. #define LG102A *(double *)L102A
  98. static short L102B[] = {
  99. 0035202,0065023,0167477,0157142
  100. };
  101. #define LG102B *(double *)L102B
  102. static short MXL[] = {
  103. 0041430,0166131,0047761,0154130,
  104. };
  105. #define MAXL10 ( *(double *)MXL )
  106. #endif
  107.  
  108. #ifdef IBMPC
  109. static short P[] = {
  110. 0x81d3,0x0038,0xa3c5,0x3fb3,
  111. 0x9590,0xefe4,0x5344,0x4018,
  112. 0xefb6,0x011c,0xa5fd,0x4044
  113. };
  114. static short Q[] = {
  115. /*0x0000,0x0000,0x0000,0x3ff0,*/
  116. 0x3510,0x5b1b,0x2161,0x4035,
  117. 0xfd0c,0xe199,0xef59,0x4041
  118. };
  119. static short L102[] = {
  120. 0x79ff,0x509f,0x4413,0x3fd3
  121. };
  122. #define LOG102 *(double *)L102
  123. static short L210[] = {
  124. 0xa371,0x0979,0x934f,0x400a
  125. };
  126. #define LOG210 *(double *)L210
  127. static short L102A[] = {
  128. 0x0000,0x0000,0x4000,0x3fd3
  129. };
  130. #define LG102A *(double *)L102A
  131. static short L102B[] = {
  132. 0xfbcc,0x7de7,0x4d42,0x3f30
  133. };
  134. #define LG102B *(double *)L102B
  135. static double MAXL10 = 308.2547155599167;
  136. #endif
  137.  
  138. #ifdef MIEEE
  139. static short P[] = {
  140. 0x3fb3,0xa3c5,0x0038,0x81d3,
  141. 0x4018,0x5344,0xefe4,0x9590,
  142. 0x4044,0xa5fd,0x011c,0xefb6
  143. };
  144. static short Q[] = {
  145. 0x4035,0x2161,0x5b1b,0x3510,
  146. 0x4041,0xef59,0xe199,0xfd0c
  147. };
  148. static short L102[] = {
  149. 0x3fd3,0x4413,0x509f,0x79ff
  150. };
  151. #define LOG102 *(double *)L102
  152. static short L210[] = {
  153. 0x400a,0x934f,0x0979,0xa371
  154. };
  155. #define LOG210 *(double *)L210
  156. static short L102A[] = {
  157. 0x3fd3,0x4000,0x0000,0x0000
  158. };
  159. #define LG102A *(double *)L102A
  160. static short L102B[] = {
  161. 0x3f30,0x4d42,0x7de7,0xfbcc
  162. };
  163. #define LG102B *(double *)L102B
  164. static double MAXL10 = 308.2547155599167;
  165. #endif
  166.  
  167.  
  168.  
  169. extern double MAXNUM;
  170.  
  171. double exp10(x)
  172. double x;
  173. {
  174. double px, qx, xx;
  175. short n;
  176. double floor(), ldexp(), polevl(), p1evl();
  177.  
  178. if( x > MAXL10 )
  179.     {
  180.     mtherr( "exp10", OVERFLOW );
  181.     return( MAXNUM );
  182.     }
  183.  
  184. if( x < -MAXL10 )    /* Would like to use MINLOG but can't */
  185.     {
  186.     mtherr( "exp10", UNDERFLOW );
  187.     return(0.0);
  188.     }
  189.  
  190. /* The following is necessary because range reduction blows up: */
  191. if( x == 0 )
  192.     return(1.0);
  193.  
  194. /* Express 10**x = 10**g 2**n
  195.  *   = 10**g 10**( n log10(2) )
  196.  *   = 10**( g + n log10(2) )
  197.  */
  198. px = x * LOG210;
  199. qx = floor( px + 0.5 );
  200. n = qx;
  201. x -= qx * LG102A;
  202. x -= qx * LG102B;
  203.  
  204. /* rational approximation for exponential
  205.  * of the fractional part:
  206.  * 10**x - 1  =  2x P(x**2)/( Q(x**2) - P(x**2) )
  207.  */
  208. xx = x * x;
  209. px = x * polevl( xx, P, 2 );
  210. x =  px/( p1evl( xx, Q, 2 ) - px );
  211. x = ldexp( x, 1 );
  212. x =  x + 1.0;
  213.  
  214. /* multiply by power of 2 */
  215. x = ldexp( x, n );
  216.  
  217. return(x);
  218. }
  219.