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

 
Output of file : POLMISC.C contained in archive : CEPHES22.ZIP
/* square root, sine, and arctangent of polynomial
* See polyn.c for data structures and discussion.
*/


/* Highest degree of polynomial to be handled
* by the polyn.c subroutine package */
#define N 16

/* Taylor series coefficients for various functions
*/
double patan[N+1] = {
0.0, 1.0, 0.0, -1.0/3.0, 0.0,
1.0/5.0, 0.0, -1.0/7.0, 0.0, 1.0/9.0, 0.0, -1.0/11.0,
0.0, 1.0/13.0, 0.0, -1.0/15.0, 0.0 };

double psin[N+1] = {
0.0, 1.0, 0.0, -1.0/6.0, 0.0, 1.0/120.0, 0.0,
-1.0/5040.0, 0.0, 1.0/362880.0, 0.0, -1.0/39916800.0,
0.0, 1.0/6227020800.0, 0.0, -1.0/1.307674368e12, 0.0};

double pcos[N+1] = {
1.0, 0.0, -1.0/2.0, 0.0, 1.0/24.0, 0.0,
-1.0/720.0, 0.0, 1.0/40320.0, 0.0, -1.0/3628800.0, 0.0,
1.0/479001600.0, 0.0, -1.0/8.7179291e10, 0.0, 1.0/2.0922789888e13};

double pasin[N+1] = {
0.0, 1.0, 0.0, 1.0/6.0, 0.0,
3.0/40.0, 0.0, 15.0/336.0, 0.0, 105.0/3456.0, 0.0, 945.0/42240.0,
0.0, 10395.0/599040.0 , 0.0, 135135.0/9676800.0 , 0.0
};


static double polt[N+1];
static double polq[N+1];
static double polu[N+1];



/* Arctangent of the ratio num/den of two polynomials.
*/
polatn( den, num, ans )
double num[], den[], ans[];
{
double a, t;
int i;
double atan2();

/*
* arctan( a + b ) = arctan(a) + arctan( b/(1 + ab + a**2) )
*/
t = num[0];
a = den[0];
if( (t == 0.0) && (a == 0.0 ) )
{
t = num[1];
a = den[1];
}
t = atan2( a, t ); /* arctan(a) */
polclr( polq, N );
i = poldiv( den, N, num, N, polq );
a = polq[0]; /* a */
polq[0] = 0.0; /* b */
polmov( polq, N, polu ); /* b */
/*
* Form the polynomial
* 1 + ab + a**2
* where a is a scalar.
*/
for( i=0; i<=N; i++ )
polu[i] *= a;
polu[0] += 1.0 + a * a;
poldiv( polu, N, polq, N, polt ); /* divide into b */
polsbt( polt, N, patan, N, polu ); /* arctan(b) */
polu[0] += t; /* plus arctan(a) */
}



/* Square root of a polynomial.
* Assumes the lowest degree nonzero term is dominant
* and of even degree. An error message is given
* if the Newton iteration does not converge.
*/
polsqt( pol, ans )
double pol[], ans[];
{
double x[N+1], y[N+1], z[N+1];
double t, u;
int i, j, n0;
double sqrt(), fabs();

polmov( pol, N, x );
polclr( y, N );
/* Initial guess is square root of lowest order
* nonzero term
*/
t = 0.0;
for( i=0; i {
if( x[i] != 0.0 )
goto nzero;
}
polmov( y, N, ans );
return;

nzero:

if( i > 0 )
{
t = x[i];
i /= 2;
y[i] = sqrt( t );
}
else
{
t = x[0];
n0 = 0;
y[0] = 1.0;
for( i=1; i<=N; i++ )
x[i] /= t;
x[0] = 0.0;
/* series development sqrt(1+x) = 1 + x / 2 - x**2 / 8 + x**3 / 16
* assumes first (constant) term is greater than what follows
*/
polmov( x, N, z );
for( i=0; i<=N; i++ )
z[i] *= 0.5;
poladd( z, N, y, N, y );
polmul( x, N, x, N, z );
for( i=0; i<=N; i++ )
z[i] *= 0.25;
polsub( z, N, y, N, y );
polmul( x, N, z, N, z );
for( i=0; i<=N; i++ )
z[i] *= 0.5;
poladd( z, N, y, N, y );
t = sqrt( t );
for( i=0; i<=N; i++ )
y[i] *= t;
}
/* Newton iterations */
for( j=0; j<10; j++ )
{
poldiv( y, N, pol, N, z );
poladd( y, N, z, N, y );
for( i=0; i<=N; i++ )
y[i] *= 0.5;
for( i=0; i<=N; i++ )
{
u = fabs( y[i] - z[i] );
if( u > 1.0e-15 )
goto more;
}
goto done;
more: ;
}
printf( "square root did not converge\n" );
done:
polmov( y, N, ans );
}



/* Sine of a polynomial.
* The computation uses
* sin(a+b) = sin(a) cos(b) + cos(a) sin(b)
* where a is the constant term of the polynomial and
* b is the sum of the rest of the terms.
* Since sin(b) and cos(b) are computed by series expansions,
* the value of b should be small.
*/
polsin( x, y )
double x[], y[];
{
double w[N+1], c[N+1];
double a, sc;
int i;
double sin(), cos();

polmov( x, N, w );
polclr( c, N );
polclr( y, N );
a = w[0];
w[0] = 0.0;
polsbt( w, N, pcos, N, c );
sc = sin(a);
for( i=0; i<=N; i++ )
c[i] *= sc;
polsbt( w, N, psin, N, y );
sc = cos(a);
for( i=0; i<=N; i++ )
y[i] *= sc;
poladd( c, N, y, N, y );
}


  3 Responses to “Category : C Source Code
Archive   : CEPHES22.ZIP
Filename : POLMISC.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/