Integer square root algorithm

sqrt.c

/*
Just like integer division, these functions round down to an integer.
Running this program should produce the following output:

sqrt(3735928559) == 61122
61122^2 == 3735898884
sqrt(244837814094590) == 15647294
15647294^2 == 244837809522436

This algorithm comes from Jack W. Crenshaw's 1998 article in Embedded:
http://www.embedded.com/electronics-blogs/programmer-s-toolbox/4219659/Integer-Square-Roots
*/

#include <stdio.h>
#include <stdint.h>

uint16_t gapsqrt32(uint32_t a);
uint32_t gapsqrt64(uint64_t a);

int main(void) {
	uint32_t small = 0xdeadbeef;
	uint32_t sqrt_small = gapsqrt32(small);
	uint32_t check_small = sqrt_small * sqrt_small;

	printf("sqrt(%u) == %u\n", small, sqrt_small);
	printf("%u^2 == %u\n", sqrt_small, check_small);

	uint64_t large = 0xdeadbeefcafe;
	uint64_t sqrt_large = gapsqrt64(large);
	uint64_t check_large = sqrt_large * sqrt_large;

	printf("sqrt(%lu) == %lu\n", large, sqrt_large);
	printf("%lu^2 == %lu\n", sqrt_large, check_large);

	return 0;
}

uint16_t gapsqrt32(uint32_t a) {
        uint32_t rem = 0, root = 0;

        for (int i = 32 / 2; i > 0; i--) {
                root <<= 1;
                rem = (rem << 2) | (a >> (32 - 2));
                a <<= 2;
                if (root < rem) {
                        rem -= root | 1;
                        root += 2;
                }
        }
        return root >> 1;
}


uint32_t gapsqrt64(uint64_t a) {
        uint64_t rem = 0, root = 0;

        for (int i = 64 / 2; i > 0; i--) {
                root <<= 1;
                rem = (rem << 2) | (a >> (64 - 2));
                a <<= 2;
                if (root < rem) {
                        rem -= root | 1;
                        root += 2;
                }
        }
        return root >> 1;
}

`// integer square root of a positive number
// http://ww1.microchip.com/downloads/en/AppNotes/91040a.pdf
uint32_t sqrt(uint64_t num) {
uint32_t res = 0; // result
uint32_t add = 1;
add <<= 31; // 'additional' bit is in position 31
uint32_t temp; // 'A'
uint64_t quad; // 'A^2'

while ( add > 0 ) { // 32x test and shift right
temp = res + add;
quad = temp;
quad *= temp;
if ( num >= quad ) {
res = temp;
}
add >>= 1; // shift right the 'additional' bit
}
return res;
}
`