数式処理システムとしてのプログラミング言語 Egison

Egison version 3.6.0から、数式処理機能が実装されました。 Egison独自のパターンマッチを数式に対して応用することにより、ユーザーによって簡単に拡張可能なシステムであることが、この数式処理システムの大きな特徴となっています。このEgison上に実装された数式処理システムについてご紹介します。

シンボリックな計算とは

シンボリックな計算とは、x + x = 2xや(x + y)^2 = x^2 + 2 x y + y^2のようにxやyといったシンボルをプログラム上で数と同様に扱い計算することをいいます。シンボリックな計算をサポートしている既存の数式処理システムとしては、Mathematicaや、Maxima、Pythonの拡張ライブラリとして実装されているSymPyなどが有名です。

最新のEgisonは以下のようにシンボリックな計算をサポートしています。

; (x + y)の2乗
(* (+ x y) (+ x y))
;=>(+ x^2 (* 2 x y) y^2)

; (1 + i)の5乗
(** (+ 1 i) 5)
;=>(+ -4 (* -4 i))

i^2を-1に変換するような書き換え規則は、Egisonのパターンマッチを用いてライブラリで定義されています。ユーザーが新たに書き換え規則を追加することも可能です。

シンボリックな計算を用いることにより、代数的数を浮動小数点数で近似することなく、扱うことができます。例えば、以下のように、二次方程式の解の公式を用いて、二次方程式の解を計算することができます。

(define $q-f
  (lambda [$a $b $c]
    [(/ (+ (* -1 b) (sqrt (- (** b 2) (* 4 a c)))) (* 2 a))
     (/ (- (* -1 b) (sqrt (- (** b 2) (* 4 a c)))) (* 2 a))]))

(q-f 1 1 -1)
;=>[(/ (+ -1 (sqrt 5)) 2) (/ (+ -1 (* -1 (sqrt 5))) 2)]

この発展的な例として、幾つかのnについて1のn乗根を計算したプログラムを用意しています。

パターンマッチによる微分の定義

Egisonによる数式処理システムの一番大きな特徴は、Egison独自のパターンマッチを数式に対して適用できることにあります。例えば、微分はEgisonのパターンマッチにより以下のように数十行のプログラムで定義できます。

