sqt.c 1.1 KB

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