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 */ 104