Snow Ha凛ちゃん!

ラブライブに関するあれこれ

float32では正しく計算されるが、float64では正しく計算されない場合のある式:ceil(1.1*x) とそのラブライブ!スクフェスのスコア計算における影響

概要

色々な計算を行っていた所*1、次のようなfloat32とfloat64における一般的な現象に直面しました。それは、

ceil(1.1 * x)

という計算は、 float32だと計算が合うが、float64だと計算が合わない場合がある という現象です。つまり、 計算精度が高いはずのfloat64の方が計算が合いにくい 、という現象が起きました。

一言で言えば浮動小数の誤差で、1.1倍している関係から、基本的に計算が合わないのは10の倍数です。しかし、 10の倍数でも、中にはfloat64で計算の合う場合もある ことがわかりました。

例えば、x

  • 230の時、float32では ceil(1.1 * x)253 となるのに対し、float64では 254 と正しく計算されない
  • しかし540の時、 float32でも 594 、float64でも 594 と、正しく計算される
  • この現象は、Python(numpy使用)、C、Javaで確認される

試しに0から1500未満の整数でこの現象が起きる(つまり、float64で計算が合わない)ような数を列挙すると、

50 90 100 110 170 180 190 200 210 220 230 330 340 350 360 370 380 390 400 410 420 430 440 450 460 650 660 670 680 690 700 710 720 730 740 750 760 770 780 790 800 810 820 830 840 850 860 870 880 890 900 910 920 930 1290 1300 1310 1320 1330 1340 1350 1360 1370 1380 1390 1400 1410 1420 1430 1440 1450 1460 1470 1480 1490

のように、一見不規則な列になります。

もっとも驚くべきことは、このような現象が 実用上の場面 で起きているということです。具体的には、後に紹介する通り、@magu_maki さんの観測データを用いて、スクフェスの点数計算を手元でシミュレーションしようとした際、float32とfloat64で計算結果が異なり、float64で計算しないと挙動が再現できなかった、という点です。

ここでは、

  • なぜfloat64で計算が合わない場合があるのか
  • 計算が合わないのは基本的に入力が10の倍数の時であるように見えるが、なぜ10の倍数の中でも、計算が合う時があるのか

について、float型のデータ表現方法や2進数展開などについて詳しく見ながら分析します。

目次

計算が合わない場合の列挙

まずは、冒頭の計算がfloat32とfloat64の場合で異なるような正の整数xを列挙します。試しに、0から1500未満の整数について調べます。

以下のコードを使います。

public class DoubleFloat {
    public static void main(String[] args) {
        float a_32 = 1.1f, b_32;
        double a_64 = 1.1d, b_64;
        for (int i = 0; i < 1500; i++) {
            b_32 = i;
            b_64 = i;
            if (Math.ceil(a_32*b_32) != Math.ceil(a_64*b_64)) {
                System.out.printf("%d %f %f\n", i, Math.ceil(a_32*b_32), Math.ceil(a_64*b_64));
            }        
        }
    }
}

すると、以下の結果を得ます。

50 90 100 110 170 180 190 200 210 220 230 330 340 350 360 370 380 390 400 410 420 430 440 450 460 650 660 670 680 690 700 710 720 730 740 750 760 770 780 790 800 810 820 830 840 850 860 870 880 890 900 910 920 930 1290 1300 1310 1320 1330 1340 1350 1360 1370 1380 1390 1400 1410 1420 1430 1440 1450 1460 1470 1480 1490

Math.ceil の結果も出力するように、9行目を

System.out.printf("%d %f %f\n", i, Math.ceil(a_32*b_32), Math.ceil(a_64*b_64));

のように変更すると、以下のような結果になります。

50 55.000000 56.000000
90 99.000000 100.000000
100 110.000000 111.000000
110 121.000000 122.000000
170 187.000000 188.000000
180 198.000000 199.000000
190 209.000000 210.000000
200 220.000000 221.000000
210 231.000000 232.000000
(中略)
1450 1595.000000 1596.000000
1460 1606.000000 1607.000000
1470 1617.000000 1618.000000
1480 1628.000000 1629.000000
1490 1639.000000 1640.000000

同様の結果が、CおよびPythonでも、以下のコードを使うことで再現できます。

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

