home *** CD-ROM | disk | FTP | other *** search
/ Power-Programmierung / CD2.mdf / c / library / dos / math / cephes / remes / remess.c < prev    next >
Encoding:
C/C++ Source or Header  |  1992-11-17  |  2.7 KB  |  171 lines
  1. /* remess.c */
  2. #include "remes.h"
  3.  
  4. /* Search subroutine */
  5.  
  6. remess()
  7. {
  8. double a, b, q, xm, ym, xn, yn, xx0, xx1;
  9. int i, j, meq, emsign, ensign, steps;
  10. double approx(), func(), geterr();
  11.  
  12. /* Search for maxima of error */
  13. eclose = 1e30;
  14. farther = 0.0;
  15. j = 0;
  16. meq = neq;
  17.  
  18. if( (d > 0) || ((config & ZER) != 0) )
  19.     {
  20.     j = 1;
  21.     meq += 1;
  22.     }
  23. xx0 = apstrt;
  24.  
  25. for( i=0; i<meq; i++ )
  26.     {
  27.     steps = 0;
  28.     if( (d > 0) || ((config & ZER) != 0) )
  29.         xx1 = xx[i]; /* Next zero */
  30.     else
  31.         xx1 = mm[i+1]; /* Next maximum */
  32.     if( i == meq-1 )
  33.         xx1 = apend;
  34.     xm = mm[i];
  35.     ym = geterr(xm);
  36.     emsign = esign; /* Sign of error */
  37.     q = step[i];
  38.     xn = xm + q;
  39. /* Cannot skip over adjacent boundaries */
  40.     if( xn < xx0 )
  41.         goto revers;
  42.     if( xn >= xx1 )
  43.         goto revers;
  44.     yn = geterr(xn);
  45.     ensign = esign;
  46.     if( yn < ym )
  47.         {
  48. revers:
  49.         q = -q;
  50.         xn = xm;
  51.         yn = ym;
  52.         ensign = emsign;
  53.         }
  54.  
  55. /* March until error becomes smaller. */
  56.  
  57.     while( yn >= ym )
  58.         {
  59.         if( ++steps > 10 )
  60.             goto tsear;
  61.         ym = yn;
  62.         xm = xn;
  63.         emsign = ensign;
  64.         a = xm + q;
  65.         if( a == xm )
  66.             goto tsear;
  67. /* Must not skip over the zeros either side. */
  68.         if( a <= xx0 )
  69.             goto tsear;
  70.         if( a >= xx1 )
  71.             goto tsear;
  72.         xn = a;
  73.         yn = geterr(xn);
  74.         ensign = esign;
  75.         }
  76.  
  77. tsear:
  78.     mm[i] = xm; /* Position of maximum */
  79.     yy[i] = ym; /* Value of maximum */
  80.  
  81.     if( ym == 0.0 )
  82.         printf( "yy[%d] = 0\n", i );
  83.     if( eclose > ym )
  84.         eclose = ym;
  85.  
  86.     if( farther < ym )
  87.         farther = ym;
  88.  
  89.  /* No denominator */
  90.     if( (d == 0) && ((config & ZER) == 0) )
  91.         xx[i] = xm;
  92. /* Walk to next zero. */
  93.     if( (d > 0) || ((config & ZER) != 0) )
  94.         xx0 = xx1;
  95.     else
  96.         xx0 = 0.5*(xm + xx1);
  97.     } /* end of search loop */
  98.  
  99.  
  100. /* Decrease step size if error spread increased. */
  101. q = (farther - eclose);
  102. /* Relative error spread */
  103. if( eclose != 0.0 )
  104.     q /= eclose;
  105. if( q >= spread )
  106.     { /* Spread is increasing; decrease step size. */
  107.     if( config & TRUNC )
  108.         delta *= 0.875;
  109.     else
  110.         delta *= 0.5;
  111.     printf( "delta = %.4E\n", delta );
  112.     }
  113.  spread = q;
  114.  
  115. printf( "peak error = %.4E, relative error spread = %.4E\n",
  116.         farther, spread );
  117.  
  118. /* Calculate new step sizes */
  119.  
  120. if( (d == 0) && ((config & ZER) == 0) )
  121.     {
  122.     if( spread < 0.25 )
  123.         q = 0.0625;
  124.     else
  125.         q = 0.5;
  126.  
  127.     for( i=0; i<neq; i++ )
  128.         step[i] *= q;
  129.     }
  130. else
  131.     {
  132.  
  133.     for( i=0; i<neq; i++ )
  134.         {
  135.         q = yy[i+1];
  136.         if( q != 0.0 )
  137.             q = yy[i]/q  - 1.0;
  138.         else
  139.             q = 0.0625;
  140.         if( q > 0.25 )
  141.             q = 0.25;
  142.         q *= mm[i+1] - mm[i];
  143.         step[i] = q * delta;
  144.         }
  145.     step[neq] = step[neq-1];
  146.  
  147.  
  148. /* Insert new locations for the zeros. */
  149.     for( i=0; i<neq; i++ )
  150.         {
  151.         xm = xx[i] - step[i];
  152.         if( xm <= apstrt )
  153.             continue;
  154.         if( xm >= apend )
  155.             continue;
  156.         if( xm <= mm[i] )
  157.             {
  158.             printf( "xx[%d] < mm\n", i );
  159.             xm = 0.5 * (mm[i] + xx[i]);
  160.             }
  161.         if( xm >= mm[i+1] )
  162.             {
  163.             printf( "xx[%d] > mm\n", i );
  164.             xm = 0.5 * (mm[i+1] + xx[i]);
  165.             }
  166.         xx[i] = xm;
  167.         }
  168.     }
  169. sdone:    ;
  170. }
  171.