LoginSignup
2
1

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?

@mikecat_mixc

浮動小数点演算で確実に真の値以下の値が欲しいなら、丸めモードを設定しろ

Posted at

記事
Integer.sqrt(4503599761588224)**2 > 4503599761588224 #Ruby - Qiita
では、指定された整数の平方根以下の最大の整数を返すはずの関数が、入力 4503599761588224 に対しては間違った値 (平方根より大きい値) を返す事例が紹介されている。
この記事によると、これはこの関数が内部で浮動小数点の sqrt を用いているために誤差が出てしまったかららしい。

間違った値が返ることを確認する

まずは、実際の Integer.sqrt の出力を見てみる。

print Integer.sqrt(4503599761588224)

実行すると、以下の出力が得られた。(ruby 3.4.1)

67108865

4503599761588224 = {67108865}^2 - 1

なので、出力するべき値は 67108864 であり、たしかに間違った値が出力された。

浮動小数点演算の結果を確認する

C言語で浮動小数点演算による 4503599761588224 の平方根を求めた結果を確認する。

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

int main(void) {
    double v = sqrt(4503599761588224.0);
    long double vl = sqrtl(4503599761588224.0);
    printf("%.40f %.20a\n", v, v);
    printf("%.40Lf %.20La\n", vl, vl);
    return 0;
}

gcc 13.2.0 で実行すると、以下の出力が得られた。

67108865.0000000000000000000000000000000000000000 0x1.00000040000000000000p+26
67108864.9999999925494194030761718750000000000000 0x8.000001ffffffc0000000p+23

clang 17.0.1 で実行すると、以下の出力が得られた。

67108865.0000000000000000000000000000000000000000 0x1.00000040000000000000p+26
67108864.9999999925494194030761718750000000000000 0x8.000001ffffffc0000000p+23

double での計算結果が、切り上げられて 67108865 になってしまっている。
long double での計算結果は、67108865 未満の値になった。

丸めモードを指定して計算する

x87 における浮動小数点数の計算では、「近い値に丸める」以外にも「0に近い値に丸める」「負の無限大に近い値に丸める」などの丸めモードを選択できたはずである。
これらの丸めモードを選択することで、double で計算しても正しい値が求まるかも知れない。

C言語では、fesetround() 関数を用いることで、丸めモードを設定できる。
設定は、以下の丸めモードを表すマクロを用いて行う。

マクロ 意味
FE_DOWNWARD 負の無限大に近い値に丸める
FE_TONEAREST 近い値に丸める
FE_TOWARDZERO 0に近い値に丸める
FE_UPWARD 正の無限大に近い値に丸める

これらのマクロは、使用可能な環境でのみ定義される。
また、fesetround() は失敗する可能性があり、失敗した場合は非0を返す。

これらを用いて丸めモードを設定した上で、再び 4503599761588224 の平方根を計算してみた。

#include <stdio.h>
#include <math.h>
#include <fenv.h>

int main(void) {
#if defined(FE_DOWNWARD)
    puts("setting rounding mode to FE_DOWNWARD");
    if (fesetround(FE_DOWNWARD)) return 1;
#elif defined(FE_TOWARDZERO)
    puts("setting rounding mode to FE_TOWARDZERO");
    if (fesetround(FE_TOWARDZERO)) return 1;
#else
#error appropriate rounding mode unavailable
#endif
    double v = sqrt(4503599761588224.0);
    long double vl = sqrtl(4503599761588224.0);
    printf("%.40f %.20a\n", v, v);
    printf("%.40Lf %.20La\n", vl, vl);
    return 0;
}

gcc 13.2.0 で実行すると、以下の結果が得られた。

setting rounding mode to FE_DOWNWARD
67108865.0000000000000000000000000000000000000000 0x1.00000040000000000000p+26
67108864.9999999925494194030761718750000000000000 0x8.000001ffffffc0000000p+23

clang 17.0.1 で実行すると、以下の結果が得られた。

setting rounding mode to FE_DOWNWARD
67108864.9999999850988388061523437500000000000000 0x1.0000003ffffff0000000p+26
67108864.9999999925494194030761718750000000000000 0x8.000001ffffffc0000000p+23

double での計算において、clang では期待通り 67108865 未満の結果が得られたが、gcc では変わらずに 67108865 になってしまった。

