市役所からのお知らせです。

＿人人人人人人人人人人人人人人人人人人人人人人人人人人＿
＞　この記事では IEEE754 偶数丸めを仮定しています！　＜
￣ＹＹＹＹＹＹＹＹＹＹＹＹＹＹＹＹＹＹＹＹＹＹＹＹＹＹ￣

述語によっては丸め方向が大事なものがあるのですが、いちいちお断りするの大変ですからね。世は大省エネ時代です。

こんなみなさまに読んでいただきたい！

ところでみなさま！初めて算術回路を学んだときの感動を憶えていらっしゃいますか？

みなさま「足し算なんて私の知っている筆算を回路に落とすだけでしょう、私は天才！人生楽勝！」

みなさま「……高さ (= 回路の最長路長) O(log w) ！？」（w はワード長）

ご存知なかった方はこの機会にそちらもチェックしていただきつつ、まずは私の一勝ということでね。

そういえば算術回路についての名著『アルゴリズムイントロダクション』の第2版第3巻（今手に入るのですかねこれ）に書いてありました。第3版で廃止されたのですが、時代ですね🙃

ともあれこういうのまだまだやりたいではありませんか！浮動小数点数 & 計算幾何 𝐼𝑠 𝐴𝑙𝑙 𝑌𝑜𝑢 𝑁𝑒𝑒𝑑（文法ミス）です。やっていきましょう！

Section 1: 適応的精度 (Adaptive Precision) の概要

適応的精度 (Adaptive Precision) は、状況に応じて計算精度を徐々に上げていくことで、最小限のコストで厳密な不等式判定を可能にする手法で、丸め誤差をうまく再利用するのがポイントです。

Adaptive Precision の詳細については Shewchuk (1997) に詳しいです。みなさまもぜひです。

適応的精度（Adaptive Precision）を知ったきっかけ

3D CAD といえば B-Rep (Boundary Representation) です。要はトポロジー（面などの接続関係）と、ジオメトリー（それらの方程式）を別で管理するということですね。

これがまた大変です。決まった形をしていませんから、機械学習がとても大変です。（実はここ最近 B-Rep 向けの機械学習もアツいのですがそれはまた別のお話🖐️）

というわけでみんな大好きルールベースになるわけですが、これアルゴリズムをある程度頑健にしておかないとデバッグが難しいのですよね。そのへんでいろいろ調べていたら知ってしまった†深淵†のひとつです。私の右目が疼く前に次のセクションへ GO! です。

適応的精度（Adaptive Precision）の動機は？

計算幾何あるあるをご紹介しましょう。

たとえば、2 次元の与えられた 3 点 A, B, C の配置がこの順に反時計回り、時計回り、共線（同一直線上）のどれであるかを判定する問題を考えてみましょう。

みなさま「こういうのは行列式を使えばよいのですよね。数学に詳しすぎて照れてしまいますね」

みなさま「あれ、こ〜れ境界値周りの挙動が一貫していなかったらぶっ壊れますでございますわ〜☀️☀️☀️」

共感していただけましたでしょうかね。要は別に微妙な結果がどちらに転ぶのかをすごく気にしている訳ではなくて、大事なのはズレるのはよいとしてどちらに転ぶかが一貫していなかったら全体がぐちゃぐちゃになってしまうということです。

反時計回りであるかどうかなど、Yes/No が定まるものを述語(predicate) といい（例：反時計回りか？円の内側か？など）、これを厳密に判定する手法は厳密幾何計算 (Exact Geometric Computation, EGC) と呼ばれていて、数値そのものを正確に計算する厳密算術計算 (Exact Arithmetic Computation) とは区別されます。符号などの判定だけを厳密にするので、数値全体を厳密化するより軽いのがポイントです。

せっかくですから、Kettner, Lutz, et al. (2008) から印象的な図表をいくつか引用しましょう。

印象的な図表 1

まずは反時計回り・時計回り・共線判定がいかに数値的に不安定であるのかを見ていきましょう。

キャプションにある 3 点 A, B, C を、ほんのちょっとだけ動かした点 A', B', C' が、どちら向きに並んでいるかを浮動小数点数演算で判定した結果です。

反時計回りになっている 👉 青色
時計回りになっている 👉 赤色
共線になっている 👉 黄色

ちなみに黒色の細い対角線は、真に共線になるべきところです。え！？！？

もはや絨毯として販売できそうなカオスっぷりです。

