sqrt.c 782 B

12345678910111213141516171819202122232425262728293031323334353637383940414243
  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. #define NITER 5
  12. double
  13. sqrt(double x)
  14. {
  15. int exponent;
  16. double val;
  17. if (__IsNan(x)) {
  18. errno = EDOM;
  19. return x;
  20. }
  21. if (x <= 0) {
  22. if (x < 0) errno = EDOM;
  23. return 0;
  24. }
  25. if (x > DBL_MAX) return x; /* for infinity */
  26. val = frexp(x, &exponent);
  27. if (exponent & 1) {
  28. exponent--;
  29. val *= 2;
  30. }
  31. val = ldexp(val + 1.0, exponent/2 - 1);
  32. /* was: val = (val + 1.0)/2.0; val = ldexp(val, exponent/2); */
  33. for (exponent = NITER - 1; exponent >= 0; exponent--) {
  34. val = (val + x / val) / 2.0;
  35. }
  36. return val;
  37. }