ref: 97f4b9b97c470501d25ac82cad771bdff75e9b59
dir: /sqrt.c/
#include <u.h>
#include <libc.h>
int iters;
/*
* Heron's method to compute the √
*
* iteratively do
* x1 = ½(x0 + n/x0)
* since
* lim M→∞ (xM) = √n
*/
//double
//√(double n)
//{
// int i;
// double x;
//
// x = 2;
// for(i = 0; i < iters; i++)
// x = 0.5*(x + n/x);
// return x;
//}
double
√(double n)
{
double x0, x;
if(n == 0)
return 0;
x0 = -1;
x = n > 1? n/2: 1;
while(x0 != x){
x0 = x;
x = 0.5*(x0 + n/x0);
iters++;
}
return x;
}
void
usage(void)
{
fprint(2, "usage: %s number [prec]\n", argv0);
exits("usage");
}
void
main(int argc, char *argv[])
{
int prec;
double n;
prec = 10;
ARGBEGIN{
default: usage();
}ARGEND
if(argc < 1)
usage();
n = strtod(argv[0], nil);
if(n < 0)
sysfatal("too complex");
if(argc > 2)
prec = strtoul(argv[1], nil, 10);
print("√%g = %.*f (took %d iterations)\n", n, prec, √(n), iters);
exits(nil);
}