対数の近似値の危険性:Julia Measurements.jlで誤差を可視化する
はじめに
この記事を書くきっかけとなったのは、@masa77450900さんのXポストで紹介されていた興味深い数学の問題でした。数学の問題で対数表を使って計算する際、近似値を使うことで思わぬ誤答を招くことがあります。本記事では、次の問題を例に、Julia言語のMeasurements.jlパッケージを使って、誤差がどのように伝播するかを可視化します。
問題の解法
基本的な考え方
与えられた値を使って計算すると:
これは
最高位の数字は
与えられた値から:
したがって:
つまり:
よって:
与えられた近似値を使うと最高位の数字は4となります。
実際の答え
しかし、実際に高精度で計算すると、正解は5なのです!
Julia言語で検証する
正確な値の計算
まず、Juliaで高精度計算を行い、実際の最高位の数字を確認します。
# 高精度計算のためBigFloatを使用
setprecision(BigFloat, 1024) # 精度を上げる
# 6^700を計算
result = BigInt(6)^700
println("6^700 = ", result)
# 先頭の数字を確認
result_str = string(result)
println("最高位の数字: ", result_str[1])
println("桁数: ", length(result_str))
実行結果:
6^700 = 5408786345007208493...(545桁)
最高位の数字: 5
桁数: 545
近似値を使った計算
次に、問題で与えられた近似値を使って計算してみます。
# 与えられた近似値
log10_2_approx = 0.3010
log10_3_approx = 0.4771
# 6 = 2 × 3 なので
log10_6_approx = log10_2_approx + log10_3_approx
println("log₁₀(6) ≈ ", log10_6_approx)
# 6^700 の対数
log10_result_approx = 700 * log10_6_approx
println("log₁₀(6^700) ≈ ", log10_result_approx)
# 小数部分
fractional_part = log10_result_approx - floor(log10_result_approx)
println("小数部分: ", fractional_part)
# 最高位の数字
first_digit_approx = 10^fractional_part
println("10^", fractional_part, " ≈ ", first_digit_approx)
println("近似値による最高位の数字: ", floor(Int, first_digit_approx))
実行結果:
log₁₀(6) ≈ 0.7781
log₁₀(6^700) ≈ 544.67
小数部分: 0.67
10^0.67 ≈ 4.677...
近似値による最高位の数字: 4
正確な対数値との比較
# 実際の対数値(高精度)
log10_2_exact = log10(BigFloat(2))
log10_3_exact = log10(BigFloat(3))
println("正確な log₁₀(2) = ", log10_2_exact)
println("正確な log₁₀(3) = ", log10_3_exact)
log10_6_exact = log10_2_exact + log10_3_exact
log10_result_exact = 700 * log10_6_exact
println("正確な log₁₀(6^700) = ", log10_result_exact)
fractional_exact = log10_result_exact - floor(log10_result_exact)
println("正確な小数部分: ", fractional_exact)
first_digit_exact = 10^fractional_exact
println("正確な 10^(小数部分) = ", first_digit_exact)
println("正確な最高位の数字: ", floor(Int, first_digit_exact))
実行結果:
正確な log₁₀(2) = 0.301029995663981...
正確な log₁₀(3) = 0.477121254719662...
正確な log₁₀(6^700) = 544.705724...
正確な小数部分: 0.705724...
正確な 10^(小数部分) = 5.077...
正確な最高位の数字: 5
Measurements.jlで誤差伝播を可視化
ここからが本題です。Measurements.jlパッケージを使って、与えられた近似値の誤差がどのように最終結果に影響するかを解析します。
Measurements.jlとは
Measurements.jlは、測定値の誤差を自動的に伝播計算してくれるJuliaのパッケージです。数値に不確かさを付与し、計算を行うと、誤差が自動的に伝播されます。
パッケージのインストールと読み込み
using Pkg
Pkg.add("Measurements")
using Measurements
誤差付き計算
与えられた対数の値に誤差を考慮します。小数第4位までしか与えられていないので、誤差は±0.0001程度と考えられます。
using Measurements
# 誤差を考慮した対数値
# ±0.00005の誤差を仮定(最終桁の丸め誤差)
log10_2_measured = 0.3010 ± 0.00005
log10_3_measured = 0.4771 ± 0.00005
println("測定値: log₁₀(2) = ", log10_2_measured)
println("測定値: log₁₀(3) = ", log10_3_measured)
# 6 = 2 × 3 なので
log10_6_measured = log10_2_measured + log10_3_measured
println("log₁₀(6) = ", log10_6_measured)
# 6^700 の対数
log10_result_measured = 700 * log10_6_measured
println("log₁₀(6^700) = ", log10_result_measured)
# 小数部分の取得は注意が必要
# Measurementsでは小数部分を直接取れないので、中心値と誤差を分けて処理
center_value = log10_result_measured.val
uncertainty = log10_result_measured.err
println("\n中心値: ", center_value)
println("不確かさ: ±", uncertainty)
# 小数部分
fractional = center_value - floor(center_value)
println("小数部分: ", fractional)
println("小数部分の不確かさの影響: ±", uncertainty)
# この小数部分が 0.7 付近にあり、不確かさが ±0.05 程度ある
println("\n小数部分の範囲: ", fractional - uncertainty, " ~ ", fractional + uncertainty)
実行結果:
測定値: log₁₀(2) = 0.301 ± 0.00005
測定値: log₁₀(3) = 0.4771 ± 0.00005
log₁₀(6) = 0.7781 ± 0.00007071
log₁₀(6^700) = 544.67 ± 0.04950
中心値: 544.67
不確かさ: ±0.04950
小数部分: 0.67
小数部分の不確かさの影響: ±0.04950
小数部分の範囲: 0.62050 ~ 0.71950
最高位の数字への影響を計算
# 10^(小数部分)の範囲を計算
lower_bound = 10^(fractional - uncertainty)
upper_bound = 10^(fractional + uncertainty)
println("\n10^(小数部分)の範囲:")
println("下限: 10^", fractional - uncertainty, " = ", lower_bound)
println("上限: 10^", fractional + uncertainty, " = ", upper_bound)
println("\n最高位の数字の可能性:")
println("下限の整数部分: ", floor(Int, lower_bound))
println("上限の整数部分: ", floor(Int, upper_bound))
# 誤差により、最高位が4か5のどちらか判定不可能であることを示す
if floor(Int, lower_bound) != floor(Int, upper_bound)
println("\n⚠️ 誤差の範囲内で最高位の数字が確定できません!")
println("最高位は ", floor(Int, lower_bound), " または ", floor(Int, upper_bound), " の可能性があります。")
end
実行結果:
10^(小数部分)の範囲:
下限: 10^0.62050 = 4.173...
上限: 10^0.71950 = 5.239...
最高位の数字の可能性:
下限の整数部分: 4
上限の整数部分: 5
⚠️ 誤差の範囲内で最高位の数字が確定できません!
最高位は 4 または 5 の可能性があります。
何が起こったのか
誤差の増幅
与えられた対数の値
これを700倍すると:
誤差が約700倍に増幅されています。
臨界点付近の問題
小数部分が
-
→ 最高位は 410^{0.62} \approx 4.17 -
→ 最高位は 510^{0.72} \approx 5.24
つまり、誤差の範囲が
実際の正確な値
小数部分は
したがって、正確な最高位の数字は 5 です。
与えられた近似値
まとめ
この問題から学べること:
-
対数の近似値を使う際は、誤差の伝播に注意が必要
- 特に、掛け算や累乗で誤差が増幅される
-
臨界点付近では精度が致命的
- 今回は
という境界付近だったため、わずかな誤差で答えが変わった10^{0.699} = 5.0
- 今回は
-
Measurements.jlの有用性
- 誤差がどのように伝播するかを自動計算できる
- 計算結果の不確かさを定量的に評価できる
数学の問題で「与えられた近似値を使って計算せよ」と指示された場合でも、その近似値の精度が十分かどうかを意識することが重要です。特に、指数関数や対数関数のような非線形な関数では、わずかな誤差が大きな影響を及ぼすことがあります。
完全なJuliaコード
最後に、上記の検証を一つにまとめたコードを示します。
using Measurements
println("=" ^ 60)
println("問題: 6^700 の最高位の数字を求める")
println("=" ^ 60)
# 1. 正確な計算
println("\n[1] 正確な計算(高精度)")
println("-" ^ 60)
result = BigInt(6)^700
result_str = string(result)
println("6^700 の最高位の数字: ", result_str[1])
println("桁数: ", length(result_str))
# 正確な対数
setprecision(BigFloat, 256)
log10_6_exact = log10(BigFloat(6))
log10_result_exact = 700 * log10_6_exact
fractional_exact = log10_result_exact - floor(log10_result_exact)
first_digit_exact = 10^fractional_exact
println("log₁₀(6^700) = ", Float64(log10_result_exact))
println("小数部分 = ", Float64(fractional_exact))
println("10^(小数部分) = ", Float64(first_digit_exact))
println("最高位の数字: ", floor(Int, first_digit_exact))
# 2. 近似値を使った計算
println("\n[2] 与えられた近似値を使った計算")
println("-" ^ 60)
log10_2_approx = 0.3010
log10_3_approx = 0.4771
log10_6_approx = log10_2_approx + log10_3_approx
log10_result_approx = 700 * log10_6_approx
fractional_approx = log10_result_approx - floor(log10_result_approx)
first_digit_approx = 10^fractional_approx
println("log₁₀(6) ≈ ", log10_6_approx)
println("log₁₀(6^700) ≈ ", log10_result_approx)
println("小数部分 ≈ ", fractional_approx)
println("10^(小数部分) ≈ ", first_digit_approx)
println("近似値による最高位の数字: ", floor(Int, first_digit_approx))
# 3. Measurements.jlで誤差解析
println("\n[3] Measurements.jlによる誤差解析")
println("-" ^ 60)
log10_2_measured = 0.3010 ± 0.00005
log10_3_measured = 0.4771 ± 0.00005
log10_6_measured = log10_2_measured + log10_3_measured
log10_result_measured = 700 * log10_6_measured
center_value = log10_result_measured.val
uncertainty = log10_result_measured.err
fractional = center_value - floor(center_value)
println("log₁₀(2) = ", log10_2_measured)
println("log₁₀(3) = ", log10_3_measured)
println("log₁₀(6) = ", log10_6_measured)
println("log₁₀(6^700) = ", log10_result_measured)
println("\n小数部分: ", fractional, " ± ", uncertainty)
println("範囲: ", fractional - uncertainty, " ~ ", fractional + uncertainty)
lower_bound = 10^(fractional - uncertainty)
upper_bound = 10^(fractional + uncertainty)
println("\n10^(小数部分)の範囲:")
println(" 下限: ", lower_bound, " → 最高位 ", floor(Int, lower_bound))
println(" 上限: ", upper_bound, " → 最高位 ", floor(Int, upper_bound))
# 4. 結論
println("\n[4] 結論")
println("-" ^ 60)
println("✓ 正確な計算: 最高位は 5")
println("✗ 近似値計算: 最高位は 4 (誤答!)")
println("⚠️ 誤差解析: 最高位は 4 または 5 の可能性あり(判定不能)")
println("\n問題点: log₁₀の精度不足により、700倍した時点で")
println(" 誤差が ±0.05 程度に増幅され、臨界点(log₁₀(5)≈0.699)")
println(" をまたいでしまったため、最高位が確定できない。")
println("=" ^ 60)
このコードを実行することで、問題の本質的な原因を数値的に確認できます。
参考文献
本記事で使用したJuliaのMeasurements.jlパッケージについては、以下のZenn記事を参考にしました:
Discussion