int main(int argc, char const *argv[])
{
    float a_32 = 1.1, b_32;
    double a_64 = 1.1, b_64;

    for (int i = 0; i < 1500; i++) {
        b_32 = i;
        b_64 = i;
        if (ceil(a_32*b_32) != ceil(a_64*b_64)) {
            printf("%d %f %f\n", i, ceil(a_32*b_32), ceil(a_64*b_64));
        }
    }

    return 0;
}
import numpy as np
for x in range(1499):
    a_32 = np.array(1.1, dtype=np.float32)
    b_32 = np.array(x, dtype=np.float32)
    ceil_32 = np.ceil(a_32*b_32)

    a_64 = np.array(1.1, dtype=np.float64)
    b_64 = np.array(x, dtype=np.float64)
    ceil_64 = np.ceil(a_64*b_64)
    
    if ceil_32 != ceil_64:
        print(x, ceil_32, ceil_64)

さて、先程の実行結果を見ると、基本的に10の倍数において、float32とfloat64で結果が1違う、ということがわかります。しかし、 10の倍数の中でも、計算結果が合う場合がある ことがわかります。

以下では、その原因を分析していきます。

1.1のfloat32とfloat64における表現

まず、float32とfloat64における1.1の、浮動小数型としての実際の2進数表現を出力します。

以下のJavaのコードを使います。

public class Sif {

    public static String float2string (float x) {
        return Integer.toBinaryString(Float.floatToRawIntBits(x));
    }

    public static String double2string (double x) {
        return Long.toBinaryString(Double.doubleToRawLongBits(x));
    }

    public static void main(String[] args) {
        float a_32 = 1.1f;
        double a_64 = 1.1d;
        System.out.printf("float32 (binary):\n");
        System.out.printf("1.1      : 00%s\n", float2string(a_32));
        System.out.printf("           ^~~~~~~~~=======================\n");
        System.out.printf("float64 (binary):\n");
        System.out.printf("1.1      : 00%s\n", double2string(a_64));
        System.out.printf("           ^~~~~~~~~~~~====================================================\n");
        System.out.printf("           ^: 符号部\n");
        System.out.printf("           ~: 指数部\n");
        System.out.printf("           =: 仮数部\n");
    }
}

すると、以下のような結果を得ます。

float32 (binary):
1.1      : 00111111100011001100110011001101
           ^~~~~~~~~=======================
float64 (binary):
1.1      : 0011111111110001100110011001100110011001100110011001100110011010
           ^~~~~~~~~~~~====================================================
           ^: 符号部
           ~: 指数部
           =: 仮数部

ここで、浮動小数の規格*2を見ると、float32、float64では、それぞれ後ろ23ビット、52ビットが仮数部ビットになっていることがわかります。つまり、float32, float64における仮数部ビットは、それぞれ

  • 00011001100110011001101
  • 0001100110011001100110011001100110011001100110011010

です。

ここで、浮動小数では、 仮数部の先頭に1をつけたものが、実際の数値としての仮数部になります。 これは、2進数では0以外の数の先頭の桁は必ず1になることから、その桁を省略しても必ず数を復元することができることに由来しています。そのため、実際の仮数部は

  • 100011001100110011001101
  • 10001100110011001100110011001100110011001100110011010

になります。

続いて、これらをわかりやすく扱うため10進数の整数に直します。そのために、以下のPythonコードを使います。

print(int("100011001100110011001101", 2))
print(int("10001100110011001100110011001100110011001100110011010", 2))

すると、以下の結果を得ます。

9227469
4953959590107546

一方、指数部の方を見てみると、それぞれ 0111111101111111111 になっています。これらをunsigned intとして解釈して10進数に直すと127、1023になります。しかし、浮動小数の規格では、数値としての指数部の値に バイアス を加えたものを指数部ビットに入れています。float32、float64ではバイアスはそれぞれ127、1023なので、これらを差し引くと、数値としての指数部はいずれもちょうど0になることがわかります。

つまり、float32、float64では、それぞれ


1.1 \approx 9227469 * 2^{-23} = 1.10000002384185791015625


1.1 \approx 4953959590107546 * 2^{-52} = 1.100000000000000088817841970012523233890533447265625

という近似を使って1.1を表していることがわかります。ここで、2^{-23}2^{-52}倍するのは、仮数部がそれぞれ23ビット、52ビットで表現されているためです。

このような誤差が出る原因は、1.1(10進表記)は10進数では循環小数にならないが、2進数では循環小数になることにあります。実際、1.1の仮数部を見てみると、00011001100…のように、000の後に1100が続いていることがわかります。浮動小数における1.1は、この循環小数を途中で打ち切ったものとして表現されているのが、この誤差を生む原因になっています。

