病的な計算式?(『ソフト・エッジ』から)

4621053833中島 震・みわ よしこ『ソフト・エッジ: ソフトウェア開発の科学を求めて』(丸善出版, 2013)で、数値計算の面白い例が紹介されていました。(参考文献リストあり。索引なし。)

精度を向上させれば近似はよくなるというのが、自然な考え方でしょう。でも、この直観に当てはまらない「病的な計算式」の例があります。(p.132)

ということで、f(a, b)=(333.75-a**2)*b**6+a**2*(11*a**2*b**2-121*b**4-2)+5.5*b**8+a/(2*b)としたときの、f(77617, 33096)が紹介されています。単精度・倍精度・4倍精度の浮動小数点数での計算結果は約1.17で、正しい値(約-0.827)とは符号も違ってしまうとのことです。

これだけ読むと、精度を向上させても近似はよくならないと思うかもしれませんが、そんなことはないはずです。実際、精度をある程度高くすれば、正しく計算できます。

Pythonで試します(ideone.comでも動きます)。

from decimal import getcontext, Decimal
for n in range(16, 40):
  getcontext().prec = n
  a = Decimal("77617")
  b = Decimal("33096")
  print(n, (Decimal("333.75")-a**2)*b**6+a**2*(11*a**2*b**2-121*b**4-2)+Decimal("5.5")*b**8+a/(2*b))
16 1.000000000000000E+21
17 -4.0000000000000000E+20
18 -2.00000000000000000E+19
19 2000000000000000001
20 -99999999999999998.827
21 1.17260394005317863186
22 -999999999999998.8273961
23 100000000000001.17260394
24 10000000000001.1726039401
25 -3999999999998.827396059947
26 -99999999998.827396059946821
27 -19999999998.8273960599468214
28 -999999998.8273960599468213681
29 100000001.17260394005317863186
30 20000001.1726039400531786318588
31 -999998.8273960599468213681411651
32 300001.17260394005317863185883490
33 -9998.82739605994682136814116509548
34 -1998.827396059946821368141165095480
35 -198.82739605994682136814116509547982
36 21.1726039400531786318588349045201837
37 -0.827396059946821368141165095479816292
38 -0.8273960599468213681411650954798162920
39 -0.82739605994682136814116509547981629200

このようにPythonでは、getcontext().precの値が37以上なら適切な結果が得られます。

Rubyでも、ちょっと工夫すれば適切な結果が得られます(ideone.comで確認)

Windows付属の電卓では、何の工夫も要りません(WIndows 7で確認)。xyのショートカットキーが「Y」なので、「(333.75-77617Y2)*33096Y6+77617Y2*(11*77617Y2*33096Y2-121*33096Y4-2)+5.5*33096Y8+77617/(2*33096)=」と打てば、適切な結果が得られます。

Mac OS Xでも計算できます(10.8.4で確認。できるようになったのは比較的最近のことだと思いますが)。「計算機」で計算する方法はよくわかりませんが、Spotlightに「(333.75-77617^2)*33096^6+77617^2*(11*77617^2*33096^2-121*33096^4-2)+5.5*33096^8+77617/(2*33096)」と入力すれば、適切な結果が得られます。

Windowsの電卓は32桁までの入力にしか対応していないので、230928922463004537535516678003369の平方根を求めたいというようなときに不便ですが、Spotlightなら、sqrt(230928922463004537535516678003369)とすればとりあえず計算できます。

参考:電卓に求められるコト

本稿執筆時点のWolframAlphaでは、まったく違う結果が2つ返ってきて、片方はあっていました

Googleはだめでした

Googleの例は別として、精度を上げてもうまくいかないようにするには、もっと「病的」でなければなりません。

Excelを使っていると採点ミスで謝罪することになる?

事例が見つかるたびに騒ぎになっているようです。今回は、「8.2 – 7.2」の結果が1より小さくなってしまうというものです。最初に非難されたのは例によってExcelのようですが、他の環境でもこういう問題は起こります(C言語の例Rubyの例)。

#include <stdio.h>

int main() {
  double x = 8.2 - 7.2;
  if (x >= 1) printf("x >= 1");
  else printf("x < 1");
  return 0;
}

「コンピュータは数値を2進数としてしか表せなくて・・・」という説明をよく見ますが、この説明は間違いです。これでは、「コンピュータは正しい計算ができない」という誤解を生んでしまいます。

10進数の1, 2, 3, …を2進数の1, 10, 11, …で表すという具合に、数を2進数で直接表現しようとするとこういう問題が起きますが、別の表現方法を用いれば、このような問題は起きません。実際、Mathematicaは正しく計算します。

