(インテルとAMDのCPUの)三角関数は不要


インテルとAMDのCPUには,三角関数を数値的に計算するための命令が備えられていますが,その結果は正しくありません。たとえば,CPUの命令FCOSの結果とstd::cosでソフトウェア的に計算した結果を比べると,後者の方が真の値に近くなります。

再現のためのコード(FCOSを使う部分だけインラインアセンブラを使っています。)

Core i5 6600とg++ 4.8.4で試した結果はこんな感じです。

x = 0.98233490762409514385
c1 = fcos(x)     = 0.55508189550073283591
c2 = std::cos(x) = 0.55508189550073294694
c  = mp::cos(x)  = 0.555081895500732891425063460345544265145531916268
abs(c - c1) = 5.5512e-17
abs(c - c2) = 5.551e-17

c1 is better: 0
c2 is better: 25
total: 100000

(その正しさをどうやって確かめるのかという問題はありますが)Boost.Multiprecisionで計算した結果と比較すると,命令`FCOS`より`std::cos`の方が正確なようです(Visual Studioの場合はまた違った結果になるようです)。

Pentiumが割り算を間違うといって大騒ぎになったことがありましたが(Pentium FDIV バグ),数値計算はどうせ近似値だからでしょうか,この問題は,ドキュメントに注意書きが追記されただけで,根本的な解決はされないままのようです。

参考:サイン、コサインをインテルの CPU で計算すると少しバグっているらしい

MathematicaのMaxValueとMinValueのバグ(11.2で解決)


11.2で解決

Mathematica 10.1, 10.2, 10.3, 10.4.1, 11.1.1

Mathematicaの最大化MaxValueと最小化MinValueには,簡約のためのSimplifyの中では使えないというバグがあります(製造元には報告済みです)。

例として,0 <= t <= 1, 0 <= p <= 2 Pi, 0 <= q <= 2 Piという制約のもとで,f = Abs[t Sin[p] Sin[q]]という関数の最大値と最小値を求めます。

f = Abs[t Sin[p] Sin[q]];
cond = And[0 <= t <= 1, 0 <= p <= 2 Pi, 0 <= q <= 2 Pi];

当たりを付けるために,まずは数値的に求めます。

{NMaxValue[{f, cond}, {t, p, q}], NMinValue[{f, cond}, {t, p, q}]}

結果が{1., 0.}になることに,特に問題はないでしょう。

解析的に求めようとすると,うまく行きません。(11.2ではうまく行きます。)

MaxValue[{f, cond}, {t, p, q}]
(*結果は割愛。入力がそのまま返される*)

こういう関数の最大・最小は,そのままではうまくいかないとしたものです。残念ですが,これはしょうがない。

しかし,Simplifyの中で使うと計算が進むことがあります。制約条件を仮定して結果を整理することを試みます。

Simplify[MaxValue[{f, cond}, {t, p, q}], cond]

得られる結果は「∞」,もちろん間違いです。MinValueの場合も同様で,「-∞」という間違った結果が得られます。

Simplifyの仮定(外側のcond)が,MaxValueの制約(内側のcond)を簡約しているのが一因だと思います。そういうことがあるのは,次の例でわかります。

Simplify[MaxValue[{1/x^2, 0 < x}, x, Integers], 0 < x]

結果はMaxValue[{x^(-2), True}, x, Integers]になります(11.2では正しい結果「1」が得られます)。制約がTrueになっていますが,これはいけません。

Tweet-a-Programではちょっと違うことが起きているみたいです(こういうのは初めて見ました)。

他に原因があるかどうかはよくわかりません。

MaximizeMinimaizeでは,この問題は発生しません。一般に,MaximizeMinimizeMaxValueMinValueより遅いということになっているのですが,速さよりは正確さが大事なので,こちらを使った方がいいのかもしれません。

MathematicaのIntegrateのバグ


Mathematica 10.1, 10.2, 10.3, 10.4.1, 11.2の積分(Integrate)には,仮定の利用に関して,9.0や10.0には無かったバグがあります。(製造元には報告済みです。)

例として,Abs[(x + 1) (x - a) (x - 1)]-1 <= x <= 1で積分することを考えます。

Absの中にパラメータがaがありますが,これは実数だと仮定します。Integrateには,そういう仮定を与えるための便利なオプションAssumptionsがあります。

Integrate[Abs[(x + 1) (x - a) (x - 1)], {x, -1, 1}, 
 Assumptions -> Element[a, Reals]]