1.1を10倍した時

さて、冒頭の計算が合わないのは、計算した範囲だと全て入力が10の倍数である場合でした。そこで、先程計算した仮数部が10倍された時の結果を見ます。そのために、仮数部の2倍、8倍、10倍を求めます。

以下のコードを使います。

a_32 = 9227469
print("a_32:   {0:30b}".format(a_32))
print()
print("a_32*2: {0:30b}".format(a_32*2))
print("a_32*8: {0:30b}".format(a_32*8))
print("a_32*10:{0:30b}".format(a_32*10))
print(" "*(11+24) + "^")
print(" "*11 + "".join(map(str, range(10)))*3)
print(" "*11 + (" "*9).join(map(str, range(3))))
print()

a_64 = 4953959590107546
print("a_64:   {0:60b}".format(a_64))
print()
print("a_64*2: {0:60b}".format(a_64*2))
print("a_64*8: {0:60b}".format(a_64*8))
print("a_64*10:{0:60b}".format(a_64*10))
print(" "*(12+53) + "^")
print(" "*12 + "".join(map(str, range(10)))*6)
print(" "*12 + (" "*9).join(map(str, range(6))))

すると、以下のようになります。

a_32:         100011001100110011001101

a_32*2:      1000110011001100110011010
a_32*8:    100011001100110011001101000
a_32*10:   101100000000000000000000010
                                   ^
           012345678901234567890123456789
           0         1         2

a_64:          10001100110011001100110011001100110011001100110011010

a_64*2:       100011001100110011001100110011001100110011001100110100
a_64*8:     10001100110011001100110011001100110011001100110011010000
a_64*10:    10110000000000000000000000000000000000000000000000000100
                                                                 ^
            012345678901234567890123456789012345678901234567890123456789
            0         1         2         3         4         5

上記において、 a_32a_64 が、それぞれ1.1のfloat32、float64における仮数部を整数に直したものを表しています。

まず、2進数において2倍するという操作は0を末尾に1つ付け加える事に対応し、8倍するという操作は0を末尾に3つ付け加えるという操作に対応します。そして、10倍するという操作はそれらの結果を足す事に相当します。

それらを踏まえた上で、上記のように1.1のfloat32、float64における仮数部を10倍すると、循環部分の1100が綺麗に打ち消し合い、末尾付近の桁のみ誤差が残ることがわかります。上記の結果において、 ^ で示した位置が、最終的に仮数部ビットに含まれない最後の桁になっています。末尾の誤差の位と仮数部ビットに含まれない位の位置関係を見ると、 float64の場合の方が、float32に比べて、誤差が仮数部ビットに含まれない桁に近い ことがわかります。

これについて、float64の場合により詳しい計算を行っていきます。

1.1を一般の10の倍数倍した時

さて、1.1を一般の10の倍数倍した時に、float64の仮数部の途中計算の様子を見るために、1.1のfloat64における仮数部を10の倍数倍した結果の2進数表記を見てみましょう。

以下のコードを使います。

import numpy as np
a = 4953959590107546

def ceil32_and_64_are_different(x):
    a_32 = np.array(1.1, dtype=np.float32)
    b_32 = np.array(x, dtype=np.float32)
    ceil_32 = np.ceil(a_32*b_32)

    a_64 = np.array(1.1, dtype=np.float64)
    b_64 = np.array(x, dtype=np.float64)
    ceil_64 = np.ceil(a_64*b_64)
    
    return ceil_32 != ceil_64

for i in [10,20,30,40,50,60,70,80,90,100,110,120,130,140,150,160,170,180,310,320,330,340,630,640,650,660,920,930,940,950,1270,1280,1290,1300]:
    o = ("x" if ceil32_and_64_are_different(i) else " ") + " "
    binstr = "{0:b}".format(a*i)
    print(o + "{0:4d}:{1} {2}".format(i, binstr[:54], binstr[54:]))
print(" "*(7+53) + "^")
print(" "*7 + ("".join(map(str, range(10)))*6)[:54])
print(" "*7 + (" "*9).join(map(str, range(6))))

すると、以下のような結果を得ます。

    10:101100000000000000000000000000000000000000000000000001 00
    20:101100000000000000000000000000000000000000000000000001 000
    30:100001000000000000000000000000000000000000000000000000 1100
    40:101100000000000000000000000000000000000000000000000001 0000
x   50:110111000000000000000000000000000000000000000000000001 0100
    60:100001000000000000000000000000000000000000000000000000 11000
    70:100110100000000000000000000000000000000000000000000000 11100
    80:101100000000000000000000000000000000000000000000000001 00000
