/* * (c) copyright 1988 by the Vrije Universiteit, Amsterdam, The Netherlands. * See the copyright notice in the ACK home directory, in the file "Copyright". * * Author: Ceriel J.H. Jacobs */ /* $Header$ */ #include #define NITER 5 static double ldexp(fl,exp) double fl; int exp; { extern double _fef(); int sign = 1; int currexp; if (fl<0) { fl = -fl; sign = -1; } fl = _fef(fl,&currexp); exp += currexp; if (exp > 0) { while (exp>30) { fl *= (double) (1L << 30); exp -= 30; } fl *= (double) (1L << exp); } else { while (exp<-30) { fl /= (double) (1L << 30); exp += 30; } fl /= (double) (1L << -exp); } return sign * fl; } double _sqt(x) double x; { extern double _fef(); int exponent; double val; if (x <= 0) { if (x < 0) error(3); return 0; } val = _fef(x, &exponent); if (exponent & 1) { exponent--; val *= 2; } val = ldexp(val + 1.0, exponent/2 - 1); /* was: val = (val + 1.0)/2.0; val = ldexp(val, exponent/2); */ for (exponent = NITER - 1; exponent >= 0; exponent--) { val = (val + x / val) / 2.0; } return val; }