これで,a <= -1なら-4a/3a >= 1なら4a/3,それ以外つまり-1 < a < 1なら(-a^4 + 6a^2 + 3)/6と,aの値で場合分けされた結果が得られます。すばらしい。

しかし,仮定がうまく取り入れられない,取り入れられないならまだしも,仮定を使ったせいで間違った答えが出てきてしまうことがあります。

Integrate[Abs[(x + 1) (x - 2 a + b) (x - 1)], {x, -1, 1}, 
 Assumptions -> And[-1 < a < 1, b == a]]

Mathematicaの一部のバージョンでは1/2という結果になりますが,これは間違いです。b == aという仮定を考慮すると被積分関数は最初の例と同じになるので,積分結果は(-a^4 + 6a^2 + 3)/6にならなければなりません。

どういうわけか,Integrateの中に直接数式を書かなければうまくいきます。

f = Abs[(x + 1) (x - a) (x - 1)];
Integrate[f, {x, -1, 1}, 
 Assumptions -> And[-1 < a < 1, b == a]]

原因はよくわかりませんが,Assumptionsは怖くて使えません。

この積分は,Mathematica 9.0(Windows)や10.0(Raspberry Pi)では正しく行えただけに,10.1でできなくなって残念です。機能追加と品質保持の優先順位が間違っている気がします。

Mathematica の新バージョンには,常に多くの新しい機能が含まれている.しかし当初からの周到なデザインにより,すべてのバージョン間でほぼ完全な互換性が保たれている.その結果,例えば1988年のバージョン1用に書かれたほとんどすべてのプログラムは,Mathematica バージョン7でもそのまま変更なしで走り,かつ実行速度も大いに向上している.(Mathematica バージョン1以降の非互換変更

この理想を,バージョン8以降でも大切にしてほしいものです。

MathematicaのSolveとReduceのバグ(10.2で修正)


10.2で修正されました。

Mathematica 10, 10.1のSolveとReduceには,ベクトルや行列が等しくないという条件を正しく扱えないという,かなり深刻なバグがあります。(製造元には報告済みです。)

例として,
x^2+y^2 == 1かつy == 1/Sqrt[2]
という方程式を
{x, y} != {-1/Sqrt[2], 1/Sqrt[2]}
という条件のもとで解きます。

Solve[And[
  x^2 + y^2 == 1,
  y == 1/Sqrt[2],
  {x, y} != {-1/Sqrt[2], 1/Sqrt[2]}],
 {x, y}]

Mathematica 9.0.1では
{x -> 1/Sqrt[2], y -> 1/Sqrt[2]}
という正しい解が得られます。

Mathematica 10.1では
解なし
という間違った結果になります。

いろいろ試してみると,Solveの中に書いた
{x, y} != {-1/Sqrt[2], 1/Sqrt[2]}
という条件が,
x != -1/Sqrt[2] || y != 1/Sqrt[2]
ではなく,
x != -1/Sqrt[2] && y != 1/Sqrt[2]
と解釈されているようです。

Reduceでも同じ問題が起こります。

Reduce[{x, y} != {0, 0}, {x, y}]

Mathematica 9.0では
x!=0 || y!=0
という正しい解が得られます。

Mathematica 10.1では
x!=0 && y!=0
という間違った解が得られます。

というわけで,Mathematica 10.1では,SolveやReduceの中に,ベクトルや行列が等しくないという条件は書けません。とりあえずの解決策は,Or[x != 0, y != 0]のように成分ごとに書くことですが,これでは2次元限定ですし,そもそもMathematicaで成分を書いたら負けな気がするので,「例えば,Mathematica 10.1を避ける」でも仕方ないでしょう。

製造元のウェブサイトには,「世界で最も信頼できる最新技術計算システム」とありますが,それはないと思います。

MathematicaのTeXFormのバグ


Mathematica 10.1, 10.2, 10.3, 10.4.1, 11.2

Mathematicaには,表現をTeX形式に変換するTeXFormがあります。たとえば,TeXForm[Integrate[f[x], x]]とすると「\int f(x) \, dx」を返してくれます。便利です。

しかし,バグがあります。

TeXForm[Abs[r + 1/2]]の結果が「\left\left| r+\frac{1}{2}\right\right|」になりますが,これはTeXで処理できません。正しくは「\left| r+\frac{1}{2}\right|」です。TeXForm[Abs[r]]なら「\left| r\right|」という正しい結果が返ります。

開発元には報告済みですが,今のところ手で直すしかありません。計算結果から自動的に論文を書かせるようなどと思うと困ったことになるでしょう。