x   90:110001100000000000000000000000000000000000000000000001 00100
x  100:110111000000000000000000000000000000000000000000000001 01000
x  110:111100100000000000000000000000000000000000000000000001 01100
   120:100001000000000000000000000000000000000000000000000000 110000
   130:100011110000000000000000000000000000000000000000000000 110100
   140:100110100000000000000000000000000000000000000000000000 111000
   150:101001010000000000000000000000000000000000000000000000 111100
   160:101100000000000000000000000000000000000000000000000001 000000
x  170:101110110000000000000000000000000000000000000000000001 000100
x  180:110001100000000000000000000000000000000000000000000001 001000
   310:101010101000000000000000000000000000000000000000000000 1111100
   320:101100000000000000000000000000000000000000000000000001 0000000
x  330:101101011000000000000000000000000000000000000000000001 0000100
x  340:101110110000000000000000000000000000000000000000000001 0001000
   630:101011010100000000000000000000000000000000000000000000 11111100
   640:101100000000000000000000000000000000000000000000000001 00000000
x  650:101100101100000000000000000000000000000000000000000001 00000100
x  660:101101011000000000000000000000000000000000000000000001 00001000
x  920:111111010000000000000000000000000000000000000000000001 01110000
x  930:111111111100000000000000000000000000000000000000000001 01110100
   940:100000010100000000000000000000000000000000000000000000 101111000
   950:100000101010000000000000000000000000000000000000000000 101111100
  1270:101011101010000000000000000000000000000000000000000000 111111100
  1280:101100000000000000000000000000000000000000000000000001 000000000
x 1290:101100010110000000000000000000000000000000000000000001 000000100
x 1300:101100101100000000000000000000000000000000000000000001 000001000
                                                            ^
       012345678901234567890123456789012345678901234567890123
       0         1         2         3         4         5

ここで、 x はfloat32とfloat64で結果が合わない場合、 ^ はfloat64において仮数部に含まれない最初の桁、つまり仮数部ビットの53ビット目に相当する位を表しています。浮動小数の規格では仮数部の先頭の位の1は省略されるので、実際の位としては上から数えて54桁目に相当します。読みやすさのため、53ビット目以降をスペースで区切っています。

さて、実際の乗算の計算では、こうして得られた仮数部をfloat64表現に直すために、最後にこの仮数部の丸め込み処理を行います。仮数部は52ビットあるので、 ^ のある53ビット目で丸め込み をすることで52ビット目を決定します。(正確には、基本は四捨五入を行うが、丁度中間の値のみ切り下げを行う「最近接丸め」を使っているようです。これについては、次の節で説明します。また、仮数部の冒頭の1を省略した状態での53ビット目を表しています。)つまり、 53ビット目が1ならば52ビット目が1に、53ビット目が0ならば52ビット目が0 になるはずです。

試しに930、940の場合にどのように丸め込まれるかを見てみると、

x  930:111111111100000000000000000000000000000000000000000001 01110100
x  930:11111111110000000000000000000000000000000000000000001

   940:100000010100000000000000000000000000000000000000000000 101111000
   940:10000001010000000000000000000000000000000000000000000
                                                            ^
       012345678901234567890123456789012345678901234567890123
       0         1         2         3         4         5

のようになります。すると、 930の場合は仮数部の一番低い位に誤差1が乗っているのに対し、940の場合は仮数部には誤差が乗っていない ことがわかります。これが、 float64において ceil(1.1 * x) の計算が合わない場合の原因 、および float32とfloat64で計算が合う場合と合わない場合の違い になっています。

四捨五入でも計算が合わない場合

さて、^ のある53ビット目で四捨五入 を行って52ビット目を決定すると書きましたが、これだけでは説明できない場合が存在します。それは、10,20,40,80,160,320,640,1280の場合で、これらの場合では53ビット目が1になっているにもかかわらず、ceilにおいて正しく結果が計算されていることがわかります。

これは、おそらく浮動小数の規格において、丸め方として 最近接丸め を使っているためであると予想されます。Wikipediaの浮動小数の規格の記事における「浮動小数点数の丸め」の項目を見ると、最近接丸めについて以下のように説明されています。

最も近くの表現できる値へ丸める。表現可能な2つの値の中間の値であったら、一番低い仮数ビット(桁)が0になるほうを採用する。これは二進での標準動作かつ十進でも推奨となっている。