(define $∂/∂
  (lambda [$f $x]
    (match f math-expr
      {; symbol
       [,x 1]
       [?symbol? 0]
       ; function application
       [(,exp $g) (* (exp g) (∂/∂ g x))]
       [(,log $g) (* (/ 1 g) (∂/∂ g x))]
       [(,cos $g) (* (* -1 (sin g)) (∂/∂ g x))]
       [(,sin $g) (* (cos g) (∂/∂ g x))]
       [(,sqrt $g) (* (/ 1 (* 2 (sqrt g))) (∂/∂ g x))]
       [(,** $g $h) (* f (∂/∂ (* (log g) h) x))]
       [<apply $g $args>
        (sum (map 2#(* (capply `g|%1 args) (∂/∂ %2 x))
                  (zip nats args)))]
       ; quote
       [<quote $g> (∂/∂ g x)]
       ; term (constant)
       [,0 0]
       [(* _ ,1) 0]
       ; term (multiplication)
       [(* ,1 $fx^$n) (* n (** fx (- n 1)) (∂/∂ fx x))]
       [(* $a $fx^$n $r)
        (+ (* a (∂/∂ (**' fx n) x) r)
           (* a (**' fx n) (∂/∂ r x)))]
       ; polynomial
       [<poly $ts> (sum (map (∂/∂ $ x) ts))]
       ; quotient
       [(/ $p1 $p2)
        (let {[$p1' (∂/∂ p1 x)]
              [$p2' (∂/∂ p2 x)]}
          (/ (- (* p1' p2) (* p2' p1)) (** p2 2)))]
       })))

(define $d/d ∂/∂)

この微分関数は以下のように動作します。

(d/d (** x 3) x)
;=>(* 3 x^2)

(d/d (** e (* i x)) x)
;=>(* i (** e (* i x)))

(d/d (* (cos x) (sin x)) x)
;=>(+ (* -1 (sin x)^2) (cos x)^2)

(∂/∂ (f x y) x)
;=>(f|1 x y)

(∂/∂ (f x y) y)
;=>(f|2 x y)

(d/d (r (u t) (v t)) t)
;=> (+ (* (r|1 (u t) (v t)) (u|1 t))
;      (* (r|2 (u t) (v t)) (v|1 t)))

この微分関数を利用して、入力された関数のテイラー展開を計算するプログラムは以下のように定義できます。

(define $taylor-expansion
  (lambda [$f $x $a]
    (map2 *
          (map 1#(/ (** (- x a) %1) (fact %1)) nats0)
          (map (substitute {[x a]} $) (iterate (d/d $ x) f)))))

(take 4 (taylor-expansion (f x) x 0))
;=>{(f 0) (* x (f|1 0)) (/ (* x^2 (f|1|1 0)) 2) (/ (* x^3 (f|1|1|1 0)) 6)}

このプログラムを利用してオイラーの公式を確かめることができます。

(take 5 (taylor-expansion (** e (* i x)) x 0))
;=>{1 (* x i) (/ (* -1 x^2) 2) (/ (* -1 x^3 i) 6) (/ x^4 24)}

(take 5 (taylor-expansion (cos x) x 0))
;=>{1 0 (/ (* -1 x^2) 2) 0 (/ x^4 24)}

(take 5 (taylor-expansion (* i (sin x)) x 0))
;=>{0 (* x i) 0 (/ (* -1 x^3 i) 6) 0}

またもう少し複雑な微分の計算例として、極座標のラプラシアンの検算と導出をしています。

このようにパターンマッチの表現力を活かして、ここでデモした微分以外の色々な数学の理論についても、数十行のプログラムで実装し、実験できる処理系を目指して開発を進めています。

クオート式による式展開の制御

数式処理システムの設計において、式展開をどのように制御するかは非常に重要な問題です。式展開の制御について、Egisonによる数式処理システムには、以下に述べるような動作をするクオート式の導入というオリジナルの工夫が実装されています。

'が先頭に付いている式をクオート式と呼びます。 'に続く式はクオートされた式と呼ばれます。クオートされた式は評価され展開されます。しかし、その展開された式は1つのまとまった式として扱われ、他の式と足したり、掛けたりしても、それ以上は展開されません。

(** '(+ x y) 2)
;=>'(+ x y)^2

'(** (+ x y) 2)
;=>'(+ x^2 (* 2 x y) y^2)

(+ '(+ x y) x y)
;=>(+ '(+ x y) x y)

(+ '(+ x y) '(+ x y))
;=>(* 2 '(+ x y))

(+ '(+ x y) '(+ y x))
;=>(* 2 '(+ x y))

意外にも、このクオート式を用いるだけでも、まるでノートで計算するように式の展開を制御できることがあります。以下の例はクオート式が効果的に使われているある程度大きな規模の計算の例です。計算の要所々々にクオート(')を挿入することにより、長い計算でも自動で紙の上で計算するように制御されています。

組み込みデータ型としてテンソルのサポート

テンソルの扱いについていくつか細かい工夫をしています。 Egisonは、タプルや配列と同様に、組み込みデータ型としてテンソルをサポートしています。そうすることにより、ベクトル解析や、テンソル解析で使われているテンソル計算の表記を、ほぼそのままプログラミングに持ち込んでいます。

いくつか具体例を見てみましょう。

2引数関数の引数の両方が、ベクトルである場合をまず見ていきます。ベクトルは[|、|]で要素を囲むことにより表現されます。またベクトルの後に、_とそれに続いて数をつなげることにより、添字付きのベクトルを表現できます。添字は自然数であることもあれば、シンボルであることもあります。それぞれのベクトルに対する添字が同じシンボルである場合、以下のように、対応する要素ごとに計算が行われます。

(+ [| 1 2 3 |]_i [| 1 2 3 |]_i)
;=>[|2 4 6|]_i

添字が指定されていない場合、テンソル積の形で計算されます。それぞれの添字が異なる場合も、テンソル積の形で計算されます。この場合、ベクトルとベクトルからは行列が生成されます。

(+ [| 10 20 30 |] [| 1 2 3 |])
;=>[| [| 11 12 13 |] [| 21 22 23 |] [| 31 32 33 |] |]

(+ [| 10 20 30 |]_i [| 1 2 3 |]_j)
;=>[| [| 11 12 13 |] [| 21 22 23 |] [| 31 32 33 |] |]_i_j

スカラーとベクトル・ベクトルと行列のように階数が異なるテンソルの演算も以下のように可能になっています。

(+ 1 [| 1 2 3 |])
;=>[|2 3 4|]

(+ [| 100 200 300 |]_i
   [|[| 11 12 |]
     [| 21 22 |]
     [| 31 32 |]|]_i_j)
;=>[| [| 111 112 |] [| 221 222 |] [| 331 332 |] |]_i_j

上記の例は、関数の引数がテンソルである場合でしたが、関数がテンソルである場合にも対応する自然な処理があります。

([| (∂/∂ $ x) (∂/∂ $ y) |] (f x y))
;=>[| (f|1 x y) (f|2 x y) |]

内積を計算するには、contract式を用います。アインシュタインの縮約記法を実現するようにこのcontract式の挙動は設計されており、同じシンボルの上添字と下添字がある場合、その箇所をcontract式の第一引数で指定された演算により縮約します。

(contract + (* [| 1 2 3 |]~i [| 10 20 30 |]_i))
;=>
140

これらの記法を組み合わせると、ベクトル解析において重要な演算子、ナブラやdiv(発散)、rot(回転)なども簡潔に記述できます。

(define $∇ ∂/∂)

(define $div
  (lambda [%A %xs]
    (trace (∇ A xs))))

(define $rot
  (lambda [%A %xs]
    (cross-product/fn ∂/∂ A xs)))

ある程度の規模になる計算の例として、リーマン曲率テンソル(空間の曲がり具合を表現する4階のテンソル)を計算するプログラムがあります。リーマン曲率テンソルは一般相対性理論の記述の際に重要な役目を果たすテンソルです。具体的にこのテンソルの値を計算するのは少々煩雑なのですが、Egisonを使うと、教科書に載っている定義を記述するだけで、計算することができます。

今後の展開

方程式のサポート

上に挙げた例はすべて方程式の解を導出するプログラムでしたが、処理対象のデータとして方程式もサポートし、方程式自体も導出できるようになる予定です。また、方程式に対するパターンマッチを適用することにより、微分方程式を分類したり、変形し解くプログラムも容易に記述することができるようになります。

グラフィックのサポート

既存の数式処理システムと同様にグラフの描画もサポートする予定です。

次に見てほしいもの...

Egison数学ノートトップに戻る