区間代数と無限小と無限大

区間代数で無限大と無限小を統一的に扱う
======================================

このドキュメントのコードのテストに必要な拡張とライブラリは以下の通りです。

1. numeric-prelude
2. QuickCheck

```haskell
{-# LANGUAGE TypeOperators        #-}
{-# LANGUAGE TypeFamilies         #-}
{-# LANGUAGE DataKinds            #-}
{-# LANGUAGE UndecidableInstances #-}
{-# LANGUAGE GeneralizedNewtypeDeriving #-}
{-# LANGUAGE RebindableSyntax     #-}
{-# LANGUAGE NoImplicitPrelude    #-}
import           Test.QuickCheck
import           NumericPrelude
import qualified Algebra.Field      as Field
import qualified Algebra.Ring       as Ring
import qualified Algebra.Additive   as Additive
import qualified Algebra.ZeroTestable as ZeroTestable
import qualified Number.Ratio       as Ratio
import qualified MathObj.Polynomial as Polynomial
import qualified MathObj.LaurentPolynomial as Laurent
import           Data.List
```

何をしたいか
------------

どんな言語でも区間の計算で必ず現れる面倒に、閉区間と開区間の区別があります。ノートでは簡単に(a,b),(a,b\],\[a,∞)などと書けても、これほど柔軟な文法をサポートする実用プログラミング言語はなさそうです。[^1]

現状、区間代数を扱うHaskellパッケージはいくつかありますが、ほぼ全て次のような直和型を使います。

```haskell
data Interval a = Open   a a
                | Closed a a
                | OpenClosed a a
                | ClosedOpen a a
```

でもこれだと`case`があちこちに必要になって、コードも重複しがちだしめんどくさい、そもそも無限区間の扱い考えたらパターンがまた膨れあがって、憂鬱だ・・・というのが本記事のモチベーションです。

パフォーマンスを検討するために、GHC-8.2.1で導入された`-XUnboxedSums`拡張も調べてみました。

無限小と無限大
--------------

"無限小"(infinitesimal)`d`というシンボルを用意して,

```haskell
open   = (0   , 1  )  -- 開区間
closed = (0-d , 1+d)  -- 閉区間
half   = (-1/d, 1  )  -- 無限区間 (-∞, 1)
```

のように書けたら嬉しくないでしょうか？これならつまらないパターンマッチを省けますが、いくつかの疑問も浮かびます。

1.  ポリモーフィックに定義できるのか？
    どんな型なら無限小dを加えてももともとの型の意味を変えず拡張できるのか？

2.  無限小や無限大を四則演算することの意味は？
    例えば`d+d`や`1/d^2`は何を意味するのか(意味すべきか)？

3.  もとの型に全順序がある場合、
    dで拡張された型にも、全順序は拡張できるのか？

モデル
------

公理を満たす構造をモデルといいます。任意の全順序集合は、全順序の公理を満たすモデルです。Haskellでは`Ord`クラスとして順序比較のインターフェースが用意されています。 QuickCheckで公理を書きくだすと：

```haskell
prop_reflexibity :: Ord a
                 => a -> Bool
prop_reflexibity a = a <= a

prop_antisymmetry :: Ord a
                   => a -> a -> Property
prop_antisymmetry a b = a <= b && b <= a
                   ==> a == b

prop_transitivity :: (Ord a)
                     => a -> a -> a -> Property
prop_transitivity a b c = a <= b && b <= c
                          ==> a <= c
-- totallityはOrdインスタンスで保証される
```

のようになります。これらの公理を満たす`Ord a => a`が全順序のモデルです。

∞を加えてみる
-------------

公理が満たされることを強調するため、以下`Ord`インスタンスではなく`Ord`モデルと呼びましょう。任意のOrdモデルに、「最大値」「最小値」を加えて、全順序を拡張することはできるでしょうか。

```haskell
data BoundedOrd a = Min | X a | Max deriving (Eq, Ord)

instance Bounded  (BoundedOrd a) where
  minBound = Min
  maxBound = Max
```

これで任意の全順序の拡張ができました。この拡張は何度でも繰り返すことができます\[^fn3]。

```haskell
-- Here Nat is promoted to a kind by -XDataKinds
data Nat = Z | S Nat
type family   Power (f :: * -> *) (n :: Nat) (a :: *) :: *
type instance Power f Z a     = a
type instance Power f (S n) a = f (Power f n a)

-- Min < X Min < X (X Min) < .. < X (..(X a) ..) < .. < X Max < Max
type BoundedOrd_n n a = Power BoundedOrd n a
```

