sqrt.c 722 B

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