Subversion Repositories planix.SVN

Rev

Details | Last modification | View Log | RSS feed

Rev Author Line No. Line
2 - 1 
#include <u.h>
2 
#include <libc.h>
3 
 
4 
static	long	sqtab[64] =
5 
{
6 
	0x6cdb2, 0x726d4, 0x77ea3, 0x7d52f, 0x82a85, 0x87eb1, 0x8d1c0, 0x923bd,
7 
	0x974b2, 0x9c4a8, 0xa13a9, 0xa61be, 0xaaeee, 0xafb41, 0xb46bf, 0xb916e,
8 
	0xbdb55, 0xc247a, 0xc6ce3, 0xcb495, 0xcfb95, 0xd41ea, 0xd8796, 0xdcca0,
9 
	0xe110c, 0xe54dd, 0xe9818, 0xedac0, 0xf1cd9, 0xf5e67, 0xf9f6e, 0xfdfef,
10 
	0x01fe0, 0x05ee6, 0x09cfd, 0x0da30, 0x11687, 0x1520c, 0x18cc8, 0x1c6c1,
11 
	0x20000, 0x2388a, 0x27068, 0x2a79e, 0x2de32, 0x3142b, 0x3498c, 0x37e5b,
12 
	0x3b29d, 0x3e655, 0x41989, 0x44c3b, 0x47e70, 0x4b02b, 0x4e16f, 0x51241,
13 
	0x542a2, 0x57296, 0x5a220, 0x5d142, 0x60000, 0x62e5a, 0x65c55, 0x689f2,
14 
};
15 
 
16 
double
17 
sqrt(double arg)
18 
{
19 
	int e, ms;
20 
	double a, t;
21 
	union
22 
	{
23 
		double	d;
24 
		struct
25 
		{
26 
			long	ms;
27 
			long	ls;
28 
		};
29 
	} u;
30 
 
31 
	u.d = arg;
32 
	ms = u.ms;
33 
 
34 
	/*
35 
	 * sign extend the mantissa with
36 
	 * exponent. result should be > 0 for
37 
	 * normal case.
38 
	 */
39 
	e = ms >> 20;
40 
	if(e <= 0) {
41 
		if(e == 0)
42 
			return 0;
43 
		return NaN();
44 
	}
45 
 
46 
	/*
47 
	 * pick up arg/4 by adjusting exponent
48 
	 */
49 
	u.ms = ms - (2 << 20);
50 
	a = u.d;
51 
 
52 
	/*
53 
	 * use 5 bits of mantissa and 1 bit
54 
	 * of exponent to form table index.
55 
	 * insert exponent/2 - 1.
56 
	 */
57 
	e = (((e - 1023) >> 1) + 1022) << 20;
58 
	u.ms = *(long*)((char*)sqtab + ((ms >> 13) & 0xfc)) | e;
59 
	u.ls = 0;
60 
 
61 
	/*
62 
	 * three laps of newton
63 
	 */
64 
	e = 1 << 20;
65 
	t = u.d;
66 
	u.d = t + a/t;
67 
	u.ms -= e;		/* u.d /= 2; */
68 
	t = u.d;
69 
	u.d = t + a/t;
70 
	u.ms -= e;		/* u.d /= 2; */
71 
	t = u.d;
72 
 
73 
	return t + a/t;
74 
}
75 
 
76 
/*
77 
 * this is the program that generated the table.
78 
 * it calls sqrt by some other means.
79 
 *
80 
 * void
81 
 * main(void)
82 
 * {
83 
 * 	int i;
84 
 * 	union	U
85 
 * 	{
86 
 * 		double	d;
87 
 * 		struct
88 
 * 		{
89 
 * 			long	ms;
90 
 * 			long	ls;
91 
 * 		};
92 
 * 	} u;
93 
 *
94 
 * 	for(i=0; i<64; i++) {
95 
 * 		u.ms = (i<<15) | 0x3fe04000;
96 
 * 		u.ls = 0;
97 
 * 		u.d = sqrt(u.d);
98 
 * 		print(" 0x%.5lux,", u.ms & 0xfffff);
99 
 * 	}
100 
 * 	print("\n");
101 
 * 	exits(0);
102 
 * }
103 
 */