Category : C Source Code
Archive   : CEPHES22.ZIP
Filename : REMESA.C

 
Output of file : REMESA.C contained in archive : CEPHES22.ZIP
/* remesa.c */

#include "remes.h"
#define PI 3.14159265358979323846

/* Subroutine to get operator instructions. */

remesa()
{
char s[40];
double a, b, c, q, r;
int i, k, ncheb;
double func(), cos(), sqrt();

/* Ask for error criterion */
printf( "Relative error (y or n) ? " );
gets( s ); /* Read in operator's response. */
relerr = 0;
if( s[0] == 'y' ) /* "y" means yes. */
relerr = 1;

/* Enable Automatic search. */
search = 1;

getconf:

printf( "Configuration word = " );
pconfig();
printf( " ? " );
gets(s);
if( s[0] != '\0' )
{
/* Display a menu of configuration bits */
if( s[0] == '?' )
{
k = config;
config = 0xffff;
pconfig();
config = k;
goto getconf;
}
else
{
sscanf( s, "%x", &config );
pconfig();
}
printf( "\n" );
}

printf( "Degree of numerator polynomial? " );
gets( s ); /* read line */
sscanf( s, "%d", &n ); /* decode characters */

printf( "Degree of denominator polynomial? " );
gets( s );
sscanf( s, "%d", &d );

printf ( "Start of approximation interval ? " );
gets( s );
sscanf( s, "%lf", &apstrt );

getwid:
printf ( "Width of approximation interval ? " );
gets( s );
sscanf( s, "%lf", &apwidt );
if( apwidt <= 0.0 )
{
puts( "? must be > 0" );
goto getwid;
}
apend = apstrt + apwidt;
nd1 = n + d + 1;
spread = 1.0e37;
iter = 0;
askitr = 1;

/* Supply initial guesses for solution points. */

if( (d == 0) && ((config & ZER) == 0) )
{ /* there is no denominator polynomial */
neq = n + 2; /* The number of equations to solve */
mm[neq] = apend;
}
else
{
neq = nd1;
}
ncheb = nd1; /* Degree of Chebyshev error estimate */

/* Find ncheb+1 extrema of Chebyshev polynomial */
a = ncheb;
mm[0] = apstrt;
for( i=1; i {
r = 0.5 * (1.0 - cos( (PI * i) / a ) );
if( config & SQL )
r = r * sqrt(r);
else if( config & SQH )
r = sqrt(r);
mm[i] = apstrt + r * apwidt;
}
mm[ncheb] = apend;

if( (d == 0) && ((config & ZER) == 0) )
{ /* The operative locations are maxima. */
for( i=0; i<=neq; i++ )
xx[i] = mm[i];
}
else
{ /* Zeros of Chebyshev polynomial */
a = 2.0 * ncheb;
for( i=0; i<=ncheb; i++ )
{
r = 0.5 * (1.0 - cos( PI * (2*i+1) / a ) );
/* Squeeze toward low end of approximation interval. */
if( config & SQL )
r = r * sqrt(r);
/* Squeeze toward high end. */
else if( config & SQH )
r = sqrt(r);
xx[i] = apstrt + r * apwidt;
}
/* Deviation at solution points */
dev = 0.0;
}

/* If form is xR(x), avoid x = 0 */
if( config & (XPX | X2PX) )
{
if( config & CW )
qx = 1.0e-15 + xx[0]/2.0;
else
qx = 1.0e-15 + apstrt * (1.0 + 1.0e-15);
mm[0] = qx;
printf( "mm[0] = %.15E\n", qx );
if( (d == 0) && ((config & ZER) == 0) )
xx[0] = qx;
}
/* Initialize step sizes */
stpini();
}



/* Subroutine to initialize step sizes. */
stpini()
{
int i;

if( search )
{
xx[neq+1] = apend;
delta = 0.25;
if( (d > 0) || ((config & ZER) != 0) )
{
step[0] = xx[0] - apstrt;
for( i=1; i step[i] = xx[i] - xx[i-1];
}
else
{
step[0] = 0.0625 * (xx[1] - xx[0]);
for( i=1; i step[i] = 0.0625 * (xx[i+1] - xx[i]);
}
step[neq] = step[neq-1];
}
}


/* Bits of configuration word */
#define NCNFG 11

static char *cnfmsg[NCNFG] = {
"PXSQ",
"XPX",
"PADE",
"CW",
"SQL",
"ZER",
"X2PX",
"SQH",
"SPECIAL",
"TRUNC",
"PXCU"
};

/* Subroutine to Display the configuration word */

pconfig()
{
int i, k;

k = 1;
for( i=0; i {
if( k & config )
printf( "%x:%s ", k, cnfmsg[i] );
k <<= 1;
}
printf( "hex value = %x ", config );
}


  3 Responses to “Category : C Source Code
Archive   : CEPHES22.ZIP
Filename : REMESA.C

  1. Very nice! Thank you for this wonderful archive. I wonder why I found it only now. Long live the BBS file archives!

  2. This is so awesome! 😀 I’d be cool if you could download an entire archive of this at once, though.

  3. But one thing that puzzles me is the “mtswslnkmcjklsdlsbdmMICROSOFT” string. There is an article about it here. It is definitely worth a read: http://www.os2museum.com/wp/mtswslnk/