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