home *** CD-ROM | disk | FTP | other *** search
/ Power-Programmierung / CD2.mdf / c / library / dos / math / cephes / remes / remese.c < prev    next >
Encoding:
C/C++ Source or Header  |  1992-11-17  |  2.1 KB  |  120 lines
  1. /*        remese.c    */
  2.  
  3. /* Solve equations */
  4.  
  5. #include "remes.h"
  6.  
  7. remese()
  8. {
  9. double x, y, z, gx, gxsq, pade;
  10. float sprec;
  11. int i, j, ip;
  12. double gofx(), func(), special();
  13.  
  14. if( (d == 0) && ((config & ZER) == 0) )
  15.     { /* no denominator: */
  16. /* In this case the deviation is an unknown variable */
  17. /* adjoined to the other unknowns; it has a coefficient */
  18. /* of plus or minus 1. */
  19.     dev = 1.0;
  20.     }
  21. else
  22.     {
  23. /* Otherwise dev = 0, since the solution sought is */
  24. /* for the locations of zero error. */
  25.     dev = 0.0;
  26.     }
  27.  
  28. /* Set up the equations for solution by simq() */
  29. for( i=0; i<neq; i++ )
  30.     {
  31. /* Offset to 1st element of this row of matrix */
  32.     ip = neq * i;
  33. /* The guess for this row */
  34.     x = xx[i];
  35. /* Right-hand-side vector */
  36.     y = func(x);
  37.     gx = gofx(x);
  38.     if( config & PXCU )
  39.         gxsq = gx * gx * gx;
  40.     if( config & PXSQ )
  41.         gxsq = gx * gx;
  42.     if( config & PADE )
  43.         {
  44.         pade = 2.0 + y;
  45.         }
  46.     if( d > 0 )
  47.         { /* add the deviation if rational form */
  48. /* Relative error criterion */
  49.         if( relerr )
  50.             y = y * (1.0+dev);
  51. /* Absolute error criterion */
  52.         else
  53.             y = y + dev;
  54.         }
  55.     if( config & CW )
  56.         { /* y(1+dev) = z + z^2 P/Q */
  57.         y = (y - gx)/(gx*gx);
  58.         }
  59.     if( config & SPECIAL )
  60.         {
  61.         y = special( y, gx );
  62.         }
  63. /* Insert powers of x[i] for numerator coefficients. */
  64.     if( config & XPX )
  65.         z = gx;
  66.     else if( config & X2PX )
  67.         z = gx * gx;
  68.     else
  69.         z = 1.0;
  70.     for( j=0; j<=n; j++ )
  71.         {
  72.         if( config & PADE )
  73.             AA[ip+j] = pade * z;
  74.         else
  75.             AA[ip+j] = z;
  76.         if( config & (PXSQ | PXCU) )
  77.             z = z * gxsq;
  78.         else
  79.             z = z * gx;
  80.         }
  81. /* Insert denominator terms, if any. */
  82.     if( d > 0 )
  83.         {
  84.         z = 1.0;
  85.         for( j=0; j<d; j++ )
  86.             {
  87.             AA[ip+n+1+j] = -y * z;
  88.             if( config & (PXSQ | PXCU) )
  89.                 z = z * gxsq;
  90.             else
  91.                 z = z * gx;
  92.             }
  93. /* Right hand side vector */
  94.         BB[i] = y * z;
  95.         }
  96.     else
  97.         { /* No denominator */
  98. /* Right hand side vector */
  99.         BB[i] = y;
  100.         z = dev;
  101.         if( relerr )
  102.             z = z * y;
  103.         AA[ip+n+1] = z;
  104.         }
  105. /* Switch sign of deviation for next row. */
  106.     dev = -1.0 * dev;
  107.     }
  108.  
  109. /* Solve the simultaneous linear equations. */
  110. simq( &AA[0], &BB[0], ¶m[0], neq, 0, &IPS[0] );
  111. if( config & TRUNC )
  112.     {
  113.     for( i=0; i<neq; i++ )
  114.         {
  115.         sprec = (float )param[i];
  116.         param[i] = sprec;
  117.         }
  118.     }
  119. }
  120.