atn.c 1.2 KB

1234567891011121314151617181920212223242526272829303132333435363738394041424344454647484950515253545556575859606162636465666768
  1. /*
  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. /* $Id$ */
  8. #define __NO_DEFS
  9. #include <math.h>
  10. double
  11. _atn(x)
  12. double x;
  13. {
  14. /* Algorithm and coefficients from:
  15. "Software manual for the elementary functions"
  16. by W.J. Cody and W. Waite, Prentice-Hall, 1980
  17. */
  18. static double p[] = {
  19. -0.13688768894191926929e+2,
  20. -0.20505855195861651981e+2,
  21. -0.84946240351320683534e+1,
  22. -0.83758299368150059274e+0
  23. };
  24. static double q[] = {
  25. 0.41066306682575781263e+2,
  26. 0.86157349597130242515e+2,
  27. 0.59578436142597344465e+2,
  28. 0.15024001160028576121e+2,
  29. 1.0
  30. };
  31. static double a[] = {
  32. 0.0,
  33. 0.52359877559829887307710723554658381, /* pi/6 */
  34. M_PI_2,
  35. 1.04719755119659774615421446109316763 /* pi/3 */
  36. };
  37. int neg = x < 0;
  38. int n;
  39. double g;
  40. if (neg) {
  41. x = -x;
  42. }
  43. if (x > 1.0) {
  44. x = 1.0/x;
  45. n = 2;
  46. }
  47. else n = 0;
  48. if (x > 0.26794919243112270647) { /* 2-sqtr(3) */
  49. n = n + 1;
  50. x = (((0.73205080756887729353*x-0.5)-0.5)+x)/
  51. (1.73205080756887729353+x);
  52. }
  53. /* ??? avoid underflow ??? */
  54. g = x * x;
  55. x += x * g * POLYNOM3(g, p) / POLYNOM4(g, q);
  56. if (n > 1) x = -x;
  57. x += a[n];
  58. return neg ? -x : x;
  59. }