ここで、これが「二進での標準動作かつ十進でも推奨となっている」と説明されていることから、浮動小数の掛け算の結果においても同様の処理が行われていると考えることができます。

さて、上記の結果において10,20,40,80,160,320,640,1280の場合のみを見てみると、仮数部の 54桁目以降が全て0 になっていることがわかります。

    10:101100000000000000000000000000000000000000000000000001 00
    20:101100000000000000000000000000000000000000000000000001 000
    40:101100000000000000000000000000000000000000000000000001 0000
    80:101100000000000000000000000000000000000000000000000001 00000
   160:101100000000000000000000000000000000000000000000000001 000000
   320:101100000000000000000000000000000000000000000000000001 0000000
   640:101100000000000000000000000000000000000000000000000001 00000000
  1280:101100000000000000000000000000000000000000000000000001 000000000
                                                            ^
       012345678901234567890123456789012345678901234567890123
       0         1         2         3         4         5

四捨五入を行う桁である53桁目が1で、それ以降の54桁目以降が全て0ということは、これらの場合において、仮数部は 丁度四捨五入の中間の位置 にあることになります。最近接丸みの説明では、この場合は 仮数部が0となるように丸める 、つまり 切り捨てを行う と書いてあるので、53桁目が1であるにも関わらず、切り捨てを行うことで最終的な結果を得ます。

その結果、 これらの場合では仮数部の最終桁は0になる ことから、浮動小数由来の誤差が乗らず、 ceilが正しく計算される ことになります。

つまり、丸め込みによる仮数部の52桁目の決定方法は、

  • 53桁目が0ならば、52桁目は0
  • 53桁目が1の場合、
    • 54桁目以降が全て0ならば、52桁目は0
    • そうでない場合は、通常の四捨五入の通り、52桁目は1

という処理を行うことで決定できると予想されます。

この方法を用いると、上記の0から1499の場合の ceil(1.1 * x) の計算結果を全て説明することができます。

いつ誤差が発生するのか

float64における1.1の仮数部の表現を10倍した時の数の表現を見ると、

a_64*10:    10110000000000000000000000000000000000000000000000000100
                                                                 ^
            012345678901234567890123456789012345678901234567890123456789
            0         1         2         3         4         5

のように、誤差は1桁のみ存在することがわかります。 そのため、1.1 * m * 10 という数の誤差の桁数は、mの桁数そのものになることがわかります。

また、float32とfloat64において結果が合わない場合の一覧を見ると、計算結果の整数部の桁が繰り上がった時に、 誤差の桁が52ビット目から押し出されることで計算が合うようになっている ことがわかります。

以上と先程の丸め込みの方法を踏まえると、float64で ceil(1.1 * x) が正しく計算される場合は、以下のようにして概ね説明できます。

  1. x = 2n * 10 (10進表現。n=0,1,2,…)の場合は、計算が合う。
  2. x = 2n * 10 + 10m (10進表現。m=1,2,…)を計算して行く。
  3. この計算結果の整数部を2進数表示した時の桁数を見ていく。桁数が繰り上がるまでは、仮数部に誤差が乗っていることから ceil(1.1 * x) の計算は合わない。桁が繰り上がった後は、繰り上がったことによって浮動小数による誤差が52ビット目から押し出され、丸め込みによって誤差が無くなり、 ceil(1.1 * x) の計算が正しく行われるようになる。
  4. 2n * 10 + 10m == 2n+1 * 10 (10進表現)に到達したら、1.に戻る。

float32の場合

最後に、同様の計算をfloat32において行った場合について見ます。

以下のコードを使います。

import numpy as np
a = 9227469

for i in [10,20,30,40,50,60,70,80,90,100,110,120,130,140,150,160,170,180,310,320,330,340,630,640,650,660,920,930,940,950,1270,1280,1290,1300]:
    binstr = "{0:b}".format(a*i)
    print("  {0:4d}:{1} {2}".format(i, binstr[:25], binstr[25:]))
print(" "*(7+24) + "^")
print(" "*7 + ("".join(map(str, range(10)))*6)[:25])
print(" "*7 + (" "*9).join(map(str, range(3))))