With[{x = 8.2 - 7.2}, If[x >= 1, "x >= 1", "x < 1"]]

x >= 1

コンパイルしても大丈夫です。(UMMでは動かないようです。)

f = Compile[{a, b}, With[{x = a - b}, Print[If[x >= 1, "x >= 1", "x < 1"]]]];
f[8.2, 7.2]

x >= 1

C言語やRubyではどうすればいいのかというと、許容可能な誤差を決めておいて、その範囲内はよしとするのです。

さらに、この許容可能な誤差は、扱う数によって変えなければならないことに注意してください。以下の例では、許容可能な誤差を1より大きくしなければなりません(実行結果)。

#include <stdio.h>

int main() {
  double a = 9007199254740992;
  double b = a + 1;
  if (a == b) printf("a==b");
  return 0;
}

4891006269こういうことは、『Microsoft Visual C++入門』の3.2.11項「浮動小数点についての注意」にちゃんと書いてあったりします。ですから、入門を終えたプログラマはこういうことは知っています。

Excelのような、こうことは知らないであろう人が使うソフトウェアでは、こういう工夫をしなくても正しく計算できるようにしておいてほしいものです。フェルマーの最終定理の「反例」のような冗談を言っているうちはいいのですが、「x >= 0.6 なら合格」というようなことって、よくやられているのではないでしょうか(ボーダーの学生が間違って不合格になる危険があります)。間違いが公になるときには、「採点ミス」になるのでしょうか。

数独で見るRuby(とMathematica)のパワーと表現力

