home *** CD-ROM | disk | FTP | other *** search
/ Power-Programmierung / CD2.mdf / c / library / dos / math / cephes / ldouble / exp10l.c < prev    next >
Encoding:
C/C++ Source or Header  |  1992-11-17  |  3.9 KB  |  170 lines
  1. /*                            exp10l.c
  2.  *
  3.  *    Base 10 exponential function, long double precision
  4.  *      (Common antilogarithm)
  5.  *
  6.  *
  7.  *
  8.  * SYNOPSIS:
  9.  *
  10.  * long double x, y, exp10l()
  11.  *
  12.  * y = exp10l( 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.  *     1 + 2x P(x**2)/( Q(x**2) - P(x**2) )
  25.  *
  26.  * is used to approximate 10**f.
  27.  *
  28.  *
  29.  *
  30.  * ACCURACY:
  31.  *
  32.  *                      Relative error:
  33.  * arithmetic   domain     # trials      peak         rms
  34.  *    IEEE      +-4900      30000       1.0e-19     2.7e-20
  35.  *
  36.  * ERROR MESSAGES:
  37.  *
  38.  *   message         condition      value returned
  39.  * exp10l underflow    x < -MAXL10        0.0
  40.  * exp10l overflow     x > MAXL10       MAXNUM
  41.  *
  42.  * IEEE arithmetic: MAXL10 = 4932.0754489586679023819
  43.  *
  44.  */
  45.  
  46. /*
  47. Cephes Math Library Release 2.2:  January, 1991
  48. Copyright 1984, 1991 by Stephen L. Moshier
  49. Direct inquiries to 30 Frost Street, Cambridge, MA 02140
  50. */
  51.  
  52.  
  53. #include "mconf.h"
  54.  
  55. #ifdef UNK
  56. static long double P[] = {
  57.  3.1341179396892496811523E1L,
  58.  4.5618283154904699073999E3L,
  59.  1.3433113468542797218610E5L,
  60.  7.6025447914440301593592E5L,
  61. };
  62. static long double Q[] = {
  63. /* 1.0000000000000000000000E0,*/
  64.  4.7705440288425157637739E2L,
  65.  2.9732606548049614870598E4L,
  66.  4.0843697951001026189583E5L,
  67.  6.6034865026929015925608E5L,
  68. };
  69. static long double LOG102 = 3.0102999566398119521373889e-1L;
  70. static long double LOG210 = 3.3219280948873623478703L;
  71. static long double LG102A = 3.01025390625e-1L;
  72. static long double LG102B = 4.6050389811952137388947e-6L;
  73. #endif
  74.  
  75.  
  76. #ifdef IBMPC
  77. static short P[] = {
  78. 0x399a,0x7dc7,0xbc43,0xfaba,0x4003,
  79. 0xb526,0xdf32,0xa063,0x8e8e,0x400b,
  80. 0x18da,0xafa1,0xc89e,0x832e,0x4010,
  81. 0x503d,0x9352,0xe7aa,0xb99b,0x4012,
  82. };
  83. static short Q[] = {
  84. /*0x0000,0x0000,0x0000,0x8000,0x3fff,*/
  85. 0x947d,0x7855,0xf6ac,0xee86,0x4007,
  86. 0x18cf,0x7749,0x368d,0xe849,0x400d,
  87. 0x85be,0x2560,0x9f58,0xc76e,0x4011,
  88. 0x6d3c,0x80c5,0xca67,0xa137,0x4012,
  89. };
  90. static short L102[] = {0xf799,0xfbcf,0x9a84,0x9a20,0x3ffd};
  91. #define LOG102 *(long double *)L102
  92. static short L210[] = {0x8afe,0xcd1b,0x784b,0xd49a,0x4000};
  93. #define LOG210 *(long double *)L210
  94. static short L102A[] = {0x0000,0x0000,0x0000,0x9a20,0x3ffd};
  95. #define LG102A *(long double *)L102A
  96. static short L102B[] = {0x8f89,0xf798,0xfbcf,0x9a84,0x3fed};
  97. #define LG102B *(long double *)L102B
  98. #endif
  99.  
  100. #ifdef MIEEE
  101. static long P[] = {
  102. 0x40030000,0xfababc43,0x7dc7399a,
  103. 0x400b0000,0x8e8ea063,0xdf32b526,
  104. 0x40100000,0x832ec89e,0xafa118da,
  105. 0x40120000,0xb99be7aa,0x9352503d,
  106. };
  107. static long Q[] = {
  108. 0x3fff0000,0x80000000,0x00000000,
  109. 0x40070000,0xee86f6ac,0x7855947d,
  110. 0x400d0000,0xe849368d,0x774918cf,
  111. 0x40110000,0xc76e9f58,0x256085be,
  112. 0x40120000,0xa137ca67,0x80c56d3c,
  113. };
  114. static long L102[] = {0x3ffd0000,0x9a209a84,0xfbcff799};
  115. #define LOG102 *(long double *)L102
  116. static long L210[] = {0x40000000,0xd49a784b,0xcd1b8afe};
  117. #define LOG210 *(long double *)L210
  118. static long L102A[] = {0x3ffd0000,0x9a200000,0x00000000};
  119. #define LG102A *(long double *)L102A
  120. static long L102B[] = {0x3fed0000,0x9a84fbcf,0xf7988f89};
  121. #define LG102B *(long double *)L102B
  122. #endif
  123.  
  124. static long double MAXL10 = 4.9320754489586679023819e3L;
  125. extern long double MAXNUML;
  126.  
  127.  
  128.  
  129. long double exp10l(x)
  130. long double x;
  131. {
  132. long double px, xx;
  133. short n;
  134. long double floorl(), ldexpl(), polevll(), p1evll();
  135.  
  136. if( x > MAXL10 )
  137.     {
  138.     mtherr( "exp10l", OVERFLOW );
  139.     return( MAXNUML );
  140.     }
  141.  
  142. if( x < -MAXL10 )    /* Would like to use MINLOG but can't */
  143.     {
  144.     mtherr( "exp10l", UNDERFLOW );
  145.     return(0.0L);
  146.     }
  147.  
  148. /* Express 10**x = 10**g 2**n
  149.  *   = 10**g 10**( n log10(2) )
  150.  *   = 10**( g + n log10(2) )
  151.  */
  152. px = floorl( LOG210 * x + 0.5L );
  153. n = px;
  154. x -= px * LG102A;
  155. x -= px * LG102B;
  156.  
  157. /* rational approximation for exponential
  158.  * of the fractional part:
  159.  * 10**x = 1 + 2x P(x**2)/( Q(x**2) - P(x**2) )
  160.  */
  161. xx = x * x;
  162. px = x * polevll( xx, P, 3 );
  163. x =  px/( p1evll( xx, Q, 4 ) - px );
  164. x = 1.0 + ldexpl( x, 1 );
  165.  
  166. /* multiply by power of 2 */
  167. x = ldexpl( x, n );
  168. return(x);
  169. }
  170.