blob: 281e73e801361e71667762e222400bfe87722988 (
plain)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
|
#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; /* initial estimate */
/*
* take advantage of the computer's discreteness
* to get the most accurate result.
*/
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);
}
|