このようなモデルの拡張は、もっと複雑な公理系でも可能でしょうか？区間代数で使いたい`(Ord a, Num a) => a`を`Ord`や`Num`の公理を壊さないよう無限大を加えて拡張するには、どうしたらよいでしょうか。

順序体
------

区間代数では、四則演算と順序の関係が大事です。[順序体](https://en.wikipedia.org/wiki/Ordered_field)`OrdField`クラスを定義してみましょう。
numeric-preludeの`Algebra.Field.C`クラスを使うと便利です。

```haskell
class (Ord a, ZeroTestable.C a, Field.C a)
    => OrdField a where
newtype  OF a = OF a
   deriving (Eq, Ord,
             Field.C,
             Ring.C,
             Additive.C,
             ZeroTestable.C)
instance (Ord a, Field.C a, ZeroTestable.C a)
        => OrdField (OF a) where

prop_OrderedAddition :: OrdField a
                      => a -> a -> a -> Property
prop_OrderedAddition a b c =
   a <= b     ==>   a + c <= b + c

prop_LawOfSign :: OrdField a
                  => a -> a -> a -> Property
prop_LawOfSign a b c =
  0 <= a && 0 <= b    ==>   0 <= a * b
```

`OrdField`のモデル`a`に、無限小d(`0 < x ==> 0 < d < x)`を加えたモデルは多項式環 `MathObj.Polynomial a`の商体に、適当な順序を入れたものになっているはずです。適当なライブラリが見つからなかったので、便宜的にクラス`QuotField a`と書きます。

どのような順序か、調べてみましょう。(0/=1な)`OrdField`は整数(1 \< 2=1+1
\< ..)を含みます(つまり標数0)。上の公理から導きだせる、幾つかの性質を見てみましょう。詳しくは[こちら](http://people.reed.edu/~mayer/math112.html/html1/node19.html)を参照してください。

```haskell
-- | 命題1
prop1 :: OrdField a => a -> a -> a -> Property
prop1 a b c =  (a < b) && (c > 0) ==> a*c < b*c
```

`0 < d < 1`から帰納法ですぐに次が言えます。

```
 0 < .. < d^n < .. d^3 < d^2 < d < .. < 1/n < ...  <  1/2  <  1
```

```haskell
-- | 命題2
prop2 :: OrdField a => a -> Property
prop2 a =  a /= 0   ==>   0 < a^2
```

`1/d`は必ず無限大でしょうか？

```haskell
-- | 命題3
-- 0 < 1/(a*b) に注意
prop3 :: OrdField a => a -> a -> Property
prop3 a b =
   sorted [0, a, b] ==> sorted [0, 1/b, 1/a]
```

なので、

```
0 < 1 < 2 .. < n < .. < 1/d < 1/d^2 < 1/d^3 < .. < 1/d^n ...
```

も言えます。

### 定理1

`OrdField`のモデル`a`、a係数の多項式環`Polynomial a`とその商体`QuotField a`に対し、次がなりたつ。

-   `a`の順序を以下のように`QuotField a`に拡張し、`QuotField a`
    を`OrdField`のモデルにできる。

```haskell
instance OrdField a => Ord (Polynomial.T a) where
  x <= y  =  lc (y-x) > 0
    where
      -- the leading coefficient of f
      lc f = last (Polynomial.coeffs f)

newtype QuotField a = QF (Ratio.T (Polynomial.T a))
    deriving(Eq, Additive.C, ZeroTestable.C, Ring.C, Field.C)

instance OrdField a => Ord (QuotField a) where
  QF x <= QF y = let
     (a, b) = (numerator x, denominator x)
     (c, d) = (numerator y, denominator y)
     in  a * d <= b * c
instance OrdField a => OrdField (QuotField a)
```

-   このような拡張はただ一つしかない。

これで任意の`OrdField`のモデルaに、無限小dを加えて`OrdField`の拡張モデルを作ることができました。`a`にシンボル`d`を加えて多項式環を作り、その商体に上のように順序を拡張すればできあがりです。

ローラン環
----------

`Rational`などの順序体でなく、`Integer`などの順序整域`OrdIntegralDomain`を考えるとどうなるでしょうか。numeric-preludeには`Algebra.IntegralDomain.C`という型クラスが用意されています。順序整域をI、その商体をQと書くと、次の環拡大の図式を書けます。

```
  I
  /\
Q  I[d,1/d]
 \  /
 Q(d)
```

ここから上の定理1と同様に、I\[d,1/d\]の順序はただ一つ定まることが言えます(分母を払って比較するものしかない)。 I\[d,1/d\]をローラン環といいます。ローラン環の計算は商体の計算(ユークリッドの互助法またはグレブナー基底などを使う)よりも格段に易しく高速なので、区間の計算には望ましいでしょう。

```haskell
-- | 0 < .. < d^n < .. d^3 < d^2 < d < .. < 1/n < ...  <  1/2  <  1
prop_laurenOrderInfinitesimal n =
   n > 0   ==>   sorted [0,  d^(n+1), d^n, 1/n]
-- | 0 < 1 < 2 .. < n < .. < 1/d < 1/d^2 < 1/d^3 < .. < 1/d^n ...
prop_laurenOrderInfinitude    n =
   n > 0   ==>   sorted [n, 1/d, 1/d^n, 1/d^(n+1)]

sorted :: Ord a => [a] -> Bool
sorted xs = xs == sort xs
```

dを一次の無限小、d\^nをn次の無限小と呼びましょう。 二次の無限小`d^2`は、どんな整数nに対しても,`n * d^2 < d`が成りたちます。n次の無限小は、(n-1)次の無限小より「無限に小さい」というわけです。同様に、n次の無限大は(n-1)次の無限大より、「無限に大きい」と言えます。

Haskellの`Rational`や`Integer`で考えると、`Laurent Rational`や`Laurent Integer`はそれぞれ順序体、順序整域になり、これはいわゆる超実数(hyperreal)の一部分になっています。

これで晴れて、無限小と無限大の正しい扱いができました。

例えば、任意の区間の補集合の、「無限に大きい」上下二つの部分区間を定義してみます。

```haskell
type HyperRational = Laurent.T Rational
type HyperInterval = (HyperRational, HyperRational)

complement :: Integer
           -> HyperInterval
           -> (HyperInterval, HyperInterval)
complement n (a,b) = let
  ω = 1/d
  in ( (a - ω^n , a - d^n)
     , (b + d^n  , b + ω^n) )
```

これで、どんな`x::HyperRational`も適当な`complement n`でカバーできますし、有限部分にしか興味がないなら`complement 1`で間に合います。

-   Ex.1
    (`Ord a => Maybe a`)の公理を書きだして、無限小と無限大を加えて拡張してください。

-   Ex.2
    順序リスト(`Ord a => [a]`)の公理を書きだして、無限小と無限大を加えて拡張してください。

-   Ex.3
    Haskellの任意の順序モデルに対して、無限小と無限大を加えた順序モデルは存在するでしょうか。(単に`Ord`というだけでなく、他のクラスと順序の関係を保存すること)

パフォーマンスの考察
--------------------

numeric-preludeによるローラン環の定義では、ボックスの使用は避けられず、`Int64`や`Float32`に比べてパフォーマンスの低下が想定されます。

GHC 8.2.1から`-XUnboxedSums`という拡張がサポートされました。この拡張を使って、パフォーマンスを向上させられるか、考えてみます

### GHC 8.2.2のインストール

11月に、GHC8.2.2がリリースされ、各ディストリビューション向けにバイナリが用意されています。

### Debian 9

GHC-8.2.2のDebian用ビルドは[ここ](http://downloads.haskell.org/debian/)にあります。

### Ubuntu Zesty

```bash
$ # Should print the messages showing where the binaries are to be installed
$ sudo add-apt-repository ppa:hvr/ghc
$ sudo apt-get update
$ sudo apt-get install ghc-8.2.2
$ PATH=/opt/ghc/bin:$PATH
```

バイナリは`/opt/ghc/bin`以下にインストールされます。
インストール時のメッセージに注意してください。

## `UnboxedSums`拡張

ghc-8.2.1でサポートされた`UnboxedSums`拡張を使うと、直和型の不要なポインタ使用を省けます。文法は`UnboxedTuple`の類似で、

```
 data T a1 .. ak =  (# t_1 | t_2 | ... | t_N #)
```

のように書けます。`t_i`は任意の型で、ボックス型でもアンボックス型でもどちらでも大丈夫です。C/C++/機械語プログラマであれば、"タグ付構造体"と言えばすぐにパフォーマンスのイメージが掴めるでしょう。タグは32-bitです。ユニオンではなく構造体を使うのは、ガーベッジコレクタとの連携のためと思われます。不要な重複フィールドを省略するアルゴリズムを実装してある、とドキュメントにはありますが、未確認です。

最小例を示します。

```haskell
{-# LANGUAGE UnboxedSums #-}

data T a = (# Int | a #)

x0 :: T String
x0 = (# 0 | #) -- 最初の選択肢

x1 :: T String
x1 = (# | "abc" #) -- 二番目の選択肢

f :: T a -> Either Int String
f x = case x of
  (# i |   #) -> Left  i
  (#   | s #) -> Right s

main = do
    print $ f x0
    print $ f x1
```

ドキュメントは[こちら](https://downloads.haskell.org/~ghc/latest/docs/html/users_guide/glasgow_exts.html)。

この拡張を使い、ローラン環`HyperInteger` を最適化してみましょう。

```haskell
type HyperInteger = (#
{- |  Default value for 'small' enough finite values. 
      For simplicity we will restrict ourselves to 32-bit
      arithmetics (except the most significant bit), so  
      there will be no loss of information in arithmetic
      operations. 
-}
   Int64# 
{- | Fallback for infinite or infinitesimal values -}

   | Laurent Integer 
   #) 
```

この定義によると、HyperIntegerのメモリ上の表現は

```c
// This is a C code
struct HyperInteger {
  uint32_t tag;
  int64_t  first;
  box*     second;
};
```

のようになり、amd64/x86_64アーキテクチャでは(アラインメントによって)64*3ビットの構造体で表されます。メモリのサイズとしては最適ではないものの、IntelのCPUでは分岐予測が効きやすく、パフォーマンスが期待できそうです。

欠点としては、unboxed型の使用によってパラメトリックポリモーフィズムが失われることがあげられます。AssociatedTypes/TypeFamiliesなどの型関数を使い、アドホックに必要なケースを実装する必要があります。

もう一つの注意として、UnboxedSumは`union`ではなく`struct`として表されるので、使用頻度の低い選択肢を含めてしまうと、キャッシュを浪費してしまうことです。選択肢は可能な限り少数にすることが望ましいでしょう。

整数、実数を部分として含む構造は無数にあります。そのような構造の計算では、大半の時間を`Int#`や`Float#`として計算、オーバーフロー/アンダーフローやエラー時のみ他のボックスデータとして扱う、というのは自然で適用範囲の広い最適化といえます。将来的にはもっと使いやすい文法でGHCに実装されることを期待したいと思います。同様に`Integer`(巡回群の多項式)や`Rational`の高速化も期待できそうです。

結論
----

コンピュータの計算では、メモリは常に有限で、実数計算の精度は限られます。浮動点小数の計算では、NaNやエプシロン、無限大などの特別の値と性質がアドホックに定められ、足し算の結合法則すら保証されません。

ここでは、Haskellでは割合簡単に、またパフォーマンスへの犠牲を少なく、「正しい」無限小・無限大の扱いができることを示唆してみました。

ローラン環は大変高速で、実数の近似にも役立ちます。極限をとって形式的ローラン級数を考えれば体にもなります。`Rational`ほか商体の代替案として、プログラミング上もっと注目されてもいいのではないでしょうか。

[fn1]: 残念ながら`Double`や`Int`は、体や整域ではないので、任意の精度の計算が求められる場合は、 そのままは使えません。

[fn2]: 形式級数の体を考えることもできますが、コンピュータでは常に計算可能性の問題を考えないといけないので、別の用語が必要になります。残念ながら私は知りません(計算級数？)

[fn3]: モデル理論のコンパクト性定理を使うと、この操作を無限に繰替えして公理を無限に増やしても、なおその無限の公理を満たすモデルが存在することが言えます。これはとても強い定理で、無限の操作を許すと、有限の場合なら成り立つ性質(対称性など)が成り立たなくなることが多いのです。

## 補足：Haskellと型レベルの文法解析

Haskellでは、型クラスを使って型レベルのプッシュダウンオートマトンを実装するテクニックがあって、`openPar openBracket a b closePar closeBracket`などの文脈自由なシンボルの並びを処理することはできますが、残念ながら記号は変数として扱えません(エラーメッセージも理解不能なくらい複雑になりがち)。