すると、以下のような結果を得ます。

    10:1011000000000000000000000 10
    20:1011000000000000000000000 100
    30:1000010000000000000000000 0110
    40:1011000000000000000000000 1000
    50:1101110000000000000000000 1010
    60:1000010000000000000000000 01100
    70:1001101000000000000000000 01110
    80:1011000000000000000000000 10000
    90:1100011000000000000000000 10010
   100:1101110000000000000000000 10100
   110:1111001000000000000000000 10110
   120:1000010000000000000000000 011000
   130:1000111100000000000000000 011010
   140:1001101000000000000000000 011100
   150:1010010100000000000000000 011110
   160:1011000000000000000000000 100000
   170:1011101100000000000000000 100010
   180:1100011000000000000000000 100100
   310:1010101010000000000000000 0111110
   320:1011000000000000000000000 1000000
   330:1011010110000000000000000 1000010
   340:1011101100000000000000000 1000100
   630:1010110101000000000000000 01111110
   640:1011000000000000000000000 10000000
   650:1011001011000000000000000 10000010
   660:1011010110000000000000000 10000100
   920:1111110100000000000000000 10111000
   930:1111111111000000000000000 10111010
   940:1000000101000000000000000 010111100
   950:1000001010100000000000000 010111110
  1270:1010111010100000000000000 011111110
  1280:1011000000000000000000000 100000000
  1290:1011000101100000000000000 100000010
  1300:1011001011000000000000000 100000100
                               ^
       0123456789012345678901234
       0         1         2

float32では仮数部は23ビットあるので、 ^ で示した24ビット目(冒頭の1を省略して数えた場合)に対して丸め込みを行った結果が最終的な計算結果となります。しかし、上記の場合では 23ビット目はいずれも0 なので、 これらの場合では常にceilの計算が正しく行われる ことがわかります。

実際、他の場合でもこれが成り立つか調べるため、

import numpy as np
for x in range(200000):
    a_32 = np.array(1.1, dtype=np.float32)
    b_32 = np.array(x, dtype=np.float32)
    ceil_32 = np.ceil(a_32*b_32)
    ceil = x + x // 10 + (1 if x % 10 > 0 else 0)

    if ceil_32 != ceil:
        print(x, ceil_32, ceil)

のようなコードを実行すると、0から200000未満の全ての x において、float32では ceil(1.1 * x) が正しく計算されていることがわかります。

これは、この記事の冒頭付近の「1.1のfloat32とfloat64における表現」の節で見た通り、float32における1.1の表現では、循環小数を丸め上げたことによる仮数部の誤差が、最終的に丸め上げる位に対して小さくなっているため起きていると考えられます。

実応用:ラブライブ!スクフェスにおけるスコア計算

この現象はスクフェスにおけるスコア計算という実用的な場面において発生しています。

この@magu_makiさんの解析の概要を説明すると、

  • スクフェスにおいて、1回タップした時のスコアの理論式の解析を行っている
  • 理論式の中に、ある整数(ベースとなるタップスコア)を1.1倍した上で切り上げをする という操作がある
  • この操作が、入力が820の時は計算が合わない(上記ツイートの3枚目の画像の赤枠の数値)
  • しかし、同じ10の倍数である970では計算が合っている(同画像の赤枠の1つ右、2つ下の枠の数値)

分析の結果、以下の事がわかりました。

結果的に、スクフェスのゲーム中における実際のスコアの測定値を説明するにはfloat64で切り上げ計算を行う必要があることがわかりました。

このようなビット単位で発生している現象が、2019年の今、ファミコンなどではなく、実際にiOSAndroidで配布されているアプリの挙動を再現するという実用的な場面で起きているということに、今回は非常に驚かされました。

結論

  • ceil(1.1 * x) の計算が、float32では正しく計算されるが、float64では正しく計算されない場合があることを確認した。この現象がPythonJava、Cで発生することを確認した。
  • この現象は基本的に入力が10の倍数の時に発生するようである、と予想できることを見た。しかし、10の倍数の中でも正しく計算される場合が存在することを確認した。
  • 10の倍数の中でも正しく計算される場合がなぜ存在するかについて、float32とfloat64における 1.1 のデータ表現における仮数部の計算を分析することで、詳しい説明を行った。特に、float64で x が10,20,40,80,160,320,640,1280の場合は、浮動小数の最近接丸めを考慮する必要がある可能性があることを指摘した。
  • float32の場合で同様の計算を行うと、0から200000未満の整数の入力では、浮動小数由来の誤差が発生せず、 ceil(1.1 * x) が正しく計算されることを確認した。
  • 以上のような浮動小数の誤差を考慮して計算を行った結果、@magu_maki によるスクフェスのスコア計算の解析において唯一説明できなかったデータ点説明できるようになり実際に観測されたデータが全て説明できるようになった