sqt.c 1.1 KB

1234567891011121314151617181920212223242526272829303132333435363738394041424344454647484950515253545556575859606162636465666768697071
  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. #define NITER 5
  11. static double
  12. ldexp(fl,exp)
  13. double fl;
  14. int exp;
  15. {
  16. extern double _fef();
  17. int sign = 1;
  18. int currexp;
  19. if (fl<0) {
  20. fl = -fl;
  21. sign = -1;
  22. }
  23. fl = _fef(fl,&currexp);
  24. exp += currexp;
  25. if (exp > 0) {
  26. while (exp>30) {
  27. fl *= (double) (1L << 30);
  28. exp -= 30;
  29. }
  30. fl *= (double) (1L << exp);
  31. }
  32. else {
  33. while (exp<-30) {
  34. fl /= (double) (1L << 30);
  35. exp += 30;
  36. }
  37. fl /= (double) (1L << -exp);
  38. }
  39. return sign * fl;
  40. }
  41. double
  42. _sqt(x)
  43. double x;
  44. {
  45. extern double _fef();
  46. int exponent;
  47. double val;
  48. if (x <= 0) {
  49. if (x < 0) error(3);
  50. return 0;
  51. }
  52. val = _fef(x, &exponent);
  53. if (exponent & 1) {
  54. exponent--;
  55. val *= 2;
  56. }
  57. val = ldexp(val + 1.0, exponent/2 - 1);
  58. /* was: val = (val + 1.0)/2.0; val = ldexp(val, exponent/2); */
  59. for (exponent = NITER - 1; exponent >= 0; exponent--) {
  60. val = (val + x / val) / 2.0;
  61. }
  62. return val;
  63. }