印象的な図表 2

次は、反時計回り・時計回り・共線判定が不安定なせいで、いかに派手に幾何アルゴリズムが \ﾄﾞﾝｶﾞﾗｶﾞｯｼｬｰﾝ💥/ するのかを見ていきましょう！

むかしむかしある座標平面に、点 p₁ さんと、点 p₂ さんと、点 p₃ さんと、点 p₄ さんが、この順に反時計回りの凸四角形を構成しながら暮らしていました。

おばあさんが川で洗濯をしていると、大きな桃🍑が流れてきました。

その桃🍑を切ると、なんと点 p₁ のすぐ右に点 p₅ がおぎゃあ、おぎゃあと………

ってこれいつまでやればよいでしょうかね。止めないと私無限に続けますよ！？！？

そんなことより凸性です！点 p₅ を追加した結果は (a) の通りです。

点 p₁ と点 p₅ の差が丸め誤差に吸収されて、ちょっと凹んでしまっていまね。ウーンム……

みなさま「こんなん誤差やん。うおｗ」

私「1εを笑うものは1εに泣きますよ。」

ここで (b) の拡大図をご覧ください！点 p₆ を追加するとあらふしぎ。前提が壊れているのでぐ〜〜〜〜っちゃぐちゃです。

これでもう「うおｗ」とは言わせませんよ。

いたるところで大活躍する基本的な述語 4 選！

Orient2d

みなさまの記憶力は素晴らしいため、問題を憶えていらっしゃいます。

2 次元において、点 c が直線 ab のどちら側にあるかは次の行列式で決まります。