浮動小数点演算の結果が環境依存なのはどんなときか
では、

libcの実装によっては「sqrt に fesetround が効かない」というような状況があるようです。

という主張がある。

Compiler Explorer の x86-64 gcc 14.2 では、以下のコンパイル結果になった。

gcc でのコンパイル結果 
.LC0:
        .string "setting rounding mode to FE_DOWNWARD"
.LC3:
        .string "%.40f %.20a\n"
.LC4:
        .string "%.40Lf %.20La\n"
main:
        pushq   %rbp
        movq    %rsp, %rbp
        subq    $32, %rsp
        movl    $.LC0, %edi
        call    puts
        movl    $1024, %edi
        call    fesetround
        testl   %eax, %eax
        je      .L2
        movl    $1, %eax
        jmp     .L3
.L2:
        movsd   .LC1(%rip), %xmm0
        movsd   %xmm0, -8(%rbp)
        fldt    .LC2(%rip)
        fstpt   -32(%rbp)
        movsd   -8(%rbp), %xmm0
        movq    -8(%rbp), %rax
        movapd  %xmm0, %xmm1
        movq    %rax, %xmm0
        movl    $.LC3, %edi
        movl    $2, %eax
        call    printf
        pushq   -24(%rbp)
        pushq   -32(%rbp)
        pushq   -24(%rbp)
        pushq   -32(%rbp)
        movl    $.LC4, %edi
        movl    $0, %eax
        call    printf
        addq    $32, %rsp
        movl    $0, %eax
.L3:
        leave
        ret
.LC1:
        .long   67108864
        .long   1099956224
.LC2:
        .long   -1024
        .long   -2147483617
        .long   16409
        .long   0

計算結果を埋め込んで用いており、埋め込む値は fesetround() の有無で変わらないようであった。
さらに、コンパイルオプション -O0 を追加しても、この fesetround() を無視した値を埋め込む挙動は変わらないようだった。

一方、x86-64 clang 20.1.0 では、以下のコンパイル結果になった。

clang でのコンパイル結果 
.LCPI0_0:
        .quad   0x4330000008000000
main:
        pushq   %rbp
        movq    %rsp, %rbp
        subq    $64, %rsp
        movl    $0, -4(%rbp)
        leaq    .L.str(%rip), %rdi
        callq   puts@PLT
        movl    $1024, %edi
        callq   fesetround@PLT
        cmpl    $0, %eax
        je      .LBB0_2
        movl    $1, -4(%rbp)
        jmp     .LBB0_3
.LBB0_2:
        movsd   .LCPI0_0(%rip), %xmm0
        callq   sqrt@PLT
        movsd   %xmm0, -16(%rbp)
        movq    %rsp, %rax
        fldl    .LCPI0_0(%rip)
        fstpt   (%rax)
        callq   sqrtl@PLT
        fstpt   -32(%rbp)
        movsd   -16(%rbp), %xmm1
        leaq    .L.str.1(%rip), %rdi
        movb    $2, %al
        movaps  %xmm1, %xmm0
        callq   printf@PLT
        fldt    -32(%rbp)
        fld     %st(0)
        movq    %rsp, %rax
        fxch    %st(1)
        fstpt   16(%rax)
        fstpt   (%rax)
        leaq    .L.str.2(%rip), %rdi
        xorl    %eax, %eax
        callq   printf@PLT
        movl    $0, -4(%rbp)
.LBB0_3:
        movl    -4(%rbp), %eax
        addq    $64, %rsp
        popq    %rbp
        retq

.L.str:
        .asciz  "setting rounding mode to FE_DOWNWARD"

.L.str.1:
        .asciz  "%.40f %.20a\n"

.L.str.2:
        .asciz  "%.40Lf %.20La\n"

値を埋め込まず、きちんと sqrt 関数を用いて計算を行っていることがわかる。

x86-64 clang 20.1.0 で -O オプションを追加すると、以下のコンパイル結果になった。

clang でのコンパイル結果 (最適化あり) 
.LCPI0_0:
        .quad   0x4330000008000000
.LCPI0_1:
        .quad   0x4190000004000000
