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進数展開などについて詳しく見ながら分析します。
目次
- 概要
- 計算が合わない場合の列挙
- 1.1のfloat32とfloat64における表現
- 1.1を10倍した時
- 1.1を一般の10の倍数倍した時
- 四捨五入でも計算が合わない場合
- いつ誤差が発生するのか
- float32の場合
- 実応用:ラブライブ!スクフェスにおけるスコア計算
- 結論
計算が合わない場合の列挙
まずは、冒頭の計算が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
一方、指数部の方を見てみると、それぞれ 01111111
と 01111111111
になっています。これらをunsigned intとして解釈して10進数に直すと127、1023になります。しかし、浮動小数の規格では、数値としての指数部の値に バイアス を加えたものを指数部ビットに入れています。float32、float64ではバイアスはそれぞれ127、1023なので、これらを差し引くと、数値としての指数部はいずれもちょうど0になることがわかります。
つまり、float32、float64では、それぞれ
という近似を使って1.1を表していることがわかります。ここで、、倍するのは、仮数部がそれぞれ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_32
と a_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)
が正しく計算される場合は、以下のようにして概ね説明できます。
- x = 2n * 10 (10進表現。n=0,1,2,…)の場合は、計算が合う。
- x = 2n * 10 + 10m (10進表現。m=1,2,…)を計算して行く。
- この計算結果の整数部を2進数表示した時の桁数を見ていく。桁数が繰り上がるまでは、仮数部に誤差が乗っていることから
ceil(1.1 * x)
の計算は合わない。桁が繰り上がった後は、繰り上がったことによって浮動小数による誤差が52ビット目から押し出され、丸め込みによって誤差が無くなり、ceil(1.1 * x)
の計算が正しく行われるようになる。 - 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の表現では、循環小数を丸め上げたことによる仮数部の誤差が、最終的に丸め上げる位に対して小さくなっているため起きていると考えられます。
実応用:ラブライブ!スクフェスにおけるスコア計算
この現象はスクフェスにおけるスコア計算という実用的な場面において発生しています。
応援タップSCOREアップの計算方法はまだ解が見つかっていないけど、とりあえず回復のパラメータアップの方の検証を進め、体力41~54まで判明。なお、応援タップSCOREアップは、SCOREの繰り下げ前の小数にアップ率を掛けて繰り上げると良さそうだけど、1つだけ合わない実測値ががががが… pic.twitter.com/9bsTnL18Nm
— まぐお@フライングバード=エリコヌスⅡ世 (@magu_maki) 2019年2月22日
この@magu_makiさんの解析の概要を説明すると、
- スクフェスにおいて、1回タップした時のスコアの理論式の解析を行っている
- 理論式の中に、ある整数(ベースとなるタップスコア)を1.1倍した上で切り上げをする という操作がある
- この操作が、入力が820の時は計算が合わない(上記ツイートの3枚目の画像の赤枠の数値)
- しかし、同じ10の倍数である970では計算が合っている(同画像の赤枠の1つ右、2つ下の枠の数値)
分析の結果、以下の事がわかりました。
- 冒頭で列挙した、
ceil(1.1 * x)
がfloat32とfloat64で計算が合わない場合を見ると、820では計算が合わないが、970では合うことがわかる - 計算が行われたExcelではfloat32で計算が行われている
- Excelのマクロ、つまりVBAではDouble型、つまりfloat64での計算が使用可能
- VBAを使って切り上げ処理を実装し、切り上げ処理のみfloat64で計算するようにした所、今まで計算が合っていなかった場所を含め、他の場所でも計算が合うようになった
結果的に、スクフェスのゲーム中における実際のスコアの測定値を説明するにはfloat64で切り上げ計算を行う必要があることがわかりました。
このようなビット単位で発生している現象が、2019年の今、ファミコンなどではなく、実際にiOSやAndroidで配布されているアプリの挙動を再現するという実用的な場面で起きているということに、今回は非常に驚かされました。
結論
ceil(1.1 * x)
の計算が、float32では正しく計算されるが、float64では正しく計算されない場合があることを確認した。この現象がPython、Java、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 によるスクフェスのスコア計算の解析において唯一説明できなかったデータ点が説明できるようになり、実際に観測されたデータが全て説明できるようになった。
*1:スクフェスの2019年2月のアップデートで変更された新たなスコア計算方法の解析を @magu_maki さんが行っていた所、一つだけ説明できないデータ点が浮上しました。様々な計算方法が試みられましたが、これに対して、私が計算がずれる原因は、浮動小数点由来の誤差ではないかと予想しました。その結果、浮動小数の誤差まで考慮すると、実際に観測されたデータが全て説明できるようになりました。