ところでこれ命名の流派たくさんあって毎回ウーンムとなります。厳密計算の分野における古典的な呼び方はおそらく Orient2d, Orientation2d なのですけれども、たとえば nalgebra では perp (perpendicular product の略？) と呼ばれていたり、私はじめ競プロ経験のある方は ccw (counter-clockwise の略）を使うことが多かったり、単に det と書く方もいらっしゃったりなどです。

Orient3d

3 次元バージョンです。萌え萌えでなくて申し訳ないです。

3 次元において、点 d が平面 abc のどちら側にあるかは次の行列式で決まります。

InCircle

2 次元において、点 d が円 abc の内外どちらにあるかは次の行列式で決まります。

InSphere

3 次元において、点 e が球 abcd の内外どちらにあるかは次の行列式で決まります。

これらの主な用途を挙げておきますね。聞いたことのない方には先着で読み上げ音声をプレゼントです。

アルゴリズム	活躍する述語
凸包構築	Orient2d
線分交差判定・線分群の交差列挙 (Bentley–Ottmann algorithm)	Orient2d
ドロネー三角形分割	Orient2d, InCircle
ボロノイ図の構築 (Fortune’s Algorithm)	InCircle
メッシュの refinement (Delaunay Refinement)	Orient2d, InCircle

3 次元でも基本的には同様の問題があって、対応する述語が活躍することが多いです。

Adaptive Precision のアプローチ

お話がそれましたわ。これらの述語をどういうアプローチで「必要なときだけ厳密に」「お爆速に」行うのかをみていきましょうですわ〜〜〜

大事なことは、近似計算でも結果が 0 から大きく離れていたら信じて OK （= 推定誤差上界より十分大きい）ということです。

ある精度で近似計算を行います。
誤差上界を評価します。
1. 0 から大きく離れている場合 👉 return 結果
2. 0 に近い場合 👉 精度を上げて 1 行目にもどる or どこかの段階で任意精度計算で厳密な判定を行います。

アッ、まだ帰らないでいただいて…… わかります。

みなさま「ハイハイ、どうせ任意精度(= 多倍長)小数でしょ。┐(´д｀)┌」

私「面白いのはここからです！！！（☝ ՞ਊ ՞）☝」

これでうまくいくのは当たり前ではありませんかと、ええ、うんうんそれは彼氏さんが悪いです。

なんと Adaptive Precision の方法は動的メモリ確保なし・ほとんど条件分岐なしで行うことができます。しかも精度を上げて再計算するときには前回の計算の途中結果を使いまわすことができます。その裏にあるのは、丸め誤差をキレイに分類するための、まるで回路図のような洗練されたアルゴリズムです。

たとえば Orient2d は下図のようになります。

アッ、まだ帰らないでいただいて…… 大丈夫です、私もわかりませんでした。

目標まとめ

やりたいのはこちらでしたね。

Orient2d, Orient3d, InCircle, InSphere などの述語を厳密に判定したいです。
境界値に近くない場合は計算が早く終わっていただきたいです。
動的メモリ確保なしです。（従って普通の任意精度小数は使えません！）
条件分岐は必要最低限です。（従って算術計算なんかで場合分けは使えません！）

全部は大変ですから、みなさまがご興味を持っていただけそうなところまでご招待したいと思います。

Section 2: 足し算・掛け算の実装

2 つの浮動小数点数の和や差を、どうやって表しましょうか

2 つの浮動小数点数 a, b の和 a + b や積 ab の正確な値は、1 つの浮動小数点数では表現できません。

しかし！ここで軽率に普通の任意精度小数を使ってしまうと、たとえば 2¹⁰²³ と 2⁻¹⁰²² の和などはと〜〜〜んでもない長さになってしまいますね。長いものがお好きな方以外は悲しくなってしまいます。

そこで、和による表現です。a + b や ab を 2 つの浮動小数点数 x, y の和で表してしまいましょう。といっても何でもアリにしてしまうと後で困るので、上手い条件を課します。a + b や ab の場合は実は x が丸めた値、y が丸め誤差に一致するようにするのがツウなのですが、一旦もう少し一般のお話をいたしましょうか。

任意精度小数の一般的な表現

有限小数 x に対し、それを浮動小数点数 x₁, x₂, ..., xₙ の総和 x = x₁ + x₂ + ⋯ + xₙ で表したものを x の表現 (Expression) と呼ぶことにしましょう。

そしてしばしば次のような条件を課します。

nonoverlapping

表現 x = x₁ + x₂ + ⋯ + xₙ が nonoverlapping であるという用語を定義していきましょう。

まずは準備です。浮動小数点数 a, b が次のいずれかの条件を満たすとき、a と b は nonoverlapping であるといいます:

a の最下位ビットが、b の最上位ビットよりも真に大きい
b の最下位ビットが、a の最上位ビットよりも真に大きい
a = 0
b = 0

表現 x = x₁ + x₂ + ⋯ + xₙ は、そのどの相異なる 2 項も nonoverlapping であるときに nonoverlapping であるといいます。

「ビットの範囲が重なっていない」と言い換えるとわかりやすいかもでしょうかね。

nonadjacent

次は、表現 x = x₁ + x₂ + ⋯ + xₙ が nonadjacent であるという用語を定義していきましょう。

まずは準備です。浮動小数点数 a, b が次のいずれかの条件を満たすとき、a と b は nonadjacent であるといいます

a の最下位ビットが、2b の最上位ビットよりも真に大きい
b の最下位ビットが、2a の最上位ビットよりも真に大きい
a = 0
b = 0

表現 x = x₁ + x₂ + ⋯ + xₙ は、そのどの相異なる 2 項も nonadjacent であるときに nonadjacent であるといいます。

「ビットの範囲の間に 1 bit 以上隙間が空いている」と言い換えるとわかりやすいかもでしょうかね。

みなさま「わざわざ隙間を 1 開ける理由はなんですか？」

わかります。浮動小数点数演算 a + b や ab の結果を丸めた値 x と丸め誤差 y の関係を思い出してみましょう。計算機イプシロン ε を用いて |y| ≤ (x の最上位ビット) · ε / 2 が成り立ちますね？はい、一つ隙間が空いているため 𝑣𝑖𝑐𝑡𝑜𝑟𝑦 です。逆にこれ以上開けるのを保証するのは不可能です。

誤差をバッチリ管理して算術演算をする方法

記法

浮動小数点数 a, b に対して、a + b, a - b, ab を丸めた値をそれぞれ a ⊕ b, a ⊖ b, a ⊗ b と表します。

浮動小数点数の和 a + b

簡単のため |a| ≥ |b| を仮定しましょう。すると筆算の様子はだいたいこういう感じですね。実際には a, b, x, y のどれか負になったりもしますけれども、そのあたりはご愛嬌です。

x は浮動小数点数演算でそのままGO!です。y は真の値 x の差 (a + b) - x を丁寧に計算すれば良いですね。具体的には y = b - (x - a) と表すと途中どこにも丸め誤差が生じなくて 𝑤𝑖𝑛 です。👇 疑似コード

みなさま「え、あの、 |a| ≥ |b| ってつまりそれを保証するために条件分岐が入りますよね？？？」

ウーンム、痛いところを突かれましたね。しかし実は y = (a - (x - (x - a))) + (b - (x - a)) がそれを解決するのですよ。証明は難しくはありませんが場合わけが多いですから、みなさまのご担当ということでよろしくお願いさせていただいてですね。そのかわり明日の給食当番は私にお任せください🖐️

Expression e と浮動小数点数 b の和

結果となる expression は長さが 1 長くなります。この世界では演算をしても項の個数が減らないという地獄みたいな設定ですからね🫠

e = e₁ + e₂ + ⋯ + eₙ が nonadjacent かつ、「一部に 0 があるかもしれないことを除いて絶対値が単調増加」という前提条件を課して、次の回路で計算すると、実はその性質が保たれます。

Expression 同士の和 e + f

ドーーン！Go Bold の精神です。

浮動小数点数同士の積 ab

これまた賢いですよね。実は入力 a, b をともに上位ビットと下位ビットに分解するとうまくいきます。というかよく考えたらそうでもしないと精度が足りないのは明らかですよね。（p 桁の正整数と q 桁の正整数の積が常に p + q - 1 桁であるという事実に注意です。）

それに注意した上で、足し算のとき同様に、丸め誤差が出ないように少しずつ差を計算していく方針です。疑似コードを追っていただいてもわかるかもしれませんし、難しければぜひとも原典をです。

Expression e と浮動小数点数 b の積

次のような回路で計算できて、項の個数は 2 倍になりますね。

とはいえ論文中のどの述語の実装でも掛け算は浮動小数点数同士でしか行っていませんから、今回必要はないみたいですね。

これで一通り回路のパーツが揃いました！あとは各述語 Orient2d, Orient3d, InCircle, InSphere を実装していくだけです！

……というのがまあ地獄みたいな作業なのですよね。

あと念の為注意なのですが、出力される expression の項は 0 を含むことがあります。例えば a + b を浮動小数点数で計算しても全く誤差の出ない場合は y = 0 になりますよね。これをハンドルしようと思うと当然場合分けが出てくるため、意図的に放置して、固定演算だけで済むようにしております。賢いですね。

Section 3: 各種幾何述語の adaptive precision 実装（の方針）

目標の下方修正

Orient2d のみを実装しましょう。

えいやだっていやじゃないですか。—— みつを

Orient2d の展開

定義を展開しましょう。

これをそのまま回路に起こした結果がこちらです。

アッ、まだ帰らないでいただいて…… いまから理解できますから、たぶん、4 割くらいは。（あの？）

最初の引き算

まずは入力を引き算して、主要項と誤差項にわかれますね。このとき、主要項はオーダー O(1)を、誤差項はオーダー O(ε) を持つとみなしましょう。実際には誤差項は主要項のε倍よりもずっと小さいこともありますが、小さい分にはよいのでこのまま形式的に突っ走っていきましょう！

次の掛け算

さらに図の左半分と右半分それぞれについて {主要項, 誤差項} と {主要項, 誤差項} の 4 種類の掛け算をして、8 つずつ、16 個の項を導出しましょう。これは図の左から順にオーダー O(1), O(ε), O(ε), O(ε²), O(ε), O(ε²), O(ε²), O(ε³) （左半分だけ記述）を持ちます。あとはこれらを全て足し引きしたものが厳密な答えです。

厳密な答え D

実際に Expression の足し引きを行なったものが D です。いかにこれを計算しないかが大事ですね。以降、さまざまな近似をみていきましょう！

最も粗い近似 A

浮動小数点数でそのまま計算した結果に相当します。

A の誤差範囲

慣れている方ならおかわり🍚いただけると思いますけれども、Orient2d 唯一の危険ポイントは減算 x₅ − x₆ による桁落ちです。そこで |x₅| + |x₆| の定数倍により評価していきましょう。

こちらが温めておいた誤差限界です🍳（料理番組風）

|A| がこの値以上でしたら、sign(A) = sign(D) が保証されます。

え、私は証明を追ったのかですって？

さて、次は A と D の中間の精度の値です。

第2の近似: B'

これは正直なところ A とほとんど変わらないので、私はあまり存在意義がわからなかったのですが、なにかご存知の方いらしたら教えていただきたいです。

B’ というのは Expression B の浮動小数点数近似です。

みなさま「そんなの B の O(1) 項を取れば良いだけでしょうそんな、またまたぁ〜〜」

わかります。私もそう勘違いしたせいで無限にお時間が溶けました。ここで扱っている expression というのは、主要ビットから貪欲になるべく長くとったものとは限らないということに注意です。ここで躓かなかった方は 𝑌𝑜𝑢 𝑘𝑛𝑜𝑤 𝑡𝑜𝑜 𝑚𝑢𝑐ℎ、バッドエンドです。

実は論文中（Theorem 23）に、主要項を抽出するのに使えるアルゴリズムが記載されています。詳細はさておき、条件分岐が出てきているのはちょっと悲しいところですね。（ちなみにループは今回固定長なので消せます。）

誤差上界をはこちらです。A とあまり変わらなくて悲しいですね。

第3の近似: C

O(ε) のオーダーを持つ項を片っ端から集めていることがわかりますね。

誤差限界はこちらです:

これでもどうしようもなければ D を計算するということですね。賢いです。

Section 4: 現代の厳密計算

結局私個人が厳密幾何計算を実装する（ハメになる）機会はなかったのですが、では現代ではもう使われていない技術なのでしょうかね。だとしたら、とても寂しいですね。

最後に、この手の厳密述語 (exact predicate) が現代のライブラリでどう扱われているかを軽く見ます。

CGAL にありましたよ！！こちら、どなたが使っているのですか！？

CGAL (Computational Geometry Algorithms Library) というライブラリがあります。計算幾何のあらゆる叡智の詰まった老舗の C++ ライブラリなのですが、そこを探してみると、Predefined Kernels というクラス群が見つかります。これは幾何計算をしたり、幾何的オブジェクトを構築したりするときのカーネルです。

命名規則を見ると、exact predicate と exact construction が見えますね。Exact predicate というのは今回ご説明したように、「述語の yes/no を厳密に判定する」ことで、比較的簡単かつあまり大きなコストがかからないと言われています。（簡単ですか？？？？？？）

一方 exact construction は非常に難しいです。というのも predicate のときのように「誤差を評価してフィルタリング」のようなテクは通用せず、すべてを厳密計算する必要があります。そのうえ複雑な制約条件により与えられる点は長大な有理数や複雑な方程式で表されるからです。「構築した点を基準に構築する」とかにも対応しないといけませんからね。

OCCT にはなさそうですね？？？

CAD の超有名 OSS であるところの OCCT (Open Cascade Technology) もチェックです！一応探してみましたが、なさそうな気がします。

根拠列挙ターイム！

OCCT は他の CAD カーネル同様、tolerance ベースで管理を行っていますね。
CAD カーネルは点の構築が大前提のソフトですから、exact predicate もうまく活躍できなそうです。
Tolerance によるヒーリングの技術も発達しているので、それで実は実用上困らないのかもしれません。
実際、非常に繊細で複雑な計算である Boolean 演算でも tolerance ベースで計算していることが、Boolean Operations の Guide を参照するとわかります。

BTW: 機械学習ってなんぼのもんなんですかね

さて、ここまで読んでくださったアルゴ系エンジニアのみなさま、機械学習はお好きですか？

みなさま「機械学習がなんぼのもんじゃい✊💥」

みなさま「私たちの憧れたヒューリ🥒スティックプログラミングいずくにかです！」

え、今日日そんなこと思いませんよって、ぐぬぬぅ……

というか私が☝️でしたというだけなのですけれども、最近は変わりました。

機械学習は、やや矮小化かもしれませんけれども「書かなくていい部分が増えただけ」で、古典的手法のエッセンスは残っている気がするのですよね、少なくとも現時点では。

良いヒューリスティックは大抵「天才的な気づき」から始まるものの、その後はつらいモグラ叩きになりがちです。機械学習は、そのモグラ叩きを「最適化問題」として正面から殴れるのが最高に気持ち良いです！

Transformerのような汎用手法が全盛の今でも、何が本質かを見極めてベクトルを練り上げる工程は、極めて domain-specific です。そこには間違いなく、古典の英知が息づいています。

現実世界には「理論上解けるはずなのに、放置されている面白い最適化問題」が山積みです。その山を、古典の英知と最新の技術で一緒に崩しませんか？

そんな仲間を、我々キャディのチームは待っています😉

speakerdeck.com

まとめ

何かをするのに遅すぎるということはありません。

おととい 4/11 は私のお誕生日でした。プレゼント🎁を送れなかったよという方、あとはわかりますね。