view release on metacpan or  search on metacpan

libmd/polmisc.c  view on Meta::CPAN

  int i, n;
#if 0
  double z[N+1];
  double u;
#endif

  if (nn > N)
    {
      mtherr ("polatn", OVERFLOW);
      return;
    }
  x = (double * )malloc( (MAXPOL+1) * sizeof (double) );
  y = (double * )malloc( (MAXPOL+1) * sizeof (double) );
  polmov( pol, nn, x );
  polclr( y, MAXPOL );

  /* Find lowest degree nonzero term.  */
  t = 0.0;
  for( n=0; n<nn; n++ )
    {
      if( x[n] != 0.0 )
	goto nzero;
    }
  polmov( y, nn, ans );
  return;

nzero:

  if( n > 0 )
    {
      if (n & 1)
        {
	  printf("error, sqrt of odd polynomial\n");
	  return;
	}
      /* Divide by x^n.  */
      y[n] = x[n];
      poldiv (y, nn, pol, N, x);
    }

  t = x[0];
  for( i=1; i<=nn; i++ )
    x[i] /= t;
  x[0] = 0.0;
  /* series development sqrt(1+x) = 1  +  x / 2  -  x**2 / 8  +  x**3 / 16
     hopes that first (constant) term is greater than what follows   */
  polsbt( x, nn, psqrt, nn, y);
  t = sqrt( t );
  for( i=0; i<=nn; i++ )
    y[i] *= t;

  /* If first nonzero coefficient was at degree n > 0, multiply by
     x^(n/2).  */
  if (n > 0)
    {
      polclr (x, MAXPOL);
      x[n/2] = 1.0;
      polmul (x, nn, y, nn, y);
    }
#if 0
/* Newton iterations */
for( n=0; n<10; n++ )
	{
	poldiv( y, nn, pol, nn, z );
	poladd( y, nn, z, nn, y );
	for( i=0; i<=nn; i++ )
		y[i] *= 0.5;
	for( i=0; i<=nn; i++ )
		{
		u = md_fabs( y[i] - z[i] );
		if( u > 1.0e-15 )
			goto more;
		}
	goto done;
more:	;
	}
printf( "square root did not converge\n" );
done:
#endif /* 0 */

polmov( y, nn, ans );
free( y );
free( x );
}



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

  if (nn > N)
    {
      mtherr ("polatn", OVERFLOW);
      return;
    }
  w = (double * )malloc( (MAXPOL+1) * sizeof (double) );
  c = (double * )malloc( (MAXPOL+1) * sizeof (double) );
  polmov( x, nn, w );
  polclr( c, MAXPOL );
  polclr( y, nn );
  /* a, in the description, is x[0].  b is the polynomial x - x[0].  */
  a = w[0];
  /* c = md_cos (b) */
  w[0] = 0.0;
  polsbt( w, nn, pcos, nn, c );
  sc = md_sin(a);
  /* md_sin(a) md_cos (b) */