72 lines
1.1 KiB
C
72 lines
1.1 KiB
C
/*
|
|
* (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$ */
|
|
#define __NO_DEFS
|
|
#include <math.h>
|
|
#include <pc_err.h>
|
|
extern _trp();
|
|
|
|
#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) _trp(ESQT);
|
|
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;
|
|
}
|