pow.c 1.5 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100
  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. #include <math.h>
  9. #include <float.h>
  10. #include <errno.h>
  11. #include <limits.h>
  12. double
  13. pow(double x, double y)
  14. {
  15. double y_intpart, y_fractpart, fp;
  16. int negexp, negx;
  17. int ex, newexp;
  18. unsigned long yi;
  19. if (x == 1.0) return x;
  20. if (x == 0 && y <= 0) {
  21. errno = EDOM;
  22. return 0;
  23. }
  24. if (y == 0) return 1.0;
  25. if (y < 0) {
  26. y = -y;
  27. negexp = 1;
  28. }
  29. else negexp = 0;
  30. y_fractpart = modf(y, &y_intpart);
  31. if (y_fractpart != 0) {
  32. if (x < 0) {
  33. errno = EDOM;
  34. return 0;
  35. }
  36. }
  37. negx = 0;
  38. if (x < 0) {
  39. x = -x;
  40. negx = 1;
  41. }
  42. if (y_intpart > ULONG_MAX) {
  43. if (negx && modf(y_intpart/2.0, &y_fractpart) == 0) {
  44. negx = 0;
  45. }
  46. x = log(x);
  47. /* Beware of overflow in the multiplication */
  48. if (x > 1.0 && y > DBL_MAX/x) {
  49. errno = ERANGE;
  50. return HUGE_VAL;
  51. }
  52. if (negexp) y = -y;
  53. if (negx) return -exp(x*y);
  54. return exp(x * y);
  55. }
  56. if (y_fractpart != 0) {
  57. fp = exp(y_fractpart * log(x));
  58. }
  59. else fp = 1.0;
  60. yi = y_intpart;
  61. if (! (yi & 1)) negx = 0;
  62. x = frexp(x, &ex);
  63. newexp = 0;
  64. for (;;) {
  65. if (yi & 1) {
  66. fp *= x;
  67. newexp += ex;
  68. }
  69. yi >>= 1;
  70. if (yi == 0) break;
  71. x *= x;
  72. ex <<= 1;
  73. if (x < 0.5) {
  74. x += x;
  75. ex -= 1;
  76. }
  77. }
  78. if (negexp) {
  79. fp = 1.0/fp;
  80. newexp = -newexp;
  81. }
  82. if (negx) {
  83. return -ldexp(fp, newexp);
  84. }
  85. return ldexp(fp, newexp);
  86. }