C言語でKaratsuba開平を利用して整数の平方根と余りを求める

math.hのsqrt(f)を使わずに整数の平方根を求めてみた。(intが4byteであるマシンでのみ動作)
sqrtmod0はアルゴリズム上0x40000000以上の値の平方根を計算でようなので、0x40000000以上の値はKaratsuba開平を利用して求めている。

↓は乱数で平方根と余りを求めてsqrtで求めたものとassertをし続けるプログラムです。

#include <stdio.h>
#include <stdlib.h>
#include <assert.h>
#include <math.h>

#define BASE (sizeof (unsigned int) * 8 / 4)

int sqrtmod(unsigned int n, unsigned int *pmod);
int sqrtmod_karatsuba(unsigned int n, unsigned int *pmod);
int sqrtmod0(unsigned int n, unsigned int *pmod);

int sqrtmod(unsigned int n, unsigned int *pmod) {
	return sqrtmod_karatsuba(n, pmod);
}

int sqrtmod_karatsuba(unsigned int n, unsigned int *pmod) {
	int s, s1, s0, r, r1, r0;

	if ((0xc0000000 & n) == 0) {
		return sqrtmod0(n, pmod);
	}

	s1 = sqrtmod0(n >> (BASE * 2), &r1);

	r0 = (r1 << BASE) | ((n >> BASE) & 0x000000ff);
	s0 = r0 / (s1 << 1);
	r0 -= s0 * (s1 << 1);

	s = (s1 << BASE) + s0;
	r = ((r0 << BASE) | (n & 0x000000ff)) - s0 * s0;
	if (r < 0) {
		r += (s << 1) - 1;
		s--;
	}
	if (pmod) *pmod = r;
	assert(s * s + r == n);

	return s;
}

int sqrtmod0(unsigned int n, unsigned int *pmod) {
	int p = 0, q = 1, r = n, h = 0;

	while (q <= n) {
		q = q << 2;
	}

	while (1 != q) {
		q = q >> 2;
		h = p + q;
		p = p >> 1;
		if (r >= h) {
			p = p + q;
			r = r - h;
		}
	}
	if (pmod) *pmod = r;
	assert(p * p + r == n);

	return p;
}

int main(int argc, char **argv) {
	unsigned int n, s1, r1, s2, r2;

	srand((unsigned)time(NULL));
	while (1) {
		n = rand();
		s1 = sqrtmod(n, &r1);
		s2 = sqrt(n);
		r2 = n - s2 * s2;
		assert(s1 == s2 && r1 == r2);
	}

	return 0;
}

とりあえずassertエラーが出ないので、実装は出来ているらしい。。

参考：
2 Karatsuba$B7ONs$N%"%k%4%j%:%`(B
GNU MP 6.1.2: Square Root Algorithm