main:
        pushq   %rbx
        subq    $48, %rsp
        leaq    .L.str(%rip), %rdi
        callq   puts@PLT
        movl    $1024, %edi
        callq   fesetround@PLT
        movl    $1, %ebx
        testl   %eax, %eax
        jne     .LBB0_2
        fldl    .LCPI0_0(%rip)
        fstpt   (%rsp)
        callq   sqrtl@PLT
        fstpt   36(%rsp)
        leaq    .L.str.1(%rip), %rdi
        movsd   .LCPI0_1(%rip), %xmm0
        movaps  %xmm0, %xmm1
        movb    $2, %al
        callq   printf@PLT
        fldt    36(%rsp)
        fld     %st(0)
        fstpt   16(%rsp)
        fstpt   (%rsp)
        leaq    .L.str.2(%rip), %rdi
        xorl    %ebx, %ebx
        xorl    %eax, %eax
        callq   printf@PLT
.LBB0_2:
        movl    %ebx, %eax
        addq    $48, %rsp
        popq    %rbx
        retq

.L.str:
        .asciz  "setting rounding mode to FE_DOWNWARD"

.L.str.1:
        .asciz  "%.40f %.20a\n"

.L.str.2:
        .asciz  "%.40Lf %.20La\n"

最適化を有効にすると、clang でも間違った (fesetround() を無視した) 値が埋め込まれてしまった。

外部から読み込んだ値を計算する

前章の例では、ソースコードに埋め込んだ固定の値の平方根を求めたため、埋め込みがしやすいと判定されたかもしれない。
そこで、値を入力として受け取って計算を行うようにしてみた。

#include <stdio.h>
#include <math.h>
#include <fenv.h>

int main(void) {
    long long q;
    if (scanf("%lld", &q) != 1) return 1;
    double v = sqrt(q);
    long double vl = sqrtl(q);
    printf("%.40f %.20a\n", v, v);
    printf("%.40Lf %.20La\n", vl, vl);
#if defined(FE_DOWNWARD)
    puts("setting rounding mode to FE_DOWNWARD");
    if (fesetround(FE_DOWNWARD)) return 1;
#elif defined(FE_TOWARDZERO)
    puts("setting rounding mode to FE_TOWARDZERO");
    if (fesetround(FE_TOWARDZERO)) return 1;
#else
#error appropriate rounding mode unavailable
#endif
    double v2 = sqrt(q);
    long double vl2 = sqrtl(q);
    printf("%.40f %.20a\n", v2, v2);
    printf("%.40Lf %.20La\n", vl2, vl2);
    return 0;
}

すると、gcc 13.2.0 で -O3 -ffast-math オプションをつけてコンパイルしたとき、以下の出力が得られた。

67108865.0000000000000000000000000000000000000000 0x1.00000040000000000000p+26
67108864.9999999925494194030761718750000000000000 0x8.000001ffffffc0000000p+23
setting rounding mode to FE_DOWNWARD
67108864.9999999850988388061523437500000000000000 0x1.0000003ffffff0000000p+26
67108864.9999999925494194030761718750000000000000 0x8.000001ffffffc0000000p+23

clang 17.0.1 でも、同様のオプションで、以下の出力が得られた。

67108865.0000000000000000000000000000000000000000 0x1.00000040000000000000p+26
67108864.9999999925494194030761718750000000000000 0x8.000001ffffffc0000000p+23
setting rounding mode to FE_DOWNWARD
67108864.9999999850988388061523437500000000000000 0x1.0000003ffffff0000000p+26
67108864.9999999925494194030761718750000000000000 0x8.000001ffffffc0000000p+23

これらの環境では、外部から読み込んだ値の計算であれば、丸めモードを設定すれば double での計算でも適切な値が得られそうである。

結論

浮動小数点演算で確実に真の値以下の値が欲しいなら、丸めモードを設定しろ。
ただし、丸めモードを設定したからといって必ず真の値以下の値が得られるわけではない。
環境によっては丸めモードの設定が効かないこともあるし、そもそも丸めモードを設定できないかもしれない。

今回の用途 (平方根以下の最大の整数を求める) では、丸めモードの設定以外に、検算 (計算結果を2乗してもとの値以下になるかをチェック) して間違っていれば補正を行う方法でも、出力を改善できると考えられる。

2
1
0

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
Sign upLogin
2
1

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?