lang/cem/libcc/math/tail_m.a

   1 e�asin.c\0\0\0\0\0\0\0\0\0\0\0\0\ 2\ 2�\ 1\0\0q\ 5/*
   2  * (c) copyright 1988 by the Vrije Universiteit, Amsterdam, The Netherlands.
   3  * See the copyright notice in the ACK home directory, in the file "Copyright".
   4  *
   5  * Author: Ceriel J.H. Jacobs
   6  */
   7
   8 /* $Id: asin.c,v 1.5 1994/06/24 12:22:30 ceriel Exp $ */
   9
  10 #include <math.h>
  11 #include <errno.h>
  12
  13 extern int errno;
  14
  15 static double
  16 asin_acos(x, cosfl)
  17         double x;
  18 {
  19         int     negative = x < 0;
  20         int     i;
  21         double  g;
  22         extern double   sqrt();
  23         static double p[] = {
  24                 -0.27368494524164255994e+2,
  25                  0.57208227877891731407e+2,
  26                 -0.39688862997540877339e+2,
  27                  0.10152522233806463645e+2,
  28                 -0.69674573447350646411e+0
  29         };
  30         static double q[] = {
  31                 -0.16421096714498560795e+3,
  32                  0.41714430248260412556e+3,
  33                 -0.38186303361750149284e+3,
  34                  0.15095270841030604719e+3,
  35                 -0.23823859153670238830e+2,
  36                  1.0
  37         };
  38
  39         if (negative) {
  40                 x = -x;
  41         }
  42         if (x > 0.5) {
  43                 i = 1;
  44                 if (x > 1) {
  45                         errno = EDOM;
  46                         return 0;
  47                 }
  48                 g = 0.5 - 0.5 * x;
  49                 x = - sqrt(g);
  50                 x += x;
  51         }
  52         else {
  53                 /* ??? avoid underflow ??? */
  54                 i = 0;
  55                 g = x * x;
  56         }
  57         x += x * g * POLYNOM4(g, p) / POLYNOM5(g, q);
  58         if (cosfl) {
  59                 if (! negative) x = -x;
  60         }
  61         if ((cosfl == 0) == (i == 1)) {
  62                 x = (x + M_PI_4) + M_PI_4;
  63         }
  64         else if (cosfl && negative && i == 1) {
  65                 x = (x + M_PI_2) + M_PI_2;
  66         }
  67         if (! cosfl && negative) x = -x;
  68         return x;
  69 }
  70
  71 double
  72 asin(x)
  73         double x;
  74 {
  75         return asin_acos(x, 0);
  76 }
  77
  78 double
  79 acos(x)
  80         double x;
  81 {
  82         return asin_acos(x, 1);
  83 }
  84 \0atan2.c\0\0\0\0\0\0\0\0\0\0\0\ 2\ 2�\ 1\0\0e\ 3/*
  85  * (c) copyright 1988 by the Vrije Universiteit, Amsterdam, The Netherlands.
  86  * See the copyright notice in the ACK home directory, in the file "Copyright".
  87  *
  88  * Author: Ceriel J.H. Jacobs
  89  */
  90
  91 /* $Id: atan2.c,v 1.2 1994/06/24 12:22:36 ceriel Exp $ */
  92
  93 #include <math.h>
  94 #include <errno.h>
  95
  96 extern int errno;
  97
  98 double
  99 atan2(y, x)
 100         double x, y;
 101 {
 102         extern double atan();
 103         double absx, absy, val;
 104
 105         if (x == 0 && y == 0) {
 106                 errno = EDOM;
 107                 return 0;
 108         }
 109         absy = y < 0 ? -y : y;
 110         absx = x < 0 ? -x : x;
 111         if (absy - absx == absy) {
 112                 /* x negligible compared to y */
 113                 return y < 0 ? -M_PI_2 : M_PI_2;
 114         }
 115         if (absx - absy == absx) {
 116                 /* y negligible compared to x */
 117                 val = 0.0;
 118         }
 119         else    val = atan(y/x);
 120         if (x > 0) {
 121                 /* first or fourth quadrant; already correct */
 122                 return val;
 123         }
 124         if (y < 0) {
 125                 /* third quadrant */
 126                 return val - M_PI;
 127         }
 128         return val + M_PI;
 129 }
 130 Oatan.c\0\0\0\0\0\0\0\0\0\0\0\0\ 2\ 2�\ 1\0\0?\ 5/*
 131  * (c) copyright 1988 by the Vrije Universiteit, Amsterdam, The Netherlands.
 132  * See the copyright notice in the ACK home directory, in the file "Copyright".
 133  *
 134  * Author: Ceriel J.H. Jacobs
 135  */
 136
 137 /* $Id: atan.c,v 1.3 1994/06/24 12:22:33 ceriel Exp $ */
 138
 139 #include <math.h>
 140 #include <errno.h>
 141
 142 extern int errno;
 143
 144 double
 145 atan(x)
 146         double x;
 147 {
 148         /*      Algorithm and coefficients from:
 149                         "Software manual for the elementary functions"
 150                         by W.J. Cody and W. Waite, Prentice-Hall, 1980
 151         */
 152
 153         static double p[] = {
 154                 -0.13688768894191926929e+2,
 155                 -0.20505855195861651981e+2,
 156                 -0.84946240351320683534e+1,
 157                 -0.83758299368150059274e+0
 158         };
 159         static double q[] = {
 160                  0.41066306682575781263e+2,
 161                  0.86157349597130242515e+2,
 162                  0.59578436142597344465e+2,
 163                  0.15024001160028576121e+2,
 164                  1.0
 165         };
 166         static double a[] = {
 167                 0.0,
 168                 0.52359877559829887307710723554658381,  /* pi/6 */
 169                 M_PI_2,
 170                 1.04719755119659774615421446109316763   /* pi/3 */
 171         };
 172
 173         int     neg = x < 0;
 174         int     n;
 175         double  g;
 176
 177         if (neg) {
 178                 x = -x;
 179         }
 180         if (x > 1.0) {
 181                 x = 1.0/x;
 182                 n = 2;
 183         }
 184         else    n = 0;
 185
 186         if (x > 0.26794919243112270647) {       /* 2-sqtr(3) */
 187                 n = n + 1;
 188                 x = (((0.73205080756887729353*x-0.5)-0.5)+x)/
 189                         (1.73205080756887729353+x);
 190         }
 191
 192         /* ??? avoid underflow ??? */
 193
 194         g = x * x;
 195         x += x * g * POLYNOM3(g, p) / POLYNOM4(g, q);
 196         if (n > 1) x = -x;
 197         x += a[n];
 198         return neg ? -x : x;
 199 }
 200 eceil.c\0\0\0\0\0\0\0\0\0\0\0\0\ 2\ 2�\ 1\0\0�\ 1/*
 201  * (c) copyright 1988 by the Vrije Universiteit, Amsterdam, The Netherlands.
 202  * See the copyright notice in the ACK home directory, in the file "Copyright".
 203  *
 204  * Author: Ceriel J.H. Jacobs
 205  */
 206
 207 /* $Id: ceil.c,v 1.2 1994/06/24 12:22:39 ceriel Exp $ */
 208
 209 double
 210 ceil(x)
 211         double x;
 212 {
 213         extern double modf();
 214         double val;
 215
 216         return modf(x, &val) > 0 ? val + 1.0 : val ;
 217         /*      this also works if modf always returns a positive
 218                 fractional part
 219         */
 220 }
 221 nfabs.c\0\0\0\0\0\0\0\0\0\0\0\0\ 2\ 2�\ 1\0\07\ 1/*
 222  * (c) copyright 1988 by the Vrije Universiteit, Amsterdam, The Netherlands.
 223  * See the copyright notice in the ACK home directory, in the file "Copyright".
 224  *
 225  * Author: Ceriel J.H. Jacobs
 226  */
 227
 228 /* $Id: fabs.c,v 1.2 1994/06/24 12:22:45 ceriel Exp $ */
 229
 230 double
 231 fabs(x)
 232         double x;
 233 {
 234         return  x < 0 ? -x : x;
 235 }
 236 bgamma.c\0\0\0\0\0\0\0\0\0\0\0\ 2\ 2�\ 1\0\0�\v/*
 237  * (c) copyright 1988 by the Vrije Universiteit, Amsterdam, The Netherlands.
 238  * See the copyright notice in the ACK home directory, in the file "Copyright".
 239  *
 240  * Author: Ceriel J.H. Jacobs
 241  */
 242
 243 /* $Id: gamma.c,v 1.3 1994/06/24 12:22:51 ceriel Exp $ */
 244
 245 #include <math.h>
 246 #include <errno.h>
 247
 248 extern int errno;
 249
 250 static double
 251 smallpos_gamma(x)
 252         double x;
 253 {
 254         /*      Approximation of gamma function using
 255                 gamma(x) = P(x-2) / Q(x-2) for x in [2,3]
 256         */
 257         /* Hart & Cheney # 5251 */
 258
 259         static double p[11] = {
 260                 -0.2983543278574342138830437659e+06,
 261                 -0.2384953970018198872468734423e+06,
 262                 -0.1170494760121780688403854445e+06,
 263                 -0.3949445048301571936421824091e+05,
 264                 -0.1046699423827521405330650531e+05,
 265                 -0.2188218110071816359394795998e+04,
 266                 -0.3805112208641734657584922631e+03,
 267                 -0.5283123755635845383718978382e+02,
 268                 -0.6128571763704498306889428212e+01,
 269                 -0.5028018054416812467364198750e+00,
 270                 -0.3343060322330595274515660112e-01
 271         };
 272
 273         static double q[9] = {
 274                 -0.2983543278574342138830438524e+06,
 275                 -0.1123558608748644911342306408e+06,
 276                  0.5332716689118142157485686311e+05,
 277                  0.8571160498907043851961147763e+04,
 278                 -0.4734865977028211706556819770e+04,
 279                  0.1960497612885585838997039621e+03,
 280                  0.1257733367869888645966647426e+03,
 281                 -0.2053126153100672764513929067e+02,
 282                  0.1000000000000000000000000000e+01
 283         };
 284
 285         double result = 1.0;
 286
 287         while (x > 3) {
 288                 x -= 1.0;
 289                 result *= x;
 290         }
 291         while (x < 2) {
 292                 result /= x;
 293                 x += 1.0;
 294         }
 295
 296         x -= 2.0;
 297
 298         return result * POLYNOM10(x, p) / POLYNOM8(x, q);
 299 }
 300
 301 #define log_sqrt_2pi 0.91893853320467274178032973640561763
 302
 303 int     signgam;
 304
 305 static double
 306 bigpos_loggamma(x)
 307         double x;
 308 {
 309         /*      computes the log(gamma(x)) function for big arguments
 310                 using the Stirling form
 311                   log(gamma(x)) = (x - 0.5)log(x) - x + log(sqrt(2*pi)) + fi(x)
 312                 where fi(x) = (1/x)*P(1/(x*x))/Q(1/(x*x)) for x in [12,1000]
 313         */
 314         /* Hart & Cheney # 5468 */
 315
 316         static double p[4] = {
 317                  0.12398282342474941538685913e+00,
 318                  0.67082783834332134961461700e+00,
 319                  0.64507302912892202513890000e+00,
 320                  0.66662907040200752600000000e-01
 321         };
 322
 323         static double q[4] = {
 324                  0.14877938810969929846815600e+01,
 325                  0.80995271894897557472821400e+01,
 326                  0.79966911236636441947720000e+01,
 327                  0.10000000000000000000000000e+01
 328         };
 329
 330         double rsq = 1.0/(x*x);
 331         extern double log();
 332
 333         return (x-0.5)*log(x)-x+log_sqrt_2pi+POLYNOM3(rsq, p)/(x*POLYNOM3(rsq, q));
 334 }
 335
 336 static double
 337 neg_loggamma(x)
 338         double x;
 339 {
 340         /*      compute the log(gamma(x)) function for negative values of x,
 341                 using the rule:
 342                         -x*gamma(x)*gamma(-x) = pi/sin(z*pi)
 343         */
 344         extern double sin(), log();
 345         double sinpix;
 346
 347         x = -x;
 348         sinpix = sin(M_PI * x);
 349         if (sinpix == 0.0) {
 350                 errno = EDOM;
 351                 return HUGE;
 352         }
 353         if (sinpix < 0) sinpix = -sinpix;
 354         else signgam = -1;
 355         return log(M_PI/(x * smallpos_gamma(x) * sinpix));
 356 }
 357
 358 double
 359 gamma(x)
 360         double x;
 361 {
 362         /*      Wrong name; Actually computes log(gamma(x))
 363         */
 364         extern double log();
 365
 366         signgam = 1;
 367         if (x <= 0) {
 368                 return neg_loggamma(x);
 369         }
 370         if (x > 12.0) {
 371                 return bigpos_loggamma(x);
 372         }
 373         return log(smallpos_gamma(x));
 374 }
 375 hypot.c\0\0\0\0\0\0\0\0\0\0\0\ 2\ 2�\ 1\0\0�\ 2/*
 376  * (c) copyright 1988 by the Vrije Universiteit, Amsterdam, The Netherlands.
 377  * See the copyright notice in the ACK home directory, in the file "Copyright".
 378  *
 379  * Author: Ceriel J.H. Jacobs
 380  */
 381
 382 /* $Id: hypot.c,v 1.2 1994/06/24 12:22:54 ceriel Exp $ */
 383
 384 double
 385 hypot(x,y)
 386         double x,y;
 387 {
 388         /*      Computes sqrt(x*x+y*y), avoiding overflow */
 389
 390         extern double sqrt();
 391
 392         if (x < 0) x = -x;
 393         if (y < 0) y = -y;
 394         if (x > y) {
 395                 double t = y;
 396                 y = x;
 397                 x = t;
 398         }
 399         /* sqrt(x*x+y*y) = sqrt(y*y*(x*x/(y*y)+1.0)) = y*sqrt(x*x/(y*y)+1.0) */
 400         x /= y;
 401         return y*sqrt(x*x+1.0);
 402 }
 403
 404 struct complex {
 405         double r,i;
 406 };
 407
 408 double
 409 cabs(p_compl)
 410         struct complex p_compl;
 411 {
 412         return hypot(p_compl.r, p_compl.i);
 413 }
 414 jn.c\0.c\0\0\0\0\0\0\0\0\0\0\0\ 2\ 2�\ 1\0\0\ 2
 415 /*
 416  * (c) copyright 1988 by the Vrije Universiteit, Amsterdam, The Netherlands.
 417  * See the copyright notice in the ACK home directory, in the file "Copyright".
 418  *
 419  * Author: Ceriel J.H. Jacobs
 420  */
 421
 422 /* $Id: jn.c,v 1.3 1994/06/24 12:23:03 ceriel Exp $ */
 423
 424 #include <math.h>
 425 #include <errno.h>
 426
 427 extern int errno;
 428
 429 double
 430 yn(n, x)
 431         double x;
 432 {
 433         /*      Use y0, y1, and the recurrence relation
 434                 y(n+1,x) = 2*n*y(n,x)/x - y(n-1, x).
 435                 According to Hart & Cheney, this is stable for all
 436                 x, n.
 437                 Also use: y(-n,x) = (-1)^n * y(n, x)
 438         */
 439
 440         int negative = 0;
 441         extern double y0(), y1();
 442         double yn1, yn2;
 443         register int i;
 444
 445         if (x <= 0) {
 446                 errno = EDOM;
 447                 return -HUGE;
 448         }
 449
 450         if (n < 0) {
 451                 n = -n;
 452                 negative = (n % 2);
 453         }
 454
 455         if (n == 0) return y0(x);
 456         if (n == 1) return y1(x);
 457
 458         yn2 = y0(x);
 459         yn1 = y1(x);
 460         for (i = 1; i < n; i++) {
 461                 double tmp = yn1;
 462                 yn1 = (i*2)*yn1/x - yn2;
 463                 yn2 = tmp;
 464         }
 465         if (negative) return -yn1;
 466         return yn1;
 467 }
 468
 469 double
 470 jn(n, x)
 471         double x;
 472 {
 473         /*      Unfortunately, according to Hart & Cheney, the recurrence
 474                 j(n+1,x) = 2*n*j(n,x)/x - j(n-1,x) is unstable for
 475                 increasing n, except when x > n.
 476                 However, j(n,x)/j(n-1,x) = 2/(2*n-x*x/(2*(n+1)-x*x/( ....
 477                 (a continued fraction).
 478                 We can use this to determine KJn and KJn-1, where K is a
 479                 normalization constant not yet known. This enables us
 480                 to determine KJn-2, ...., KJ1, KJ0. Now we can use the
 481                 J0 or J1 approximation to determine K.
 482                 Use: j(-n, x) = (-1)^n * j(n, x)
 483                      j(n, -x) = (-1)^n * j(n, x)
 484         */
 485
 486         extern double j0(), j1();
 487
 488         if (n < 0) {
 489                 n = -n;
 490                 x = -x;
 491         }
 492
 493         if (n == 0) return j0(x);
 494         if (n == 1) return j1(x);
 495         if (x > n) {
 496                 /* in this case, the recurrence relation is stable for
 497                    increasing n, so we use that.
 498                 */
 499                 double jn2 = j0(x), jn1 = j1(x);
 500                 register int i;
 501
 502                 for (i = 1; i < n; i++) {
 503                         double tmp = jn1;
 504                         jn1 = (2*i)*jn1/x - jn2;
 505                         jn2 = tmp;
 506                 }
 507                 return jn1;
 508         }
 509         {
 510                 /* we first compute j(n,x)/j(n-1,x) */
 511                 register int i;
 512                 double quotient = 0.0;
 513                 double xsqr = x*x;
 514                 double jn1, jn2;
 515
 516                 for (i = 20;    /* ??? how many do we need ??? */
 517                      i > 0; i--) {
 518                         quotient = xsqr/(2*(i+n) - quotient);
 519                 }
 520                 quotient = x / (2*n - quotient);
 521
 522                 jn1 = quotient;
 523                 jn2 = 1.0;
 524                 for (i = n-1; i > 0; i--) {
 525                         /* recurrence relation is stable for decreasing n
 526                         */
 527                         double tmp = jn2;
 528                         jn2 = (2*i)*jn2/x - jn1;
 529                         jn1 = tmp;
 530                 }
 531                 /* So, now we have K*Jn = quotient and K*J0 = jn2.
 532                    Now it is easy; compute real j0, this gives K = jn2/j0,
 533                    and this then gives Jn = quotient/K = j0 * quotient / jn2.
 534                 */
 535                 return j0(x)*quotient/jn2;
 536         }
 537 }
 538 j0.c\0.c\0\0\0\0\0\0\0\0\0\0\0\ 2\ 2�\ 1\0\0�\13/*
 539  * (c) copyright 1988 by the Vrije Universiteit, Amsterdam, The Netherlands.
 540  * See the copyright notice in the ACK home directory, in the file "Copyright".
 541  *
 542  * Author: Ceriel J.H. Jacobs
 543  */
 544
 545 /* $Id: j0.c,v 1.3 1994/06/24 12:22:57 ceriel Exp $ */
 546
 547 #include <math.h>
 548 #include <errno.h>
 549
 550 extern int errno;
 551
 552 static double
 553 P0(x)
 554         double x;
 555 {
 556         /*      P0(x) = P(z*z)/Q(z*z) where z = 8/x, with x >= 8 */
 557         /*      Hart & Cheney # 6554 */
 558
 559         static double p[9] = {
 560                  0.9999999999999999999999995647e+00,
 561                  0.5638253933310769952531889297e+01,
 562                  0.1124846237418285392887270013e+02,
 563                  0.1009280644639441488899111404e+02,
 564                  0.4290591487686900980651458361e+01,
 565                  0.8374209971661497198619102718e+00,
 566                  0.6702347074465611456598882534e-01,
 567                  0.1696260729396856143084502774e-02,
 568                  0.6463970103128382090713889584e-05
 569         };
 570
 571         static double q[9] = {
 572                  0.9999999999999999999999999999e+00,
 573                  0.5639352566123269952531467562e+01,
 574                  0.1125463057106955935416066535e+02,
 575                  0.1010501892629524191262518048e+02,
 576                  0.4301396985171094350444425443e+01,
 577                  0.8418926086780046799127094223e+00,
 578                  0.6784915305473610998681570734e-01,
 579                  0.1754416614608056207958880988e-02,
 580                  0.7482977995134121064747276923e-05
 581         };
 582
 583         double zsq = 64.0/(x*x);
 584
 585         return POLYNOM8(zsq, p) / POLYNOM8(zsq, q);
 586 }
 587
 588 static double
 589 Q0(x)
 590         double x;
 591 {
 592         /*      Q0(x) = z*P(z*z)/Q(z*z) where z = 8/x, x >= 8 */
 593         /*      Hart & Cheney # 6955 */
 594         /*      Probably typerror in Hart & Cheney; it sais:
 595                 Q0(x) = x*P(z*z)/Q(z*z)
 596         */
 597
 598         static double p[9] = {
 599                 -0.1562499999999999999999995808e-01,
 600                 -0.1111285583113679178917024959e+00,
 601                 -0.2877685516355036842789761274e+00,
 602                 -0.3477683453166454475665803194e+00,
 603                 -0.2093031978191084473537206358e+00,
 604                 -0.6209520943730206312601003832e-01,
 605                 -0.8434508346572023650653353729e-02,
 606                 -0.4414848186188819989871882393e-03,
 607                 -0.5768946278415631134804064871e-05
 608         };
 609
 610         static double q[10] = {
 611                  0.9999999999999999999999999999e+00,
 612                  0.7121383005365046745065850254e+01,
 613                  0.1848194194302368046679068851e+02,
 614                  0.2242327522435983712994071530e+02,
 615                  0.1359286169255959339963319677e+02,
 616                  0.4089489268101204780080944780e+01,
 617                  0.5722140925672174525430730669e+00,
 618                  0.3219814230905924725810683346e-01,
 619                  0.5299687475496044642364124073e-03,
 620                  0.9423249021001925212258428217e-06
 621         };
 622
 623         double zsq = 64.0/(x*x);
 624
 625         return (8.0/x) * POLYNOM8(zsq, p) / POLYNOM9(zsq, q);
 626 }
 627
 628 static double
 629 smallj0(x)
 630         double x;
 631 {
 632         /*      J0(x) = P(x*x)/Q(x*x) for x in [0,8] */
 633         /*      Hart & Cheney # 5852 */
 634
 635         static double p[10] = {
 636                  0.1641556014884554385346147435e+25,
 637                 -0.3943559664767296636012616471e+24,
 638                  0.2172018385924539313982287997e+23,
 639                 -0.4814859952069817648285245941e+21,
 640                  0.5345457598841972345381674607e+19,
 641                 -0.3301538925689637686465426220e+17,
 642                  0.1187390681211042949874031474e+15,
 643                 -0.2479851167896144439689877514e+12,
 644                  0.2803148940831953934479400118e+09,
 645                 -0.1336625500481224741885945416e+06
 646         };
 647
 648         static double q[10] = {
 649                  0.1641556014884554385346137617e+25,
 650                  0.1603303724440893273539045602e+23,
 651                  0.7913043777646405204323616203e+20,
 652                  0.2613165313325153278086066185e+18,
 653                  0.6429607918826017759289213100e+15,
 654                  0.1237672982083407903483177730e+13,
 655                  0.1893012093677918995179541438e+10,
 656                  0.2263381356781110003609399116e+07,
 657                  0.1974019272727281783930443513e+04,
 658                  0.1000000000000000000000000000e+01
 659         };
 660
 661         double xsq = x*x;
 662
 663         return POLYNOM9(xsq, p) / POLYNOM9(xsq, q);
 664 }
 665
 666 double
 667 j0(x)
 668         double x;
 669 {
 670         /*      Use J0(x) = sqrt(2/(pi*x))*(P0(x)*cos(X0)-Q0(x)*sin(X0))
 671                         where X0 = x - pi/4 for |x| > 8.
 672                 Use J0(-x) = J0(x).
 673                 Use direct approximation of smallj0 for |x| <= 8.
 674         */
 675         extern double sqrt(), sin(), cos();
 676
 677         if (x < 0) x = -x;
 678         if (x > 8.0) {
 679                 double X0 = x - M_PI_4;
 680                 return sqrt(M_2_PI/x)*(P0(x)*cos(X0) - Q0(x)*sin(X0));
 681         }
 682         return smallj0(x);
 683 }
 684
 685 static double
 686 smally0_bar(x)
 687         double x;
 688 {
 689         /*      Y0(x) = Y0BAR(x)+(2/pi)*J0(x)ln(x)
 690                 Approximation of Y0BAR for 0 <= x <= 8:
 691                         Y0BAR(x) = P(x*x)/Q(x*x)
 692                 Hart & Cheney #6250
 693         */
 694
 695         static double p[14] = {
 696                 -0.2692670958801060448840356941e+14,
 697                  0.6467231173109037044444917683e+14,
 698                 -0.5563036156275660297303897296e+13,
 699                  0.1698403391975239335187832821e+12,
 700                 -0.2606282788256139370857687880e+10,
 701                  0.2352841334491277505699488812e+08,
 702                 -0.1365184412186963659690851354e+06,
 703                  0.5371538422626582142170627457e+03,
 704                 -0.1478903875146718839145348490e+01,
 705                  0.2887840299886172125955719069e-02,
 706                 -0.3977426824263991024666116123e-05,
 707                  0.3738169731655229006655176866e-08,
 708                 -0.2194460874896856106887900645e-11,
 709                  0.6208996973821484304384239393e-15
 710         };
 711
 712         static double q[6] = {
 713                  0.3648393301278364629844168660e+15,
 714                  0.1698390180526960997295118328e+13,
 715                  0.3587111679107612117789088586e+10,
 716                  0.4337760840406994515845890005e+07,
 717                  0.3037977771964348276793136205e+04,
 718                  0.1000000000000000000000000000e+01
 719         };
 720
 721         double xsq = x*x;
 722
 723         return POLYNOM13(xsq, p) / POLYNOM5(xsq, q);
 724 }
 725
 726 double
 727 y0(x)
 728         double x;
 729 {
 730         extern double sqrt(), sin(), cos(), log();
 731
 732         if (x <= 0.0) {
 733                 errno = EDOM;
 734                 return -HUGE;
 735         }
 736         if (x > 8.0) {
 737                 double X0 = x - M_PI_4;
 738                 return sqrt(M_2_PI/x) * (P0(x)*sin(X0)+Q0(x)*cos(X0));
 739         }
 740         return smally0_bar(x) + M_2_PI*j0(x)*log(x);
 741 }
 742 \0j1.c\0.c\0\0\0\0\0\0\0\0\0\0\0\ 2\ 2�\ 1\0\0H\14/*
 743  * (c) copyright 1988 by the Vrije Universiteit, Amsterdam, The Netherlands.
 744  * See the copyright notice in the ACK home directory, in the file "Copyright".
 745  *
 746  * Author: Ceriel J.H. Jacobs
 747  */
 748
 749 /* $Id: j1.c,v 1.3 1994/06/24 12:23:00 ceriel Exp $ */
 750
 751 #include <math.h>
 752 #include <errno.h>
 753
 754 extern int errno;
 755
 756 static double
 757 P1(x)
 758         double x;
 759 {
 760         /*      P1(x) = P(z*z)/Q(z*z) where z = 8/x, with x >= 8 */
 761         /*      Hart & Cheney # 6755 */
 762
 763         static double p[9] = {
 764                  0.1000000000000000000000000489e+01,
 765                  0.5581663300347182292169450071e+01,
 766                  0.1100186625131173123750501118e+02,
 767                  0.9727139359130463694593683431e+01,
 768                  0.4060011483142278994462590992e+01,
 769                  0.7742832212665311906917358099e+00,
 770                  0.6021617752811098752098248630e-01,
 771                  0.1482350677236405118074646993e-02,
 772                  0.6094215148131061431667573909e-05
 773         };
 774
 775         static double q[9] = {
 776                  0.9999999999999999999999999999e+00,
 777                  0.5579832245659682292169922224e+01,
 778                  0.1099168447731617288972771040e+02,
 779                  0.9707206835125961446797916892e+01,
 780                  0.4042610016540342097334497865e+01,
 781                  0.7671965204303836019508430169e+00,
 782                  0.5893258668794493100786371406e-01,
 783                  0.1393993644981256852404222530e-02,
 784                  0.4585597769784750669754696825e-05
 785         };
 786
 787         double zsq = 64.0/(x*x);
 788
 789         return POLYNOM8(zsq, p) / POLYNOM8(zsq, q);
 790 }
 791
 792 static double
 793 Q1(x)
 794         double x;
 795 {
 796         /*      Q1(x) = z*P(z*z)/Q(z*z) where z = 8/x, x >= 8 */
 797         /*      Hart & Cheney # 7157 */
 798         /*      Probably typerror in Hart & Cheney; it sais:
 799                 Q1(x) = x*P(z*z)/Q(z*z)
 800         */
 801
 802         static double p[9] = {
 803                 0.4687499999999999999999995275e-01,
 804                 0.3302394516691663879252493748e+00,
 805                 0.8456888491208195767613862428e+00,
 806                 0.1008551084218946085420665147e+01,
 807                 0.5973407972399900690521296181e+00,
 808                 0.1737697433393258207540273097e+00,
 809                 0.2303862814819568573893610740e-01,
 810                 0.1171224207976250587945594946e-02,
 811                 0.1486418220337492918307904804e-04
 812         };
 813
 814         static double q[10] = {
 815                 0.9999999999999999999999999999e+00,
 816                 0.7049380763213049609070823421e+01,
 817                 0.1807129960468949760845562209e+02,
 818                 0.2159171174362827330505421695e+02,
 819                 0.1283239297740546866114600499e+02,
 820                 0.3758349275324260869598403931e+01,
 821                 0.5055985453754739528620657666e+00,
 822                 0.2665604326323907148063400439e-01,
 823                 0.3821140353404633025596424652e-03,
 824                 0.3206696590241261037875154062e-06
 825         };
 826
 827         double zsq = 64.0/(x*x);
 828
 829         return (8.0/x) * POLYNOM8(zsq, p) / POLYNOM9(zsq, q);
 830 }
 831
 832 static double
 833 smallj1(x)
 834         double x;
 835 {
 836         /*      J1(x) = x*P(x*x)/Q(x*x) for x in [0,8] */
 837         /*      Hart & Cheney # 6054 */
 838
 839         static double p[10] = {
 840                  0.1921176307760798128049021316e+25,
 841                 -0.2226092031387396254771375773e+24,
 842                  0.7894463902082476734673226741e+22,
 843                 -0.1269424373753606065436561036e+21,
 844                  0.1092152214043184787101134641e+19,
 845                 -0.5454629264396819144157448868e+16,
 846                  0.1634659487571284628830445048e+14,
 847                 -0.2909662785381647825756152444e+11,
 848                  0.2853433451054763915026471449e+08,
 849                 -0.1197705712815379389149134705e+05
 850         };
 851
 852         static double q[10] = {
 853                  0.3842352615521596256098041912e+25,
 854                  0.3507567066272028105798868716e+23,
 855                  0.1611334311633414344007062889e+21,
 856                  0.4929612313959850319632645381e+18,
 857                  0.1117536965288162684489793105e+16,
 858                  0.1969278625584719037168592923e+13,
 859                  0.2735606122949877990248154504e+10,
 860                  0.2940957355049651347475558106e+07,
 861                  0.2274736606126590905134610965e+04,
 862                  0.1000000000000000000000000000e+01
 863         };
 864
 865         double xsq = x*x;
 866
 867         return x * POLYNOM9(xsq, p) / POLYNOM9(xsq, q);
 868 }
 869
 870 double
 871 j1(x)
 872         double x;
 873 {
 874         /*      Use J1(x) = sqrt(2/(pi*x))*(P1(x)*cos(X1)-Q1(x)*sin(X1))
 875                         where X1 = x - 3*pi/4 for |x| > 8.
 876                 Use J1(-x) = -J1(x).
 877                 Use direct approximation of smallj1 for |x| <= 8.
 878         */
 879         extern double sqrt(), sin(), cos();
 880         int negative = x < 0.0;
 881
 882         if (negative) x = -x;
 883         if (x > 8.0) {
 884                 double X1 = x - (M_PI - M_PI_4);
 885                 x = sqrt(M_2_PI/x)*(P1(x)*cos(X1) - Q1(x)*sin(X1));
 886         }
 887         else x = smallj1(x);
 888         if (negative) return -x;
 889         return x;
 890 }
 891
 892 static double
 893 smally1_bar(x)
 894         double x;
 895 {
 896         /*      Y1(x) = Y1BAR(x)+(2/pi)*(J1(x)ln(x) - 1/x)
 897                 Approximation of Y1BAR for 0 <= x <= 8:
 898                         Y1BAR(x) = x*P(x*x)/Q(x*x)
 899                 Hart & Cheney # 6449
 900         */
 901
 902         static double p[10] = {
 903                 -0.5862655424363443992938931700e+24,
 904                  0.1570668341992328458208364904e+24,
 905                 -0.7351681299005467428400402479e+22,
 906                  0.1390658785759080111485190942e+21,
 907                 -0.1339544201526785345938109179e+19,
 908                  0.7290257386242270629526344379e+16,
 909                 -0.2340575603057015935501295099e+14,
 910                  0.4411516199185230690878878903e+11,
 911                 -0.4542128738770213026987060358e+08,
 912                  0.1988612563465350530472715888e+05
 913         };
 914
 915         static double q[10] = {
 916                  0.2990279721605116022908679994e+25,
 917                  0.2780285010357803058127175655e+23,
 918                  0.1302687474507355553192845146e+21,
 919                  0.4071330372239164349602952937e+18,
 920                  0.9446611865086570116528399283e+15,
 921                  0.1707657951197456205887347694e+13,
 922                  0.2440358986882941823431612517e+10,
 923                  0.2708852767034077697963790196e+07,
 924                  0.2174361138333330803617969305e+04,
 925                  0.1000000000000000000000000000e+01
 926         };
 927
 928         double xsq = x*x;
 929
 930         return x * POLYNOM9(xsq, p) / POLYNOM9(xsq, q);
 931 }
 932
 933 double
 934 y1(x)
 935         double x;
 936 {
 937         extern double sqrt(), sin(), cos(), log();
 938
 939         if (x <= 0.0) {
 940                 errno = EDOM;
 941                 return -HUGE;
 942         }
 943         if (x > 8.0) {
 944                 double X1 = x - (M_PI - M_PI_4);
 945                 return sqrt(M_2_PI/x) * (P1(x)*sin(X1)+Q1(x)*cos(X1));
 946         }
 947         return smally1_bar(x) + M_2_PI*(j1(x)*log(x) - 1/x);
 948 }
 949 log10.c\0\0\0\0\0\0\0\0\0\0\0\ 2\ 2�\ 1\0\0�\ 1/*
 950  * (c) copyright 1988 by the Vrije Universiteit, Amsterdam, The Netherlands.
 951  * See the copyright notice in the ACK home directory, in the file "Copyright".
 952  *
 953  * Author: Ceriel J.H. Jacobs
 954  */
 955
 956 /* $Id: log10.c,v 1.2 1994/06/24 12:23:08 ceriel Exp $ */
 957
 958 #include <math.h>
 959 #include <errno.h>
 960
 961 extern int errno;
 962
 963 double
 964 log10(x)
 965         double x;
 966 {
 967         extern double log();
 968
 969         if (x <= 0) {
 970                 errno = EDOM;
 971                 return 0;
 972         }
 973
 974         return log(x) / M_LN10;
 975 }
 976 pow.c\0c\0\0\0\0\0\0\0\0\0\0\0\ 2\ 2�\ 1\0\0.\ 4/*
 977  * (c) copyright 1988 by the Vrije Universiteit, Amsterdam, The Netherlands.
 978  * See the copyright notice in the ACK home directory, in the file "Copyright".
 979  *
 980  * Author: Ceriel J.H. Jacobs
 981  */
 982
 983 /* $Id: pow.c,v 1.3 1994/06/24 12:23:11 ceriel Exp $ */
 984
 985 #include <math.h>
 986 #include <errno.h>
 987
 988 extern int errno;
 989 extern double modf(), exp(), log();
 990
 991 double
 992 pow(x,y)
 993         double x,y;
 994 {
 995         /*      Simple version for now. The Cody and Waite book has
 996                 a very complicated, much more precise version, but
 997                 this version has machine-dependant arrays A1 and A2,
 998                 and I don't know yet how to solve this ???
 999         */
1000         double dummy;
1001         int     result_neg = 0;
1002
1003         if ((x == 0 && y == 0) ||
1004             (x < 0 && modf(y, &dummy) != 0)) {
1005                 errno = EDOM;
1006                 return 0;
1007         }
1008
1009         if (x == 0) return x;
1010
1011         if (x < 0) {
1012                 if (modf(y/2.0, &dummy) != 0) {
1013                         /* y was odd */
1014                         result_neg = 1;
1015                 }
1016                 x = -x;
1017         }
1018         x = log(x);
1019         if (x < 0) {
1020                 x = -x;
1021                 y = -y;
1022         }
1023         if (y > M_LN_MAX_D/x) {
1024                 errno = ERANGE;
1025                 return 0;
1026         }
1027         if (y < M_LN_MIN_D/x) {
1028                 errno = ERANGE;
1029                 return 0;
1030         }
1031
1032         x = exp(x * y);
1033         return result_neg ? -x : x;
1034 }
1035 log.c\0c\0\0\0\0\0\0\0\0\0\0\0\ 2\ 2�\ 1\0\0\88\ 4/*
1036  * (c) copyright 1989 by the Vrije Universiteit, Amsterdam, The Netherlands.
1037  * See the copyright notice in the ACK home directory, in the file "Copyright".
1038  *
1039  * Author: Ceriel J.H. Jacobs
1040  */
1041
1042 /* $Id: log.c,v 1.3 1994/06/24 12:23:05 ceriel Exp $ */
1043
1044 #include <math.h>
1045 #include <errno.h>
1046
1047 extern int      errno;
1048 extern double   frexp();
1049
1050 double
1051 log(x)
1052         double  x;
1053 {
1054         /*      Algorithm and coefficients from:
1055                         "Software manual for the elementary functions"
1056                         by W.J. Cody and W. Waite, Prentice-Hall, 1980
1057         */
1058         static double a[] = {
1059                 -0.64124943423745581147e2,
1060                  0.16383943563021534222e2,
1061                 -0.78956112887491257267e0
1062         };
1063         static double b[] = {
1064                 -0.76949932108494879777e3,
1065                  0.31203222091924532844e3,
1066                 -0.35667977739034646171e2,
1067                  1.0
1068         };
1069
1070         double  znum, zden, z, w;
1071         int     exponent;
1072
1073         if (x <= 0) {
1074                 errno = EDOM;
1075                 return 0;
1076         }
1077
1078         x = frexp(x, &exponent);
1079         if (x > M_1_SQRT2) {
1080                 znum = (x - 0.5) - 0.5;
1081                 zden = x * 0.5 + 0.5;
1082         }
1083         else {
1084                 znum = x - 0.5;
1085                 zden = znum * 0.5 + 0.5;
1086                 exponent--;
1087         }
1088         z = znum/zden; w = z * z;
1089         x = z + z * w * (POLYNOM2(w,a)/POLYNOM3(w,b));
1090         z = exponent;
1091         x += z * (-2.121944400546905827679e-4);
1092         return x + z * 0.693359375;
1093 }
1094 sin.c\0c\0\0\0\0\0\0\0\0\0\0\0\ 2\ 2�\ 1\0\0�\ 6/*
1095  * (c) copyright 1988 by the Vrije Universiteit, Amsterdam, The Netherlands.
1096  * See the copyright notice in the ACK home directory, in the file "Copyright".
1097  *
1098  * Author: Ceriel J.H. Jacobs
1099  */
1100
1101 /* $Id: sin.c,v 1.3 1994/06/24 12:23:14 ceriel Exp $ */
1102
1103 #include <math.h>
1104 #include <errno.h>
1105
1106 extern int      errno;
1107 extern double   modf();
1108
1109 static double
1110 sinus(x, cos_flag)
1111         double x;
1112 {
1113         /*      Algorithm and coefficients from:
1114                         "Software manual for the elementary functions"
1115                         by W.J. Cody and W. Waite, Prentice-Hall, 1980
1116         */
1117
1118         static double r[] = {
1119                 -0.16666666666666665052e+0,
1120                  0.83333333333331650314e-2,
1121                 -0.19841269841201840457e-3,
1122                  0.27557319210152756119e-5,
1123                 -0.25052106798274584544e-7,
1124                  0.16058936490371589114e-9,
1125                 -0.76429178068910467734e-12,
1126                  0.27204790957888846175e-14
1127         };
1128
1129         double  xsqr;
1130         double  y;
1131         int     neg = 0;
1132
1133         if (x < 0) {
1134                 x = -x;
1135                 neg = 1;
1136         }
1137         if (cos_flag) {
1138                 neg = 0;
1139                 y = M_PI_2 + x;
1140         }
1141         else    y = x;
1142
1143         /* ??? avoid loss of significance, if y is too large, error ??? */
1144
1145         y = y * M_1_PI + 0.5;
1146
1147         /*      Use extended precision to calculate reduced argument.
1148                 Here we used 12 bits of the mantissa for a1.
1149                 Also split x in integer part x1 and fraction part x2.
1150         */
1151 #define A1 3.1416015625
1152 #define A2 -8.908910206761537356617e-6
1153         {
1154                 double x1, x2;
1155
1156                 modf(y, &y);
1157                 if (modf(0.5*y, &x1)) neg = !neg;
1158                 if (cos_flag) y -= 0.5;
1159                 x2 = modf(x, &x1);
1160                 x = x1 - y * A1;
1161                 x += x2;
1162                 x -= y * A2;
1163 #undef A1
1164 #undef A2
1165         }
1166
1167         if (x < 0) {
1168                 neg = !neg;
1169                 x = -x;
1170         }
1171
1172         /* ??? avoid underflow ??? */
1173
1174         y = x * x;
1175         x += x * y * POLYNOM7(y, r);
1176         return neg ? -x : x;
1177 }
1178
1179 double
1180 sin(x)
1181         double x;
1182 {
1183         return sinus(x, 0);
1184 }
1185
1186 double
1187 cos(x)
1188         double x;
1189 {
1190         if (x < 0) x = -x;
1191         return sinus(x, 1);
1192 }
1193 +sinh.c\0\0\0\0\0\0\0\0\0\0\0\0\ 2\ 2�\ 1\0\0�\ 5/*
1194  * (c) copyright 1989 by the Vrije Universiteit, Amsterdam, The Netherlands.
1195  * See the copyright notice in the ACK home directory, in the file "Copyright".
1196  *
1197  * Author: Ceriel J.H. Jacobs
1198  */
1199
1200 /* $Id: sinh.c,v 1.3 1994/06/24 12:23:17 ceriel Exp $ */
1201
1202 #include <math.h>
1203 #include <errno.h>
1204
1205 extern int      errno;
1206 extern double   exp();
1207
1208 static double
1209 sinh_cosh(x, cosh_flag)
1210         double  x;
1211 {
1212         /*      Algorithm and coefficients from:
1213                         "Software manual for the elementary functions"
1214                         by W.J. Cody and W. Waite, Prentice-Hall, 1980
1215         */
1216
1217         static double p[] = {
1218                 -0.35181283430177117881e+6,
1219                 -0.11563521196851768270e+5,
1220                 -0.16375798202630751372e+3,
1221                 -0.78966127417357099479e+0
1222         };
1223         static double q[] = {
1224                 -0.21108770058106271242e+7,
1225                  0.36162723109421836460e+5,
1226                 -0.27773523119650701167e+3,
1227                  1.0
1228         };
1229         int     negative = x < 0;
1230         double  y = negative ? -x : x;
1231
1232         if (! cosh_flag && y <= 1.0) {
1233                 /* ??? check for underflow ??? */
1234                 y = y * y;
1235                 return x + x * y * POLYNOM3(y, p)/POLYNOM3(y,q);
1236         }
1237
1238         if (y >= M_LN_MAX_D) {
1239                 /* exp(y) would cause overflow */
1240 #define LNV     0.69316101074218750000e+0
1241 #define VD2M1   0.52820835025874852469e-4
1242                 double  w = y - LNV;
1243
1244                 if (w < M_LN_MAX_D+M_LN2-LNV) {
1245                         x = exp(w);
1246                         x += VD2M1 * x;
1247                 }
1248                 else {
1249                         errno = ERANGE;
1250                         x = HUGE;
1251                 }
1252         }
1253         else {
1254                 double  z = exp(y);
1255
1256                 x = 0.5 * (z + (cosh_flag ? 1.0 : -1.0)/z);
1257         }
1258         return negative ? -x : x;
1259 }
1260
1261 double
1262 sinh(x)
1263         double x;
1264 {
1265         return sinh_cosh(x, 0);
1266 }
1267
1268 double
1269 cosh(x)
1270         double x;
1271 {
1272         if (x < 0) x = -x;
1273         return sinh_cosh(x, 1);
1274 }
1275
1276 sqrt.c\0\0\0\0\0\0\0\0\0\0\0\0\ 2\ 2�\ 1\0\0\0\ 3/*
1277  * (c) copyright 1988 by the Vrije Universiteit, Amsterdam, The Netherlands.
1278  * See the copyright notice in the ACK home directory, in the file "Copyright".
1279  *
1280  * Author: Ceriel J.H. Jacobs
1281  */
1282
1283 /* $Id: sqrt.c,v 1.2 1994/06/24 12:23:20 ceriel Exp $ */
1284
1285 #include <math.h>
1286 #include <errno.h>
1287
1288 extern int errno;
1289
1290 #define NITER   5
1291
1292 double
1293 sqrt(x)
1294         double x;
1295 {
1296         extern double frexp(), ldexp();
1297         int exponent;
1298         double val;
1299
1300         if (x <= 0) {
1301                 if (x < 0) errno = EDOM;
1302                 return 0;
1303         }
1304
1305         val = frexp(x, &exponent);
1306         if (exponent & 1) {
1307                 exponent--;
1308                 val *= 2;
1309         }
1310         val = ldexp(val + 1.0, exponent/2 - 1);
1311         /* was: val = (val + 1.0)/2.0; val = ldexp(val, exponent/2); */
1312         for (exponent = NITER - 1; exponent >= 0; exponent--) {
1313                 val = (val + x / val) / 2.0;
1314         }
1315         return val;
1316 }
1317 tan.c\0\0\0\0\0\0\0\0\0\0\0\0\0\ 2\ 2�\ 1\0\0�\ 5/*
1318  * (c) copyright 1988 by the Vrije Universiteit, Amsterdam, The Netherlands.
1319  * See the copyright notice in the ACK home directory, in the file "Copyright".
1320  *
1321  * Author: Ceriel J.H. Jacobs
1322  */
1323
1324 /* $Id: tan.c,v 1.4 1994/06/24 12:23:23 ceriel Exp $ */
1325
1326 #include <math.h>
1327 #include <errno.h>
1328
1329 extern int      errno;
1330 extern double   modf();
1331
1332 double
1333 tan(x)
1334         double x;
1335 {
1336         /*      Algorithm and coefficients from:
1337                         "Software manual for the elementary functions"
1338                         by W.J. Cody and W. Waite, Prentice-Hall, 1980
1339         */
1340
1341         int negative = x < 0;
1342         int invert = 0;
1343         double  y;
1344         static double   p[] = {
1345                  1.0,
1346                 -0.13338350006421960681e+0,
1347                  0.34248878235890589960e-2,
1348                 -0.17861707342254426711e-4
1349         };
1350         static double   q[] = {
1351                  1.0,
1352                 -0.46671683339755294240e+0,
1353                  0.25663832289440112864e-1,
1354                 -0.31181531907010027307e-3,
1355                  0.49819433993786512270e-6
1356         };
1357
1358         if (negative) x = -x;
1359
1360         /* ??? avoid loss of significance, error if x is too large ??? */
1361
1362         y = x * M_2_PI + 0.5;
1363
1364         /*      Use extended precision to calculate reduced argument.
1365                 Here we used 12 bits of the mantissa for a1.
1366                 Also split x in integer part x1 and fraction part x2.
1367         */
1368 #define A1 1.57080078125
1369 #define A2 -4.454455103380768678308e-6
1370         {
1371                 double x1, x2;
1372
1373                 modf(y, &y);
1374                 if (modf(0.5*y, &x1)) invert = 1;
1375                 x2 = modf(x, &x1);
1376                 x = x1 - y * A1;
1377                 x += x2;
1378                 x -= y * A2;
1379 #undef A1
1380 #undef A2
1381         }
1382
1383         /* ??? avoid underflow ??? */
1384         y = x * x;
1385         x += x * y * POLYNOM2(y, p+1);
1386         y = POLYNOM4(y, q);
1387         if (negative) x = -x;
1388         return invert ? -y/x : x/y;
1389 }
1390 tanh.c\0\0\0\0\0\0\0\0\0\0\0\0\ 2\ 2�\ 1\0\0D\ 4/*
1391  * (c) copyright 1989 by the Vrije Universiteit, Amsterdam, The Netherlands.
1392  * See the copyright notice in the ACK home directory, in the file "Copyright".
1393  *
1394  * Author: Ceriel J.H. Jacobs
1395  */
1396
1397 /* $Id: tanh.c,v 1.3 1994/06/24 12:23:27 ceriel Exp $ */
1398
1399 #include <math.h>
1400 #include <errno.h>
1401
1402 extern int      errno;
1403 extern double   exp();
1404
1405 double
1406 tanh(x)
1407         double  x;
1408 {
1409         /*      Algorithm and coefficients from:
1410                         "Software manual for the elementary functions"
1411                         by W.J. Cody and W. Waite, Prentice-Hall, 1980
1412         */
1413
1414         static double p[] = {
1415                 -0.16134119023996228053e+4,
1416                 -0.99225929672236083313e+2,
1417                 -0.96437492777225469787e+0
1418         };
1419         static double q[] = {
1420                  0.48402357071988688686e+4,
1421                  0.22337720718962312926e+4,
1422                  0.11274474380534949335e+3,
1423                  1.0
1424         };
1425         int     negative = x < 0;
1426
1427         if (negative) x = -x;
1428
1429         if (x >= 0.5*M_LN_MAX_D) {
1430                 x = 1.0;
1431         }
1432 #define LN3D2   0.54930614433405484570e+0       /* ln(3)/2 */
1433         else if (x > LN3D2) {
1434                 x = 0.5 - 1.0/(exp(x+x)+1.0);
1435                 x += x;
1436         }
1437         else {
1438                 /* ??? avoid underflow ??? */
1439                 double g = x*x;
1440                 x += x * g * POLYNOM2(g, p)/POLYNOM3(g, q);
1441         }
1442         return negative ? -x : x;
1443 }
1444 exp.c\0\0\0\0\0\0\0\0\0\0\0\0\0\ 2\ 2�\ 1\0\0|\ 5/*
1445  * (c) copyright 1989 by the Vrije Universiteit, Amsterdam, The Netherlands.
1446  * See the copyright notice in the ACK home directory, in the file "Copyright".
1447  *
1448  * Author: Ceriel J.H. Jacobs
1449  */
1450
1451 /* $Id: exp.c,v 1.5 1994/06/24 12:22:42 ceriel Exp $ */
1452
1453 #include <math.h>
1454 #include <errno.h>
1455
1456 extern int      errno;
1457 extern double   ldexp();
1458
1459 double
1460 exp(x)
1461         double  x;
1462 {
1463         /*      Algorithm and coefficients from:
1464                         "Software manual for the elementary functions"
1465                         by W.J. Cody and W. Waite, Prentice-Hall, 1980
1466         */
1467
1468         static double p[] = {
1469                 0.25000000000000000000e+0,
1470                 0.75753180159422776666e-2,
1471                 0.31555192765684646356e-4
1472         };
1473
1474         static double q[] = {
1475                 0.50000000000000000000e+0,
1476                 0.56817302698551221787e-1,
1477                 0.63121894374398503557e-3,
1478                 0.75104028399870046114e-6
1479         };
1480         double  xn, g;
1481         int     n;
1482         int     negative = x < 0;
1483
1484         if (x <= M_LN_MIN_D) {
1485                 if (x < M_LN_MIN_D) errno = ERANGE;
1486                 return M_MIN_D;
1487         }
1488         if (x >= M_LN_MAX_D) {
1489                 if (x > M_LN_MAX_D) errno = ERANGE;
1490                 return M_MAX_D;
1491         }
1492         if (negative) x = -x;
1493
1494         /* ??? avoid underflow ??? */
1495
1496         n = x * M_LOG2E + 0.5;  /* 1/ln(2) = log2(e), 0.5 added for rounding */
1497         xn = n;
1498         {
1499                 double  x1 = (long) x;
1500                 double  x2 = x - x1;
1501
1502                 g = ((x1-xn*0.693359375)+x2) - xn*(-2.1219444005469058277e-4);
1503         }
1504         if (negative) {
1505                 g = -g;
1506                 n = -n;
1507         }
1508         xn = g * g;
1509         x = g * POLYNOM2(xn, p);
1510         n += 1;
1511         return (ldexp(0.5 + x/(POLYNOM3(xn, q) - x), n));
1512 }
1513 floor.c\0\0\0\0\0\0\0\0\0\0\0\ 2\ 2�\ 1\0\0�\ 1/*
1514  * (c) copyright 1988 by the Vrije Universiteit, Amsterdam, The Netherlands.
1515  * See the copyright notice in the ACK home directory, in the file "Copyright".
1516  *
1517  * Author: Ceriel J.H. Jacobs
1518  */
1519
1520 /* $Id: floor.c,v 1.2 1994/06/24 12:22:48 ceriel Exp $ */
1521
1522 double
1523 floor(x)
1524         double x;
1525 {
1526         extern double modf();
1527         double val;
1528
1529         return modf(x, &val) < 0 ? val - 1.0 : val ;
1530         /*      this also works if modf always returns a positive
1531                 fractional part
1532         */
1533 }
1534 "