You can not select more than 25 topics
Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.
52 lines
675 B
52 lines
675 B
#include <soft.h>
|
|
|
|
/*
|
|
sqrt returns the square root of its floating
|
|
point argument. Newton's method.
|
|
|
|
calls frexp
|
|
need error handling
|
|
*/
|
|
|
|
double sqrt(double arg)
|
|
{
|
|
double x, temp;
|
|
int exp;
|
|
int i;
|
|
|
|
if(arg <= 0.) {
|
|
errno = EDOM;
|
|
return(0.);
|
|
}
|
|
x = frexp(arg,&exp);
|
|
while(x < 0.5) {
|
|
x *= 2;
|
|
exp--;
|
|
}
|
|
/*
|
|
* NOTE
|
|
* this wont work on 1's comp
|
|
*/
|
|
if(exp & 1) {
|
|
x *= 2;
|
|
exp--;
|
|
}
|
|
temp = 0.5*(1.0+x);
|
|
|
|
while(exp > 60) {
|
|
temp *= (1L<<30);
|
|
exp -= 60;
|
|
}
|
|
while(exp < -60) {
|
|
temp /= (1L<<30);
|
|
exp += 60;
|
|
}
|
|
if(exp >= 0)
|
|
temp *= 1L << (exp/2);
|
|
else
|
|
temp /= 1L << (-exp/2);
|
|
for(i=0; i<=4; i++)
|
|
temp = 0.5*(temp + arg/temp);
|
|
return(temp);
|
|
}
|