素数の計算のために平方根を求める関数を作成しています。ただし浮動小数の平方根関数は精度の問題で使用しませんでした。関数は以下サイトを参考に作成しました。
作成したコードは以下の通りです。
/*
* def.h
*/
#ifndef DEF_H_40D70383_4030_8346_DF89_F0F6D9210451
#define DEF_H_40D70383_4030_8346_DF89_F0F6D9210451
typedef long long int INT;
extern INT sqrt_i(INT x);
#endif
/*
* sqrt_i_c.c
*/
#include <stdio.h>
#include "def.h"
#include <stdlib.h>
INT sqrt_i(INT x) {
INT y, y_old;
if (x >= (1LL << 62) || x <= 0) {
fprintf(stderr, "ERROR: OUT OF RANGE, in sqrt_i()\n");
exit(1);
}
{
INT tmp_x = x >> 2;
INT min = 1LL;
for (; tmp_x != 0LL; tmp_x >>= 2) min <<= 1;
y = (min << 1) - 1;
}
do {
y_old = y;
y = (y + x / y) >> 1;
} while (y < y_old);
if (y * y > x) y--; /* この行を追加しないと単体テストが通らない */
return y;
}
/*
* test_sqrt_i.c
*/
#include <stdio.h>
#include "def.h"
void check(INT x, INT sqrt_x) {
if (sqrt_i(x) != sqrt_x) {
fprintf(stderr, "ASSERTION FAILED: x=%lld, sqrt_x=%lld, sqrt_i(x)=%lld\n", x, sqrt_x, sqrt_i(x));
}
}
int main() {
check(1LL, 1LL);
check(2LL, 1LL);
check(3LL, 1LL);
check(4LL, 2LL);
check(5LL, 2LL);
check(6LL, 2LL);
check(7LL, 2LL);
check(8LL, 2LL);
check(9LL, 3LL);
check(10LL, 3LL);
check(11LL, 3LL);
check(12LL, 3LL);
check(13LL, 3LL);
check(14LL, 3LL);
check(15LL, 3LL);
check(16LL, 4LL);
check(123456LL * 123456LL, 123456LL);
check(123456LL * 123456LL - 1, 123455LL);
check(123456LL * 123456LL + 1, 123456LL);
return 0;
}
平方根の計算に、ニュートン・ラフソン法を用いているのですが、if (y * y > x) y--;
というコードを追加しないと単体テストが失敗してしまいます。失敗するのは、平方数の1個前の数の場合です。
ASSERTION FAILED: x=3, sqrt_x=1, sqrt_i(x)=2
ASSERTION FAILED: x=8, sqrt_x=2, sqrt_i(x)=3
ASSERTION FAILED: x=15, sqrt_x=3, sqrt_i(x)=4
ASSERTION FAILED: x=15241383935, sqrt_x=123455, sqrt_i(x)=123456
ちなみに、二分法でもコードは書いているのですが、できるだけ簡潔で速いコードを書きたいということでニュートン・ラフソン法を使っています。もしもっといいアルゴリズムがあったら教えていただければと思います。