参考:数独の平凡な解法(C言語Mathematica

プログラミング言語 Ruby (大型本)Rubyのバイブル『プログラミング言語 Ruby』の第1.4節では、「Rubyプログラムが実際にはどのようなものかというイメージをもっとよくつかめるように(p.18)」数独を解くRubyプログラムが紹介されています(ソースコードは原著のサポートサイトにあります)。曰く、

コメントと空行を取り除くと、ちょうど129行のコードが残る。これは、単純な力任せのアルゴリズムに頼るわけではないオブジェクト指向の数独ソルバとしてはまずまずの長さだ。このサンプルは、Rubyのパワーと表現力をよく示していると思うがどうだろうか。(p.25)

どうだろうかって、このサンプルは、Rubyのパワーと表現力を示すものとしてはふさわしくないんじゃないでしょうか。100行以上書かないと数独が解けないなんて、Rubyのイメージダウンにならなければいいのですが。

Rubyが手続き型プログラミングのためのすばらしい言語であることは否定しませんが、所詮は手続き型なので、「欲しい結果を得るための処理」を書かなければならない場合には便利でも、「欲しい結果の性質」を書けばよい言語に比べれば、プログラマの負担は大きいはずです。

Ruby本のコードには入力からゴミを取り除くような本質的ではない処理も含まれているので、行数で比較することが最善の評価法というわけではありませんが、数独はSQLなら1行ですし、その1行を生成させるJavaのコードでも70行程度です。行数というのはそんなにいい指標ではないわけですが(Pythonでは1行という話もありますし)、一つの目安ということで。

適材適所という観点から言って、数独はProlog(やSQL)のような言語で解くのがいいと思うのですが、どうしても手続き型の言語を使いたいという向きのために、Ruby本よりはシンプルな方法(10行)を紹介しましょう。実装にはMathematicaを使いますが、関数型の記述ができる言語になら簡単に翻訳できるでしょう。引用箇所の「オブジェクト指向の」は無視します。

単純探索のテンプレートを使います。まず、数独に限らず、一般的な探索に使える述語を定義します。

search[x_] := Or[goal@x, deepen@x]
scanOr[f_, x_] := Catch[Scan[If[f@#, Throw@True] &, x]; False] (* 深さ優先探索 *)

数独の解を定義します。空白(0)の部分が無い(Not@MemberQ)ものが解です。

goal[x_] := And[Not@MemberQ[x, 0, 2], report@x]
report[x_] := (Sow@TableForm@x; True) (* 解の表示 *)

力任せの方法

引用箇所には「単純な力任せのアルゴリズムに頼るわけではない」とあります。Ruby本のコードは、未定のマス目に入りうる数字を数独のルール通りに決めているだけなので、「力任せ」だと言っていいと思うのですが、ここで言う「力任せ」の方法とは次のようなものでしょう。

力任せの方法:「空白部分に1から9までの数字を入れてルールをチェック」の繰り返し

述語deepenは、空白部分を調べ(Position)、1から9までの数字(Range[1,9])を入れてみて、ルール(test)に違反していなければ(Select)探索(search)を進めます。

deepen[x_] := scanOr[search, Select[ReplacePart[x, First@Position[x, 0] -> #] & /@ Range[1, 9], test]]

ルール(test)は、各行(x)と各列(Transpose@x)、各ブロック(Flatten…)を集計(Tally)し、1以上の部分({_, a_ /; a > 1})がないことを確認します。

test[x_] :=
 And @@ Map[Not[MemberQ[Tally[Select[#, (# != 0) &]], {_, a_ /; a > 1}]] &,
   Flatten[{
     x,
     Transpose@x,
     Flatten[Table[Flatten[Take[x, 3 i - {2, 0}, 3 j - {2, 0}]], {i, 1, 2}, {j, 1, 2}], 1]
     }, 1]]

再帰の深さ制限をなくしてから探索すれば、どんな問題でも解けますが、さすがに時間がかかりすぎます(「scanOr」を「Or @@ Map」に置き換えればすべての解を見つけますが、やめたほうがいいでしょう)。

$RecursionLimit = Infinity;

Reap@search@{
  {1,0,0,0,0,7,0,9,0},
  {0,3,0,0,2,0,0,0,8},
  {0,0,9,6,0,0,5,0,0},
  {0,0,5,3,0,0,9,0,0},
  {0,1,0,0,8,0,0,0,2},
  {6,0,0,0,0,4,0,0,0},
  {3,0,0,0,0,0,0,1,0},
  {0,4,0,0,0,0,0,0,7},
  {0,0,7,0,0,0,3,0,0}}

方法2(これも力任せ)

Ruby本のように、入りうる数字(candidates)をルール通りに決めて、それだけを調べるようにしてみましょう。

述語deepenは、空白の場所(pos)を、そこに入りうる数字(candidates)で置き換え(ReplacePart)、探索を進めます(入りうる数字が無いと、Or @@ Mapの結果はFalseになります)。入りうる数字だけを使うので、ルールのチェック(test)は不要です。

deepen[x_] := With[{pos = First@Position[x, 0]}, 
  Or @@ Map[search, (ReplacePart[x, pos -> #] & /@ candidates[x, pos])]]

i行j列に入りうる数字(candidates)は、1から9の数字(Range[1,9])から、i行目(board[[i]])とj列目(board[[Range[1, 9], j]])、そのマス目が属するブロック(Flatten…)に属する数字を除いたもの(Complement)です。

candidates[board_, {i_, j_}] := Complement[Range[1, 9],
  board[[i]],
  board[[Range[1, 9], j]],
  Flatten[Take[board, 3 Ceiling[i/3] - {2, 0}, 3 Ceiling[j/3] - {2, 0}]]]

今度はすぐに解が求まります(「Or @@ Map」なので、すべての解を求めます)。

Reap@search@{
  {1,0,0,0,0,7,0,9,0},
  {0,3,0,0,2,0,0,0,8},
  {0,0,9,6,0,0,5,0,0},
  {0,0,5,3,0,0,9,0,0},
  {0,1,0,0,8,0,0,0,2},
  {6,0,0,0,0,4,0,0,0},
  {3,0,0,0,0,0,0,1,0},
  {0,4,0,0,0,0,0,0,7},
  {0,0,7,0,0,0,3,0,0}}

{True, {{1 6 2 8 5 7 4 9 3}}}
         5 3 4 1 2 9 6 7 8
         7 8 9 6 4 3 5 2 1
         4 7 5 3 1 2 9 8 6
         9 1 3 5 8 6 7 4 2
         6 2 8 7 9 4 1 3 5
         3 5 6 4 7 8 2 1 9
         2 4 1 9 3 5 8 6 7
         8 9 7 2 6 1 3 5 4

結論

空行を取り除くと、ちょうど10行のコードが残ります。これは、単純な力任せのアルゴリズムに頼る関数型言語の数独ソルバとしてはまずまずの長さです。このサンプルは、Mathematicaのパワーと表現力をよく示していると思いますがどうでしょう。

search[x_] := Or[goal@x, deepen@x]
goal[x_] := And[Not@MemberQ[x, 0, 2], report@x]
report[x_] := (Sow@TableForm@x; True)
deepen[x_] := With[{pos = First@Position[x, 0]}, 
  Or @@ Map[search, (ReplacePart[x, pos -> #] & /@ candidates[x, pos])]]
candidates[board_, {i_, j_}] := Complement[Range[1, 9],
  board[[i]],
  board[[Range[1, 9], j]],
  Flatten[Take[board, 3 Ceiling[i/3] - {2, 0}, 3 Ceiling[j/3] - {2, 0}]]]
$RecursionLimit = Infinity;

補足:Mathematica 5で試す場合には、述語deepenの中のReplacePart[x, pos -> #]をReplacePart[x, #, pos]に置き換えてください。