home *** CD-ROM | disk | FTP | other *** search
/ Power-Programmierung / CD2.mdf / c / library / dos / math / cephes / remes / remesa.c < prev    next >
Encoding:
C/C++ Source or Header  |  1992-11-17  |  3.4 KB  |  203 lines
  1. /* remesa.c */
  2.  
  3. #include "remes.h"
  4. #define PI 3.14159265358979323846
  5.  
  6. /* Subroutine to get operator instructions. */
  7.  
  8. remesa()
  9. {
  10. char s[40];
  11. double a, b, c, q, r;
  12. int i, k, ncheb;
  13. double func(), cos(), sqrt();
  14.  
  15. /* Ask for error criterion */
  16. printf( "Relative error (y or n) ? " );
  17. gets( s ); /* Read in operator's response. */
  18. relerr = 0;
  19. if( s[0] == 'y' ) /* "y" means yes. */
  20.     relerr = 1;
  21.  
  22. /* Enable Automatic search. */
  23. search = 1;
  24.  
  25. getconf:
  26.  
  27. printf( "Configuration word = " );
  28. pconfig();
  29. printf( " ? " );
  30. gets(s);
  31. if( s[0] != '\0' )
  32.     {
  33. /* Display a menu of configuration bits */
  34.     if( s[0] == '?' )
  35.         {
  36.         k = config;
  37.         config = 0xffff;
  38.         pconfig();
  39.         config = k;
  40.         goto getconf;
  41.         }
  42.     else
  43.         {
  44.         sscanf( s, "%x", &config );
  45.         pconfig();
  46.         }
  47.     printf( "\n" );
  48.     }
  49.  
  50. printf( "Degree of numerator polynomial? " );
  51. gets( s );        /* read line */
  52. sscanf( s, "%d", &n );    /* decode characters */
  53.  
  54. printf( "Degree of denominator polynomial? " );
  55. gets( s );
  56. sscanf( s, "%d", &d );
  57.  
  58. printf ( "Start of approximation interval ? " );
  59. gets( s );
  60. sscanf( s, "%lf", &apstrt );
  61.  
  62. getwid:
  63. printf ( "Width of approximation interval ? " );
  64. gets( s );
  65. sscanf( s, "%lf", &apwidt );
  66. if( apwidt <= 0.0 )
  67.     {
  68.     puts( "? must be > 0" );
  69.     goto getwid;
  70.     }
  71. apend = apstrt + apwidt;
  72. nd1 = n + d + 1;
  73. spread = 1.0e37;
  74. iter = 0;
  75. askitr = 1;
  76.  
  77. /* Supply initial guesses for solution points. */
  78.  
  79. if( (d == 0) && ((config & ZER) == 0) )
  80.     { /* there is no denominator polynomial */
  81.     neq = n + 2; /* The number of equations to solve */
  82.     mm[neq] = apend;
  83.     }
  84. else
  85.     {
  86.     neq = nd1;
  87.     }
  88. ncheb = nd1; /* Degree of Chebyshev error estimate */
  89.  
  90. /* Find ncheb+1 extrema of Chebyshev polynomial */
  91. a = ncheb;
  92. mm[0] = apstrt;
  93. for( i=1; i<ncheb; i++ )
  94.     {
  95.     r = 0.5 * (1.0 - cos( (PI * i) / a ) );
  96.     if( config & SQL )
  97.         r = r * sqrt(r);
  98.     else if( config & SQH )
  99.         r = sqrt(r);
  100.     mm[i] = apstrt + r * apwidt;
  101.     }
  102. mm[ncheb] = apend;
  103.  
  104. if( (d == 0) && ((config & ZER) == 0) )
  105.     { /* The operative locations are maxima. */
  106.     for( i=0; i<=neq; i++ )
  107.         xx[i] = mm[i];
  108.     }
  109. else
  110.     { /* Zeros of Chebyshev polynomial */
  111.     a = 2.0 * ncheb;
  112.     for( i=0; i<=ncheb; i++ )
  113.         {
  114.         r = 0.5 * (1.0 - cos( PI * (2*i+1) / a ) );
  115. /* Squeeze toward low end of approximation interval. */
  116.         if( config & SQL )
  117.             r = r * sqrt(r);
  118. /* Squeeze toward high end. */
  119.         else if( config & SQH )
  120.             r = sqrt(r);
  121.         xx[i] = apstrt + r * apwidt;
  122.         }
  123. /* Deviation at solution points */
  124.     dev = 0.0;
  125.     }
  126.  
  127. /* If form is xR(x), avoid x = 0 */
  128. if( config & (XPX | X2PX) )
  129.     {
  130.     if( config & CW )
  131.         qx = 1.0e-15 + xx[0]/2.0;
  132.     else
  133.         qx =  1.0e-15 + apstrt * (1.0 + 1.0e-15);
  134.     mm[0] = qx;
  135.     printf( "mm[0] = %.15E\n", qx );
  136.     if( (d == 0) && ((config & ZER) == 0) )
  137.         xx[0] = qx;
  138.     }
  139. /* Initialize step sizes */
  140. stpini();
  141. }
  142.  
  143.  
  144.  
  145. /* Subroutine to initialize step sizes. */
  146. stpini()
  147. {
  148. int i;
  149.  
  150. if( search )
  151.     {
  152.     xx[neq+1] = apend;
  153.     delta = 0.25;
  154.     if( (d > 0) || ((config & ZER) != 0) )
  155.         {
  156.         step[0] = xx[0] - apstrt;
  157.         for( i=1; i<neq; i++ )
  158.             step[i] = xx[i] - xx[i-1];
  159.         }
  160.     else
  161.         {
  162.         step[0] = 0.0625 * (xx[1] - xx[0]);
  163.         for( i=1; i<neq; i++ )
  164.             step[i] = 0.0625 * (xx[i+1] - xx[i]);
  165.         }
  166.     step[neq] = step[neq-1];
  167.     }
  168. }
  169.  
  170.  
  171. /* Bits of configuration word */
  172. #define NCNFG 11
  173.  
  174. static char *cnfmsg[NCNFG] = {
  175. "PXSQ",
  176. "XPX",
  177. "PADE",
  178. "CW",
  179. "SQL",
  180. "ZER",
  181. "X2PX",
  182. "SQH",
  183. "SPECIAL",
  184. "TRUNC",
  185. "PXCU"
  186. };
  187.  
  188. /* Subroutine to Display the configuration word */
  189.  
  190. pconfig()
  191. {
  192. int i, k;
  193.  
  194. k = 1;
  195. for( i=0; i<NCNFG; i++ )
  196.     {
  197.     if( k & config )
  198.         printf( "%x:%s ", k, cnfmsg[i] );
  199.     k <<= 1;
  200.     }
  201. printf( "hex value = %x ", config );
  202. }
  203.