2

素数の計算のために平方根を求める関数を作成しています。ただし浮動小数の平方根関数は精度の問題で使用しませんでした。関数は以下サイトを参考に作成しました。

64bit数の素数判定

作成したコードは以下の通りです。

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

ちなみに、二分法でもコードは書いているのですが、できるだけ簡潔で速いコードを書きたいということでニュートン・ラフソン法を使っています。もしもっといいアルゴリズムがあったら教えていただければと思います。

1 件の回答 1

3

sqrt_iから返すべき値は、yではなくy_oldです。

参考サイトのソースは、以下のようになっています。

// ニュートン法による整数平方根
fn sqrt_newton(x: u64) -> u64 {
    if x <= 1 { return x; }
    let k = 32 - ((x - 1).leading_zeros() >> 1);
    let mut s = (1 as u64) << k; // s = 2**k
    let mut t = (s + (x >> k)) >> 1; // t = (s + x/s)/2
    // whileループ回数=除算回数は u64:最大6回 u32:最大5回 を想定
    // s > floor(sqrt(x)) -> floor(sqrt(x)) <= t < s
    // s == floor(sqrt(x)) -> s == floor(sqrt(x)) <= t <= floor(sqrt(x)) + 1
    while t < s { s = t; t = (s + (x / s)) >> 1; }
    s
}

参考サイトのtyに、sy_oldに対応します。参考サイトの戻り値はtではなくsとなっています。

参考サイトの説明に「(2) s_n ​= ⌊√a​⌋ のとき、 ⌊√a⌋ ≤ s_(n+1) ≤ ⌊√a⌋+1 である。」とある通り、s_(n+1) は ⌊√a⌋+1 に等しくなる(=正しい値より1だけ大きい値となる)可能性があります。

この質問に回答するには、ログインする必要があります。

求めていた回答ではありませんか? のタグが付いた他の質問を参照する。