English edition is here.

東京は芝で開催された、あるディジタル信号処理研究会での近況報告で、私はおずおずと話し始めた。

「惚気のような話かも知れなくて恐縮ですが、私の妻はローストビーフが好きで…」

私の妻はローストビーフが好きで―

「イギリス料理に美味しいものはほとんど無いけど、ローストビーフだけは評価する！」

「そこは、古き良きヒューレット・パッカードのモットー 『推察するな！測定しろ！』 だよ」

「…と、いうことで、調理用の温度計買っとくねー！」

そして、A****n で安く売られていた調理用の温度計が届いた。 試してみるとプローブの先では上手く温度が測れずに 30 mm ほど入れないとダメみたいだし応答がとても鈍い。 分解してみると、プローブのステンレス・シース内のセンサーの配線が短く、プローブ先端に届いても触れてもいないみたいだ。 足りない長さが、ちょうどシース・キャップの長さの 20 mm 程なので、きっと製造指示書の間違いだったのに違いない。

「これが安さの秘密だったみたいです」

そそくさと、センサーの配線を継ぎ足して組み立て直し、改めて温度を測ってみると思いの外正確っぽい。 氷水と、沸騰する精製水では、正確に 0.0 ℃ 及び 100.0 ℃ を示し、校正済みの K 熱電対との比較でも互いに 20.0 ℃ であったりと正確だった。 改めて調べてみると 300 ℃ 対応の DO-35 型のガラス封入サーミスタが使われてた。

「チップ型のはよく使ってますが、DO-35…1N4148 みたいなアキシャル・リードのは、見過ごしてて使ったことがありませんでした」

（sqrt と cbrt を間違えたまま、10 年以上ずっと使っていたのか…）

「組み込み用に sqrt はあっても cbrt は無いっすよ！」

cbrt は 1/3 乗なので、それの近似関数を sqrt で作って √√x√√x√x (= x^11/32 ) としたら、レンジリダクション無しで、1e±40 の範囲で誤差 -0.7〜1.6 倍程度の近似ができるのが分かった。 それに対して、解の近傍での (1, 1) Padé 近似を漸化式に適用したら 3, 4回の繰り返しで double のほぼ限界まで収束するのが分かった。 出来た cbrt の関数はそのままで複素数対応だった。 複素数対応には近似の展開点だけでなく特異点で漸近的に等しくなるように近似すると上手くいくことが多い。 でも、double の範囲で使うにはちょっと無駄な気がする

寝ている時に思いついた、ピュイズー級数的かつ制限付きのチェビシェフ近似的な実数版の cbrt 用の初期近似を、目が覚めてから試したら、 x ∈ [0.125, 8] の範囲で 3 項で 1.83 桁、4 項で 2.31 桁の精度が得られるのが分かった。 これに先の (1,1) Pade − 意味的には Hally 法と同じ −の漸化式 2 回で double の精度になるはず。

他に何かヒントはあるかなと Gnu C library (version 2.41) の cbrt のソース・コードを見てみたら、予想に反してなんだか詰めが甘い気がする。 1,000,000 個のランダムな倍精度の値で評価してみると、IEEE/dbl-64 版だと誤差が +9/-7 ulp と随分精度が悪い。 速度重視で精度はあまり気にしてないのか？ gcc の native 版でも +3/-2.5 ulp で改善の余地がありそう。

一方、私の 4 項の初期近似のコードで試したら倍精度でも拡張倍精度でも ±2 ulp の精度で数値計算の限界に到達出来てた。 コードもシンプルで移植性も高い。

「これで組み込めるよね！」

「ところで先程は修正って一言で済ませましたが、５桁ほどのダイナミックレンジの非線形で誤差を含む限られたサンプル点から滑らかな近似曲線を推定するのは…」

「それで美味しいローストビーフは食べられたんですか？」

「あ、それは

『美味しいローストビーフのためには、温度計じゃなくて美味しいお肉を買うのが正しいのよ！』

って言われました！次は遺伝子工学とかですかね？」

「それ違うー！」
「可愛い奥さんに美味しいお肉買ってあげてー！」

テクニカルノート

• 3係数による NTC サーミスタの近似式
• 平方根を用いた立方根の近似と複素数および拡張倍精度実数での実装例