lib/libm/polmisc.c

   1
   2 /* Square root, sine, cosine, and arctangent of polynomial.
   3  * See polyn.c for data structures and discussion.
   4  */
   5
   6 #include <stdio.h>
   7 #include "mconf.h"
   8 #ifdef ANSIPROT
   9 extern double atan2 ( double, double );
  10 extern double sqrt ( double );
  11 extern double fabs ( double );
  12 extern double sin ( double );
  13 extern double cos ( double );
  14 extern void polclr ( double *a, int n );
  15 extern void polmov ( double *a, int na, double *b );
  16 extern void polmul ( double a[], int na, double b[], int nb, double c[] );
  17 extern void poladd ( double a[], int na, double b[], int nb, double c[] );
  18 extern void polsub ( double a[], int na, double b[], int nb, double c[] );
  19 extern int poldiv ( double a[], int na, double b[], int nb, double c[] );
  20 extern void polsbt ( double a[], int na, double b[], int nb, double c[] );
  21 extern void * malloc ( long );
  22 extern void free ( void * );
  23 #else
  24 double atan2(), sqrt(), fabs(), sin(), cos();
  25 void polclr(), polmov(), polsbt(), poladd(), polsub(), polmul();
  26 int poldiv();
  27 void * malloc();
  28 void free ();
  29 #endif
  30
  31 /* Highest degree of polynomial to be handled
  32    by the polyn.c subroutine package.  */
  33 #define N 16
  34 /* Highest degree actually initialized at runtime.  */
  35 extern int MAXPOL;
  36
  37 /* Taylor series coefficients for various functions
  38  */
  39 double patan[N+1] = {
  40   0.0,     1.0,      0.0, -1.0/3.0,     0.0,
  41   1.0/5.0, 0.0, -1.0/7.0,      0.0, 1.0/9.0, 0.0, -1.0/11.0,
  42   0.0, 1.0/13.0, 0.0, -1.0/15.0, 0.0 };
  43
  44 double psin[N+1] = {
  45   0.0, 1.0, 0.0,   -1.0/6.0,  0.0, 1.0/120.0,  0.0,
  46   -1.0/5040.0, 0.0, 1.0/362880.0, 0.0, -1.0/39916800.0,
  47   0.0, 1.0/6227020800.0, 0.0, -1.0/1.307674368e12, 0.0};
  48
  49 double pcos[N+1] = {
  50   1.0, 0.0,   -1.0/2.0,  0.0, 1.0/24.0,  0.0,
  51   -1.0/720.0, 0.0, 1.0/40320.0, 0.0, -1.0/3628800.0, 0.0,
  52   1.0/479001600.0, 0.0, -1.0/8.7179291e10, 0.0, 1.0/2.0922789888e13};
  53
  54 double pasin[N+1] = {
  55   0.0,     1.0,  0.0, 1.0/6.0,  0.0,
  56   3.0/40.0, 0.0, 15.0/336.0, 0.0, 105.0/3456.0, 0.0, 945.0/42240.0,
  57   0.0, 10395.0/599040.0 , 0.0, 135135.0/9676800.0 , 0.0
  58 };
  59
  60 /* Square root of 1 + x.  */
  61 double psqrt[N+1] = {
  62   1.0, 1./2., -1./8., 1./16., -5./128., 7./256., -21./1024., 33./2048.,
  63   -429./32768., 715./65536., -2431./262144., 4199./524288., -29393./4194304.,
  64   52003./8388608., -185725./33554432., 334305./67108864.,
  65   -9694845./2147483648.};
  66
  67 /* Arctangent of the ratio num/den of two polynomials.
  68  */
  69 void
  70 polatn( num, den, ans, nn )
  71      double num[], den[], ans[];
  72      int nn;
  73 {
  74   double a, t;
  75   double *polq, *polu, *polt;
  76   int i;
  77
  78   if (nn > N)
  79     {
  80       mtherr ("polatn", OVERFLOW);
  81       return;
  82     }
  83   /* arctan( a + b ) = arctan(a) + arctan( b/(1 + ab + a**2) ) */
  84   t = num[0];
  85   a = den[0];
  86   if( (t == 0.0) && (a == 0.0 ) )
  87     {
  88       t = num[1];
  89       a = den[1];
  90     }
  91   t = atan2( t, a );  /* arctan(num/den), the ANSI argument order */
  92   polq = (double * )malloc( (MAXPOL+1) * sizeof (double) );
  93   polu = (double * )malloc( (MAXPOL+1) * sizeof (double) );
  94   polt = (double * )malloc( (MAXPOL+1) * sizeof (double) );
  95   polclr( polq, MAXPOL );
  96   i = poldiv( den, nn, num, nn, polq );
  97   a = polq[0]; /* a */
  98   polq[0] = 0.0; /* b */
  99   polmov( polq, nn, polu ); /* b */
 100   /* Form the polynomial
 101      1 + ab + a**2
 102      where a is a scalar.  */
 103   for( i=0; i<=nn; i++ )
 104     polu[i] *= a;
 105   polu[0] += 1.0 + a * a;
 106   poldiv( polu, nn, polq, nn, polt ); /* divide into b */
 107   polsbt( polt, nn, patan, nn, polu ); /* arctan(b)  */
 108   polu[0] += t; /* plus arctan(a) */
 109   polmov( polu, nn, ans );
 110   free( polt );
 111   free( polu );
 112   free( polq );
 113 }
 114
 115
 116
 117 /* Square root of a polynomial.
 118  * Assumes the lowest degree nonzero term is dominant
 119  * and of even degree.  An error message is given
 120  * if the Newton iteration does not converge.
 121  */
 122 void
 123 polsqt( pol, ans, nn )
 124      double pol[], ans[];
 125      int nn;
 126 {
 127   double t;
 128   double *x, *y;
 129   int i, n;
 130 #if 0
 131   double z[N+1];
 132   double u;
 133 #endif
 134
 135   if (nn > N)
 136     {
 137       mtherr ("polatn", OVERFLOW);
 138       return;
 139     }
 140   x = (double * )malloc( (MAXPOL+1) * sizeof (double) );
 141   y = (double * )malloc( (MAXPOL+1) * sizeof (double) );
 142   polmov( pol, nn, x );
 143   polclr( y, MAXPOL );
 144
 145   /* Find lowest degree nonzero term.  */
 146   t = 0.0;
 147   for( n=0; n<nn; n++ )
 148     {
 149       if( x[n] != 0.0 )
 150         goto nzero;
 151     }
 152   polmov( y, nn, ans );
 153   return;
 154
 155 nzero:
 156
 157   if( n > 0 )
 158     {
 159       if (n & 1)
 160         {
 161           printf("error, sqrt of odd polynomial\n");
 162           return;
 163         }
 164       /* Divide by x^n.  */
 165       y[n] = x[n];
 166       poldiv (y, nn, pol, N, x);
 167     }
 168
 169   t = x[0];
 170   for( i=1; i<=nn; i++ )
 171     x[i] /= t;
 172   x[0] = 0.0;
 173   /* series development sqrt(1+x) = 1  +  x / 2  -  x**2 / 8  +  x**3 / 16
 174      hopes that first (constant) term is greater than what follows   */
 175   polsbt( x, nn, psqrt, nn, y);
 176   t = sqrt( t );
 177   for( i=0; i<=nn; i++ )
 178     y[i] *= t;
 179
 180   /* If first nonzero coefficient was at degree n > 0, multiply by
 181      x^(n/2).  */
 182   if (n > 0)
 183     {
 184       polclr (x, MAXPOL);
 185       x[n/2] = 1.0;
 186       polmul (x, nn, y, nn, y);
 187     }
 188 #if 0
 189 /* Newton iterations */
 190 for( n=0; n<10; n++ )
 191         {
 192         poldiv( y, nn, pol, nn, z );
 193         poladd( y, nn, z, nn, y );
 194         for( i=0; i<=nn; i++ )
 195                 y[i] *= 0.5;
 196         for( i=0; i<=nn; i++ )
 197                 {
 198                 u = fabs( y[i] - z[i] );
 199                 if( u > 1.0e-15 )
 200                         goto more;
 201                 }
 202         goto done;
 203 more:   ;
 204         }
 205 printf( "square root did not converge\n" );
 206 done:
 207 #endif /* 0 */
 208
 209 polmov( y, nn, ans );
 210 free( y );
 211 free( x );
 212 }
 213
 214
 215
 216 /* Sine of a polynomial.
 217  * The computation uses
 218  *     sin(a+b) = sin(a) cos(b) + cos(a) sin(b)
 219  * where a is the constant term of the polynomial and
 220  * b is the sum of the rest of the terms.
 221  * Since sin(b) and cos(b) are computed by series expansions,
 222  * the value of b should be small.
 223  */
 224 void
 225 polsin( x, y, nn )
 226      double x[], y[];
 227      int nn;
 228 {
 229   double a, sc;
 230   double *w, *c;
 231   int i;
 232
 233   if (nn > N)
 234     {
 235       mtherr ("polatn", OVERFLOW);
 236       return;
 237     }
 238   w = (double * )malloc( (MAXPOL+1) * sizeof (double) );
 239   c = (double * )malloc( (MAXPOL+1) * sizeof (double) );
 240   polmov( x, nn, w );
 241   polclr( c, MAXPOL );
 242   polclr( y, nn );
 243   /* a, in the description, is x[0].  b is the polynomial x - x[0].  */
 244   a = w[0];
 245   /* c = cos (b) */
 246   w[0] = 0.0;
 247   polsbt( w, nn, pcos, nn, c );
 248   sc = sin(a);
 249   /* sin(a) cos (b) */
 250   for( i=0; i<=nn; i++ )
 251     c[i] *= sc;
 252   /* y = sin (b)  */
 253   polsbt( w, nn, psin, nn, y );
 254   sc = cos(a);
 255   /* cos(a) sin(b) */
 256   for( i=0; i<=nn; i++ )
 257     y[i] *= sc;
 258   poladd( c, nn, y, nn, y );
 259   free( c );
 260   free( w );
 261 }
 262
 263
 264 /* Cosine of a polynomial.
 265  * The computation uses
 266  *     cos(a+b) = cos(a) cos(b) - sin(a) sin(b)
 267  * where a is the constant term of the polynomial and
 268  * b is the sum of the rest of the terms.
 269  * Since sin(b) and cos(b) are computed by series expansions,
 270  * the value of b should be small.
 271  */
 272 void
 273 polcos( x, y, nn )
 274      double x[], y[];
 275      int nn;
 276 {
 277   double a, sc;
 278   double *w, *c;
 279   int i;
 280   double sin(), cos();
 281
 282   if (nn > N)
 283     {
 284       mtherr ("polatn", OVERFLOW);
 285       return;
 286     }
 287   w = (double * )malloc( (MAXPOL+1) * sizeof (double) );
 288   c = (double * )malloc( (MAXPOL+1) * sizeof (double) );
 289   polmov( x, nn, w );
 290   polclr( c, MAXPOL );
 291   polclr( y, nn );
 292   a = w[0];
 293   w[0] = 0.0;
 294   /* c = cos(b)  */
 295   polsbt( w, nn, pcos, nn, c );
 296   sc = cos(a);
 297   /* cos(a) cos(b)  */
 298   for( i=0; i<=nn; i++ )
 299     c[i] *= sc;
 300   /* y = sin(b) */
 301   polsbt( w, nn, psin, nn, y );
 302   sc = sin(a);
 303   /* sin(a) sin(b) */
 304   for( i=0; i<=nn; i++ )
 305     y[i] *= sc;
 306   polsub( y, nn, c, nn, y );
 307   free( c );
 308   free( w );
 309 }