1234567891011121314151617181920212223242526272829303132333435363738394041424344454647484950515253545556575859606162636465666768697071 |
- /*
- * (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
- */
- /* $Id$ */
- #define __NO_DEFS
- #include <math.h>
- #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;
- }
|