polevll.c 4.7 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182
  1. /* polevll.c
  2. * p1evll.c
  3. *
  4. * Evaluate polynomial
  5. *
  6. *
  7. *
  8. * SYNOPSIS:
  9. *
  10. * int N;
  11. * long double x, y, coef[N+1], polevl[];
  12. *
  13. * y = polevll( x, coef, N );
  14. *
  15. *
  16. *
  17. * DESCRIPTION:
  18. *
  19. * Evaluates polynomial of degree N:
  20. *
  21. * 2 N
  22. * y = C + C x + C x +...+ C x
  23. * 0 1 2 N
  24. *
  25. * Coefficients are stored in reverse order:
  26. *
  27. * coef[0] = C , ..., coef[N] = C .
  28. * N 0
  29. *
  30. * The function p1evll() assumes that coef[N] = 1.0 and is
  31. * omitted from the array. Its calling arguments are
  32. * otherwise the same as polevll().
  33. *
  34. * This module also contains the following globally declared constants:
  35. * MAXNUML = 1.189731495357231765021263853E4932L;
  36. * MACHEPL = 5.42101086242752217003726400434970855712890625E-20L;
  37. * MAXLOGL = 1.1356523406294143949492E4L;
  38. * MINLOGL = -1.1355137111933024058873E4L;
  39. * LOGE2L = 6.9314718055994530941723E-1L;
  40. * LOG2EL = 1.4426950408889634073599E0L;
  41. * PIL = 3.1415926535897932384626L;
  42. * PIO2L = 1.5707963267948966192313L;
  43. * PIO4L = 7.8539816339744830961566E-1L;
  44. *
  45. * SPEED:
  46. *
  47. * In the interest of speed, there are no checks for out
  48. * of bounds arithmetic. This routine is used by most of
  49. * the functions in the library. Depending on available
  50. * equipment features, the user may wish to rewrite the
  51. * program in microcode or assembly language.
  52. *
  53. */
  54. /*
  55. Cephes Math Library Release 2.2: July, 1992
  56. Copyright 1984, 1987, 1988, 1992 by Stephen L. Moshier
  57. Direct inquiries to 30 Frost Street, Cambridge, MA 02140
  58. */
  59. #include <math.h>
  60. #if UNK
  61. /* almost 2^16384 */
  62. long double MAXNUML = 1.189731495357231765021263853E4932L;
  63. /* 2^-64 */
  64. long double MACHEPL = 5.42101086242752217003726400434970855712890625E-20L;
  65. /* log( MAXNUML ) */
  66. long double MAXLOGL = 1.1356523406294143949492E4L;
  67. #ifdef DENORMAL
  68. /* log(smallest denormal number = 2^-16446) */
  69. long double MINLOGL = -1.13994985314888605586758E4L;
  70. #else
  71. /* log( underflow threshold = 2^(-16382) ) */
  72. long double MINLOGL = -1.1355137111933024058873E4L;
  73. #endif
  74. long double LOGE2L = 6.9314718055994530941723E-1L;
  75. long double LOG2EL = 1.4426950408889634073599E0L;
  76. long double PIL = 3.1415926535897932384626L;
  77. long double PIO2L = 1.5707963267948966192313L;
  78. long double PIO4L = 7.8539816339744830961566E-1L;
  79. #ifdef INFINITIES
  80. long double NANL = 0.0L / 0.0L;
  81. long double INFINITYL = 1.0L / 0.0L;
  82. #else
  83. long double INFINITYL = 1.189731495357231765021263853E4932L;
  84. long double NANL = 0.0L;
  85. #endif
  86. #endif
  87. #if IBMPC
  88. short MAXNUML[] = {0xffff,0xffff,0xffff,0xffff,0x7ffe, XPD};
  89. short MAXLOGL[] = {0x79ab,0xd1cf,0x17f7,0xb172,0x400c, XPD};
  90. #ifdef INFINITIES
  91. short INFINITYL[] = {0,0,0,0x8000,0x7fff, XPD};
  92. short NANL[] = {0,0,0,0xc000,0x7fff, XPD};
  93. #else
  94. short INFINITYL[] = {0xffff,0xffff,0xffff,0xffff,0x7ffe, XPD};
  95. long double NANL = 0.0L;
  96. #endif
  97. #ifdef DENORMAL
  98. short MINLOGL[] = {0xbaaa,0x09e2,0xfe7f,0xb21d,0xc00c, XPD};
  99. #else
  100. short MINLOGL[] = {0xeb2f,0x1210,0x8c67,0xb16c,0xc00c, XPD};
  101. #endif
  102. short MACHEPL[] = {0x0000,0x0000,0x0000,0x8000,0x3fbf, XPD};
  103. short LOGE2L[] = {0x79ac,0xd1cf,0x17f7,0xb172,0x3ffe, XPD};
  104. short LOG2EL[] = {0xf0bc,0x5c17,0x3b29,0xb8aa,0x3fff, XPD};
  105. short PIL[] = {0xc235,0x2168,0xdaa2,0xc90f,0x4000, XPD};
  106. short PIO2L[] = {0xc235,0x2168,0xdaa2,0xc90f,0x3fff, XPD};
  107. short PIO4L[] = {0xc235,0x2168,0xdaa2,0xc90f,0x3ffe, XPD};
  108. #endif
  109. #if MIEEE
  110. long MAXNUML[] = {0x7ffe0000,0xffffffff,0xffffffff};
  111. long MAXLOGL[] = {0x400c0000,0xb17217f7,0xd1cf79ab};
  112. #ifdef INFINITIES
  113. long INFINITY[] = {0x7fff0000,0x80000000,0x00000000};
  114. long NANL[] = {0x7fff0000,0xffffffff,0xffffffff};
  115. #else
  116. long INFINITYL[] = {0x7ffe0000,0xffffffff,0xffffffff};
  117. long double NANL = 0.0L;
  118. #endif
  119. #ifdef DENORMAL
  120. long MINLOGL[] = {0xc00c0000,0xb21dfe7f,0x09e2baaa};
  121. #else
  122. long MINLOGL[] = {0xc00c0000,0xb16c8c67,0x1210eb2f};
  123. #endif
  124. long MACHEPL[] = {0x3fbf0000,0x80000000,0x00000000};
  125. long LOGE2L[] = {0x3ffe0000,0xb17217f7,0xd1cf79ac};
  126. long LOG2EL[] = {0x3fff0000,0xb8aa3b29,0x5c17f0bc};
  127. long PIL[] = {0x40000000,0xc90fdaa2,0x2168c235};
  128. long PIO2L[] = {0x3fff0000,0xc90fdaa2,0x2168c235};
  129. long PIO4L[] = {0x3ffe0000,0xc90fdaa2,0x2168c235};
  130. #endif
  131. #ifdef MINUSZERO
  132. long double NEGZEROL = -0.0L;
  133. #else
  134. long double NEGZEROL = 0.0L;
  135. #endif
  136. /* Polynomial evaluator:
  137. * P[0] x^n + P[1] x^(n-1) + ... + P[n]
  138. */
  139. long double polevll( x, p, n )
  140. long double x;
  141. void *p;
  142. int n;
  143. {
  144. register long double y;
  145. register long double *P = (long double *)p;
  146. y = *P++;
  147. do
  148. {
  149. y = y * x + *P++;
  150. }
  151. while( --n );
  152. return(y);
  153. }
  154. /* Polynomial evaluator:
  155. * x^n + P[0] x^(n-1) + P[1] x^(n-2) + ... + P[n]
  156. */
  157. long double p1evll( x, p, n )
  158. long double x;
  159. void *p;
  160. int n;
  161. {
  162. register long double y;
  163. register long double *P = (long double *)p;
  164. n -= 1;
  165. y = x + *P++;
  166. do
  167. {
  168. y = y * x + *P++;
  169. }
  170. while( --n );
  171. return( y );
  172. }