home *** CD-ROM | disk | FTP | other *** search
/ Power-Programmierung / CD2.mdf / c / library / dos / math / cephes / remes / remesf.c < prev    next >
Encoding:
C/C++ Source or Header  |  1992-11-17  |  5.9 KB  |  380 lines
  1. /* remesf.c */
  2.  
  3. #include "remes.h"
  4.  
  5.  
  6. /* The flag bits for type of approximation: */
  7. /* PXSQ | XPX | X2PX | SQL | SQH | PADE | CW | ZER | SPECIAL | PXCU */
  8. int config =   ZER ;
  9.  
  10. /* Insert function name and formulas for printout */
  11. char funnam[] = {
  12.  "log(1+x) "
  13. };
  14.  
  15. char znam[] = { "x" };
  16. double sqrt(), log(), exp(), pow(), sin(), cos();
  17. double j0(), y0(), atan(), floor();
  18. double j1(), y1();
  19. double exp10(), log10();
  20. double erfc();
  21. double lx();
  22.  
  23. /* This subroutine computes the rational form P(x)/Q(x) */
  24. /* using coefficients from the solution vector param[]. */
  25.  
  26. double approx(x)
  27. double x;
  28. {
  29. double gx, z, yn, yd, q;
  30. double gofx(), speci();
  31. int i;
  32.  
  33. gx = gofx(x);
  34. if( config & PXCU )
  35.     z = gx * gx * gx;
  36. else if( config & PXSQ )
  37.     z = gx * gx;
  38. else
  39.     z = gx;
  40.  
  41. /* Highest order numerator coefficient */
  42. yn = param[n];
  43.  
  44. /* Work backwards toward the constant term. */
  45. for( i=n-1; i>=0; i-- )
  46.     yn = z * yn  +  param[i];
  47.     
  48. if( d > 0 )
  49.     {
  50. /* Highest degree coefficient = 1.0 */
  51.     yd = z + param[n+d];
  52.  
  53.     for( i=n+d-1; i>n; i-- )
  54.         yd = z * yd  +  param[i];
  55.     }
  56. else
  57. /* There is no denominator. */
  58.     yd = 1.0;
  59.  
  60. if( config & XPX )
  61.     yn = yn * gx;
  62. if( config & X2PX )
  63.     yn = yn * gx * gx;
  64. if( config & PADE )
  65.     { /* 2P/(Q-P) */
  66.     yd = yd - yn;
  67.     yn = 2.0 * yn;
  68.     }
  69. qyaprx = yn/yd;
  70. if( config & CW )
  71.     qyaprx = gx + qyaprx * gx * gx;
  72. if( config & SPECIAL )
  73.     qyaprx = speci( qyaprx, gx );
  74. return( qyaprx );
  75. }
  76.  
  77.  
  78.  
  79. /* Subroutine to compute approximation error at x */
  80.  
  81. double geterr(x)
  82. double x;
  83. {
  84. double e, f;
  85. double fabs(), approx(), func();
  86.  
  87. f = func(x);
  88. e = approx(x) - f;
  89. if( relerr )
  90.     {
  91.     if( f != 0.0 )
  92.         e /= f;
  93.     }
  94. if( e < 0.0 )
  95.     {
  96.     esign = -1;
  97.     e = -e;
  98.     }
  99. else
  100.     esign = 1;
  101.  
  102. return(e);
  103. }
  104.  
  105.  
  106.  
  107. /* Subroutine for special argument transformations */
  108.  
  109. double gofx(x)
  110. double x;
  111. {
  112.  
  113. return(x);
  114. }
  115.  
  116. /* Routine for special modifications of the approximating form.
  117.  * Example already provided by the CW flag: 
  118.  *    y(1+dev) = gx + gx^2 R(gx)
  119.  * would change y to
  120.  *    R(gx) = (y - gx)/(gx*gx)
  121.  * This function is called from remese.c.
  122.  *
  123.  * An inverse routine called speci() must also be supplied.
  124.  * This finds y from R and gx (see below).
  125.  */
  126. extern double PI, PIO4, THPIO4;
  127. #define LS2PI 0.91893853320467274178
  128. #define SQTPI 2.50662827463100050242
  129. #define JZ1   5.783185962946784521176
  130. #define JZ2  30.471262343662086399078
  131. #define JZ3  74.887006790695183444889
  132. #define YZ1  0.43221455686510834878
  133. #define YZ2 22.401876406482861405
  134. #define YZ3 64.130620282338755553298
  135. #define JO1 14.6819706421238932572
  136. #define YO1  4.66539330185668857532
  137.  
  138.  
  139. extern double TWOOPI;
  140.  
  141. double special( y, gx )
  142. double y, gx;
  143. {
  144. double a, b;
  145. /* exponential, y = exp(x) = 1 + x + .5x^2 + x^2 R(x) */
  146. /*
  147. if( gx == 0.0 )
  148.     return(1.0);
  149. b = gx * gx;
  150. a = (y - 1.0 - gx - .5*b)/b;
  151. */
  152.  
  153. /* y = exp2(x) = 1 + x R(x) */
  154. /*
  155. a = y - 1.0;
  156. */
  157.  
  158. /* y = cos(x) = 1 - .5 x^2 + x^2 x^2 P(x^2) */
  159. /*
  160. b = gx * gx;
  161. a = (y - 1.0 + 0.5*b)/b;
  162. */
  163.  
  164. /* logarithm, y = log(1+x) = x - .5x^2 + x^2 * (xP(x))
  165.  * configuration is  ZER | XPX | SPECIAL
  166.  */
  167. /*
  168. if( gx == 0.0 )
  169.     return(0.0);
  170. b = gx * gx;
  171. a = (y - gx + 0.5 * b)/b;
  172. */
  173.  
  174. /* y = log gamma(x) = q(x) + 1/x P(1/x^2)
  175.  * configufation is SQH | ZER | XPX | PXSQ | SPECIAL
  176.  */
  177. /*
  178. b = 1.0/gx;
  179. b = ( b - 0.5 ) * log(b) - b + LS2PI;
  180. a = y - b;
  181. */
  182.  
  183. /* y = gamma(x) = sqrt(2 pi) x^(x-.5) exp(-x) (1 + 1/x P(1/x)) */
  184. /*
  185. b = 1.0/gx;
  186. a = SQTPI * pow( b, b-0.5 ) * exp(-b);
  187. a = (y - a)/a;
  188. */
  189.  
  190. /* y = j0(x) = (x^2 - JZ1)(x^2-JZ2)(x^2-JZ3)P(x^2) */
  191. /*
  192. b=gx*gx;
  193. a = (b-JZ1);
  194. a = y/a;
  195. */
  196.  
  197. /* y = y0(x) = TWOOPI * log(x) * j0(x) + (x^2-YZ1)*P(x^2) */
  198. /*
  199. b=gx*gx;
  200. a = y - TWOOPI * log(gx) * j0(gx);
  201. a /= (b-YZ1);
  202. */
  203.  
  204. /* y = j1(x) = (x^2 - JO1)P(x^2) */
  205. /*
  206. b=gx*gx;
  207. a = (b-JO1);
  208. a = y/a;
  209. */
  210. /* y = y1(x) = TWOOPI * (log(x) * j1(x) - 1/x) + (x^2-YZ1)*P(x^2) */
  211. /*
  212. b = gx*gx;
  213. a = y - TWOOPI * ( j1(gx) * log(gx)  -  1.0/gx );
  214. a /= (b-YO1);
  215. */
  216.  
  217. /* y = exp2(x) = 1 + x P(x) */
  218. a = y - 1.0;
  219.  
  220. /* y = erfc(x) = exp(-x^2) P(x) */
  221. a = y * exp( gx * gx );
  222.  
  223. /* acosh() */
  224. /*
  225. if( gx == 0.0 )
  226.     return(0.0);
  227. a = y/(2.0*sqrt(gx));
  228. */
  229. return( a );
  230. }
  231.  
  232.  
  233.  
  234. double speci( R, gx )
  235. double R, gx;
  236. {
  237. double y, b;
  238.  
  239. /* exponential, y = exp(x) = 1 + x + .5x^2 + x^2 R(x) */
  240. /*
  241. b = gx * gx;
  242. y = 1.0 + gx + .5 * b;
  243. y = y + b * R;
  244. */
  245.  
  246. /* y = exp2(x) = 1 + x R(x) */
  247. /*
  248. y = R + 1.0;
  249. */
  250. /* y = cos(x) = 1 - .5 x^2 +  x^2 x^2 R(x^2) */
  251. /*
  252. b = gx * gx;
  253. y = 1.0 - 0.5*b + b * R;
  254. */
  255.  
  256. /* log(1+x) = x - .5x^2 + x^2 xR(x) */
  257. /*
  258. b = gx * gx;
  259. y = gx  - 0.5 * b  +  b * R;
  260. */
  261.  
  262. /* y = log gamma(x) = q(x) + 1/x P(1/x^2) */
  263. /*
  264. b = 1.0/gx;
  265. b = ( b - 0.5 ) * log(b) - b + LS2PI;
  266. y = b + R;
  267. */
  268.  
  269. /* y = gamma(x) = sqrt(2 pi) x^(x-.5) exp(-x) (1 + 1/x P(1/x)) */
  270. /*
  271. b = 1.0/gx;
  272. b = SQTPI * pow( b, b-0.5 ) * exp(-b);
  273. y = b + b * R;
  274. */
  275. /* y = j0(x) = (x^2 - JZ1)(x^2-JZ2)(x^2-JZ3)P(x^2) */
  276. /*
  277. b=gx*gx;
  278. y = (b-JZ1)*R;
  279. */
  280.  
  281. /* y = y0(x) = TWOOPI * log(x) * j0(x) + (x^2-YZ1)*P(x^2) */
  282. /*
  283. b=gx*gx;
  284. y = TWOOPI * log(gx) * j0(gx) + (b-YZ1)*R;
  285. */
  286.  
  287. /* y = j1(x) = (x^2 - JO1)P(x^2) */
  288. /*
  289. b=gx*gx;
  290. y = (b-JO1)*R;
  291. */
  292.  
  293. /* y = y1(x) = TWOOPI * (log(x) * j1(x) - 1/x) + (x^2-YZ1)*P(x^2) */
  294. /*
  295. b = gx*gx;
  296. y =  TWOOPI * ( j1(gx) * log(gx)  -  1.0/gx ) + (b-YO1)*R;
  297. */
  298.  
  299. /* y = exp2(x) = 1 + x P(x) */
  300. /*y = 1.0 + R;*/
  301.  
  302. /* y = erfc(x) = exp(-x^2) P(x) */
  303. y = exp( -gx * gx ) * R;
  304.  
  305. /*y = 2.0 * sqrt(gx) * R;*/
  306. return( y );
  307. }
  308.  
  309. /* Put here an accurate routine */
  310. /* for the function to be approximated. */
  311.  
  312. static int fflg = 0;
  313. static double ff = 0.0;
  314.  
  315. double func(x)
  316. double x;
  317. {
  318. double xx, y, t, u, s, c;
  319.  
  320.  
  321. qx = x;
  322.  
  323. if( fflg == 0 )
  324.     {
  325.     fflg = 1;
  326.     ff = 10.0 * log10(2.0);
  327.     }
  328.  
  329.  
  330. if( x == 0.0 )
  331.     {
  332. /*    y = ff;*/
  333.     y = 0.0;
  334.     goto done;
  335.     }
  336.  
  337. /* Bessel, phase */
  338. /*
  339. xx = 1.0/x;
  340. y = j0(xx);
  341. t = y0(xx);
  342. y = atan(t/y);
  343. t = xx - PIO4;
  344. t = t - PI * floor(t/PI + 0.5);
  345. y -= t;
  346. if( y > 0.5*PI )
  347.     y -= PI;
  348. if( y < -0.5*PI )
  349.     y += PI;
  350. */
  351. /* Bessel, modulus */
  352. /*
  353. xx = 1.0/(x*x);
  354. y = j1(xx);
  355. t = y1(xx);
  356. y = sqrt( y*y + t*t );
  357. */
  358.  
  359. /* bes1, phase */
  360. /*
  361. xx = 1.0/x;
  362. y = j1(xx);
  363. t = y1(xx);
  364. y = atan(t/y);
  365. t = xx - THPIO4;
  366. t = t - PI * floor(t/PI + 0.5);
  367. y -= t;
  368. if( y > 0.5*PI )
  369.     y -= PI;
  370. if( y < -0.5*PI )
  371.     y += PI;
  372. */
  373.  
  374. y = log(1.0+x);
  375.  
  376. done:
  377. qy = y;
  378. return( y );
  379. }
  380.