blob: 27902e10e2f0af14ac01e0b576e7e378fa666f85 [file] [log] [blame]
 /* +++Date last modified: 05-Jul-1997 */ #include #include "snipmath.h" #define BITSPERLONG 32 #define TOP2BITS(x) ((x & (3 << (BITSPERLONG-2))) >> (BITSPERLONG-2)) /* usqrt: ENTRY x: unsigned int EXIT returns floor(sqrt(x) * pow(2, BITSPERLONG/2)) Since the square root never uses more than half the bits of the input, we use the other half of the bits to contain extra bits of precision after the binary point. EXAMPLE suppose BITSPERLONG = 32 then usqrt(144) = 786432 = 12 * 65536 usqrt(32) = 370727 = 5.66 * 65536 NOTES (1) change BITSPERLONG to BITSPERLONG/2 if you do not want the answer scaled. Indeed, if you want n bits of precision after the binary point, use BITSPERLONG/2+n. The code assumes that BITSPERLONG is even. (2) This is really better off being written in assembly. The line marked below is really a "arithmetic shift left" on the double-long value with r in the upper half and x in the lower half. This operation is typically expressible in only one or two assembly instructions. (3) Unrolling this loop is probably not a bad idea. ALGORITHM The calculations are the base-two analogue of the square root algorithm we all learned in grammar school. Since we're in base 2, there is only one nontrivial trial multiplier. Notice that absolutely no multiplications or divisions are performed. This means it'll be fast on a wide range of processors. */ void usqrt(unsigned int x, struct int_sqrt *q) { unsigned int a = 0; /* accumulator */ unsigned int r = 0; /* remainder */ unsigned int e = 0; /* trial product */ int i; for (i = 0; i < BITSPERLONG; i++) /* NOTE 1 */ { r = (r << 2) + TOP2BITS(x); x <<= 2; /* NOTE 2 */ a <<= 1; e = (a << 1) + 1; if (r >= e) { r -= e; a++; } } memcpy(q, &a, sizeof(int)); } #ifdef TEST #include #include main(void) { int i; unsigned long l = 0x3fed0169L; struct int_sqrt q; for (i = 0; i < 101; ++i) { usqrt(i, &q); printf("sqrt(%3d) = %2d, remainder = %2d\n", i, q.sqrt, q.frac); } usqrt(l, &q); printf("\nsqrt(%lX) = %X, remainder = %X\n", l, q.sqrt, q.frac); return 0; } #endif /* TEST */