9数式処理システム入門

数式処理システムとは,シンボリックな計算をサポートするプログラミング言語のことをいう. 数式処理システムは,未束縛の変数をシンボルとして扱う. 数式処理システムを使うと,たとえば,\(x + x \rightarrow 2x\)\((x + y)^2 \rightarrow x^2 + 2 x y + y^2\)のような計算ができる. Egisonもこのようなシンボリックな計算をサポートしている. 本章は数式処理システムを使うとどのような計算ができるのかみていく.

9.1シンボルの宣言

数式処理システムは,未束縛の変数をシンボルとして扱う. Egisonでは,シンボルを使う前にdeclare symbol文で宣言する必要がある. declare symbol x : MathValueのように型を指定することもできる(数式と型の関係については第11章で解説する).

declare symbol x, y, z

x -- x

数式処理システムとしての機能をもつEgisonには,シンボル同士の足し算や掛け算が組み込みで定義されている. そのため,たとえば,\(x + x \rightarrow 2x\)\((x + y)^2 \rightarrow x^2 + 2 x y + y^2\)のような計算ができる.

x + x -- 2 * x
(x + y)^2 -- y^2 + x^2 + 2 * x y

Egisonは数式を自動で積和標準形に展開する. 積和標準形とは,掛け算が内側に,足し算が外側になるような形の数式のことである. そのため,\((x + y)^2\)\(x^2 + 2 x y + y^2\)のかたちに展開される. なお,出力では,係数と因子の間の掛け算は*で表示され,シンボル同士の掛け算は空白区切りで表示される.

出力の数式の項や因子の並び順は,ユーザーが入力した数式に現れる順番ではなく,処理系の定める正準的な順序に整列される. たとえば,上記の\((x + y)^2\)\(x\)\(y\)を入れ替えて\((y + x)^2\)を計算しても,出力は同じになる.

(y + x)^2 -- y^2 + x^2 + 2 * x y

9.2関数適用のシンボル化

sqrt xexp xのようにシンボルを引数にとる可能性がある関数について,関数適用を止めて,その式をシンボリックな値としてそのまま返すようにライブラリで定義されている. シンボリックな値として残った関数適用は,先頭にシングルクオート(')を付けて表示される.

declare symbol x

sqrt 4 -- 2
sqrt x -- 'sqrt x
exp 1 -- e
exp x -- 'exp x

sqrt 2のように評価を進められない関数適用についても,関数適用を止めて,その式をそのまま返すようにライブラリで定義されている. \(\sqrt{8} \rightarrow 2 \sqrt{2}\)のような簡約もライブラリで定義されている.

sqrt 2 -- 'sqrt 2
sqrt 8 -- 2 * 'sqrt 2

入れ子になった根号も,外せる場合(\(\sqrt{a+b\sqrt{c}}\)\(a^2-b^2c\)が平方数の場合)には自動的に外される.

sqrt (9 - 4 * sqrt 5) -- 'sqrt 5 - 2
sqrt (2 + sqrt 2) -- 'sqrt ('sqrt 2 + 2)

このような数学関数はdeclare mathfunc文で宣言されており,適用時の簡約はdeclare apply文で定義されている. これらの宣言文は第1111.3節で,定義に使われるシングルクオートは10.3節で解説する.

9.3数式の簡約

\(\sqrt{x} \cdot \sqrt{x} = x\)\(sin^2\theta + cos^2\theta = 1\)のような数式を紙の上で扱うときにおこなう簡約が,Egisonを含む多くの数式処理システムには実装されている.

declare symbol x, `$\theta$`

sqrt x * sqrt x -- x
(sin `$\theta$`)^2 + (cos `$\theta$`)^2 -- 1

Egisonの簡約機能は,二つの層に分かれて実装されている. 数式の内部表現と,数式を積和標準形へ正規化する処理は,Sweet Egison(Egisonのパターンマッチ機能を提供するHaskellライブラリ,第8章)を使ってHaskellによりEgison処理系の内部に実装されている. 一方,上記のような個々の簡約規則の大部分は,Haskellコードとしてではなく,Egison自身のdeclare rule文による簡約規則の宣言として,標準ライブラリで定義されている.

複素数などについての簡約規則も同様に定義されている.

i^2 -- -1
w^3 -- 1
w + w^2 -- -1

たとえば,虚数単位\(i\)と1の原始3乗根\(w\)についての上記の簡約は,標準ライブラリのlib/math/normalize.egihttps://github.com/egison/egison/blob/master/lib/math/normalize.egi)にある以下の宣言により実現されている.

declare rule auto term i^2 = -1
declare rule auto term w^3 = 1
declare rule auto term w^2 = -1 - w

ユーザーもこれとまったく同じdeclare rule文を使って,簡約規則を自分で宣言できる(第1111.1節). なお,微分形式の計算のように性能が特に重要になる一部の簡約規則だけは,例外的にHaskell側(hs-src/Language/Egison/Math/Rewrite.hs)に実装されている.

多項式の分数の約分は,多項式GCD(単変数はユークリッドの互除法,多変数はsubresultant PRS)により自動でおこなわれる.

(x^2 - 1) / (x - 1) -- x + 1
(x^2 - y^2) / (x - y) -- y + x

イデアルを法とする正規形も,グレブナー基底の理論に基づいて計算できる. 関係式の集合から完備な簡約規則の集合を機械生成するdeclare ideal文と,正規形を明示的に計算するpolyNF関数が用意されている(第1111.2節). 一方,数式の簡約の完全な自動化は非常に難しい問題であり,Egisonが複雑な数式を適切に簡約できないこともある.

Wolfram言語(Mathematica)は,このような高度な数式の簡約ができる. “-M mathematica”オプションを追加してEgison処理系を実行するとWolfram言語で読み込めるかたちで数式を出力する. この出力をWolfram言語の処理系に入力すれば,Wolfram言語の高度な簡約機能を利用できる.

$ egison -M mathematica
> declare symbol x
> x + sqrt 2
#mathematica|x + Sqrt[2]|#

簡約規則や数学関数をユーザーが宣言する方法や,型注釈による正規形の指定,型の昇格タワーとその拡張,商の型など,簡約システムの詳細は第11章で解説する.

9.4二次方程式を解くプログラム - 数式を処理するアルゴリズムの記述

本節では,二次方程式を解くアルゴリズムをプログラムとして書いていく. 以下のように二次方程式を引数にとり,解を返す関数qFを定義していく.

declare symbol x, a, b, c

qF (x ^ 2 + x + 1) x
-- ((1 / 2) * ('sqrt 3) * i + -1 / 2,
--  (-1 / 2) * ('sqrt 3) * i + -1 / 2)

qF (x ^ 2 + b * x + c) x
-- ((-1 / 2) * b + (1 / 2) * 'sqrt (b^2 - 4 * c),
--  (-1 / 2) * b + (-1 / 2) * 'sqrt (b^2 - 4 * c))

qF (a * x ^ 2 + b * x + c) x
-- ((-1 / 2) * a^-1 b + (1 / 2) * ('sqrt (b^2 - 4 * a c)) * a^-1,
--  (-1 / 2) * a^-1 b + (-1 / 2) * ('sqrt (b^2 - 4 * a c)) * a^-1)

qF (a * x ^ 2 + 2 * b * x + c) x
-- (- a^-1 b + ('sqrt (b^2 - a c)) * a^-1,
--  - a^-1 b - ('sqrt (b^2 - a c)) * a^-1)

出力は,解の公式\(\frac{-b \pm \sqrt{b^2 - 4 a c}}{2 a}\)の分数を各項に分配した積和標準形で表示される. また,\(\frac{1}{a}\)\(a^{-1}\)という負の冪として表示される(Egisonの数式の内部表現は,負の冪指数を許すローラン多項式になっている). qFは解の公式をそのまま書くだけでも定義することができる. 以下のプログラムはqF'という補助関数を使ってqFを定義している. qF'は第一引数に二次の係数,第二引数に一次の係数,第三引数に定数項をとる. たとえば,qF' 1 0 2は2次方程式\(x^2 + 2 = 0\)の解を返す. coefficientsは多項式とシンボルを引数にとり,そのシンボルについて係数のリストを返す関数である. qF'は単純に\(\frac{-b + \sqrt{b^2 - 4 a c}}{2 a}\)\(\frac{-b - \sqrt{b^2 - 4 a c}}{2 a}\)からなるタプルを返す関数として定義されている.

def qF f x :=
  match coefficients f x as list mathValue with
    | [$a_0, $a_1, $a_2] -> qF' a_2 a_1 a_0

def qF' a b c :=
  ( ((- b) + sqrt (b ^ 2 - 4 * a * c)) / (2 * a)
  , ((- b) - sqrt (b ^ 2 - 4 * a * c)) / (2 * a) )

しかし,上記の記述ではアルゴリズムっぽさがないため,今回は平方完成によりこの解の公式を求めるアルゴリズムを書いていく. 平方完成とは,\(x^2 + b x + c\)のような二次式を\((x + b / 2)^2 - b^2 / 4 + c\)のように変形する操作のことである.

def qF' a b c :=
  match (a, b, c) as (mathValue, mathValue, mathValue) with
    | (#1, #0, _) -> (sqrt (- c), - sqrt (- c))
    | (#1, _, _) ->
      let (r1, r2) := withSymbols [x, y]
                        qF (substitute [(x, y - b / 2)] (x ^ 2 + b * x + c)) y
       in ((- (b / 2)) + r1, (- (b / 2)) + r2)
    | (_, _, _) -> qF' 1 (b / a) (c / a)

平方完成を使ってを二次方程式を解くアルゴリズムを記述するとqF'は3つのマッチ節からなる関数となる. 3行目のマッチ節は,\(x^2 + c = 0\)のような形の二次方程式を処理する. このマッチ節は,\(\sqrt{-c}\)\(-\sqrt{-c}\)という結果を返す. 4-7行目のマッチ節は,\(x^2 + b x + c = 0\)のような二次の係数が\(1\)である形の二次方程式を処理する. このマッチ節が平方完成の操作をするこの関数の本体である. substituteはEgisonのライブラリ関数である. substituteは第一引数で指定された代入を第二引数の式に対しておこなう. ここでは第一引数で,(x, y - b / 2)という代入が指定されているため,\(x = y - b / 2\)という代入を式\(x^2 + b x + c\)に対しておこなう. 同時に複数の代入を指定するために,substituteの第一引数はリストを引数に取るようになっている. \(x = y - b / 2\)という代入をするとこの方程式は3行目のマッチ節でマッチする\(y\)についての二次方程式に変形される. そのため,qF関数にこの変形された方程式を渡せば,\(y\)についてこの二次方程式の解を得られる. この\(y\)についての解に\(\frac{-b}{2}\)を加えることにより,もとの\(x\)についての解を得ることができる. withSymbols式は,局所シンボルを生成する構文である. 第一引数のリストで指定されたシンボル(上記の場合,xy)を第二引数の式の中でシンボルとして使うことができる. withSymbols式のこの機能のおかげで,もしプログラムの別の箇所で,xyにたとえば具体的な整数が束縛されていたとしても,withSymbols式の内部では,その影響を考えなくてよくなる. 8行目のマッチ節は,二次の係数が\(1\)でない二次方程式を処理する. 二次の係数で方程式を割ることにより,二次の係数が\(1\)の方程式に変換し,ふたたびqF関数にわたす処理をしている.

二次方程式より一段と複雑になるが,三次方程式・四次方程式を解くプログラムを書いてみると面白い. これらのプログラムはEgisonのソースコードのsample/math/algebraディレクトリ以下に公開されている.

9.5微分するプログラム - 数式に対するパターンマッチ

Egisonには微分計算をするための関数d/dがライブラリ関数として実装されている. この関数は第一引数の関数を第二引数のシンボルについて微分した結果を返す. 関数の名前をd/dのようにすることで数式に近い形で微分を表現できている. Egisonは/を変数名に使うことを許している.

d/d x x -- 1
d/d (x^2) x -- 2 * x
d/d (exp x) x -- 'exp x
d/d (log x) x -- x^-1
d/d (x * log x) x -- 'log x + 1
d/d (1 / log x) x -- - 'log x^-2 * x^-1

ライブラリのd/dは,テンソル(第12章)にもそのまま適用できるように型クラスを使って構成されており,新しい数学関数の微分公式は11.3節で紹介するdeclare derivative文で追加できるようになっている. しかし,その中核にあるのは数式に対するパターンマッチである. 本節では,この中核部分を切り出して,d/dと同じように動作する関数d/d'を定義してみる.

d/d'は以下のように定義できる.

def d/d' (f: MathValue) (x: MathValue) : MathValue :=
  match f as mathValue with
    -- Differentiation of symbols
    | #x -> 1
    | symbol _ _ -> 0
    -- Differentiation of function applications
    | apply1 #exp $g -> exp g * d/d' g x
    | apply1 #log $g -> 1 / g * d/d' g x
    | apply1 #sqrt $g -> 1 / (2 * sqrt g) * d/d' g x
    | apply2 #(^) $g $h -> f * d/d' (log g * h) x
    | apply1 #cos $g -> (- sin g) * d/d' g x
    | apply1 #sin $g -> cos g * d/d' g x
    -- Differentiation of constants
    | #0 -> 0
    | _ * #1 -> 0
    -- Differentiation of terms
    | #1 * $fx ^ $n -> n * fx ^ (n - 1) * d/d' fx x
    | $a * $fx ^ $n * $r -> a * d/d' (fx ^' n) x * r + a * fx ^' n * d/d' r x
    -- Differentiation of polynomials
    | poly $ts -> sum (map 1#(d/d' $1 x) ts)
    -- Differentiation of quotients
    | $p1 / $p2 ->
      let p1' := d/d' p1 x
          p2' := d/d' p2 x
       in (p1' * p2 - p2' * p1) / p2 ^ 2

第一引数の\(f\)を数式としてパターンマッチしている. 4-5行目はシンボルの微分を定義している. もし,fxそのものであったら,\(1\)を返し,fx以外のシンボルであったら\(0\)を返すように定義されている. 7-12行目はいくつかの基本的な関数の適用について合成関数の微分\(\frac{df(g(x))}{dx} = \frac{df(g(x))}{dg(x)} \frac{dg(x)}{x}\)を記述している. たとえば,7行目は指数関数の微分をしている. このマッチ節の本体は\(\frac{de^{g(x)}}{dx} = e^{g(x)} \frac{dg(x)}{dx}\)を表現している. また,8行目は対数関数の微分をしている. このマッチ節の本体は\(\frac{dlog(g(x))}{dx} = \frac{1}{g(x)} \frac{dg(x)}{dx}\)を表現している. 14-15行目は定数項に対する微分を定義している. 定数項は微分すると\(0\)になる. 17-18行目は,ここまでに微分を定義した項を掛け合わせた項についての微分を定義している. \(x^n\)\(n x^{n-1}\)に微分されることと,\(f g\)\(f' g + f g'\)に微分されることを利用して定義している. ここでも合成関数の微分の公式が使われている. 20行目は多項式の微分が定義されている. ここまでで単項式についての微分が定義されているので,それらを足し合わせたものについては,それぞれの項について微分して足し合わせてやればよい. 22-24行目は多項式を分母と分子にもつ式に対する微分が定義されている. 商の微分の公式を使って定義されている.

9.6ベクトルや行列の計算 - 添字記法

Egisonの数式処理システムとしての特徴に,ベクトルや行列,そしてそれらを一般化したテンソルを扱う計算を簡潔に記述できるというものがある. これは二十世紀初頭に微分幾何や相対性理論の研究の過程で発明されたテンソルの添字記法をプログラミング言語の機能としてEgisonがサポートしているためである.

添字記法は,テンソルに添字を付加することによって,さまざまな種類のテンソルの掛け算を表現する. たとえば,ベクトル同士の掛け算には,テンソル積・アダマール積・内積の三種類があるが,これらは以下のように添字を使い分けることによって表現される.

declare symbol x1, x2, y1, y2, i, j

[| x1, x2 |]_i .[| y1, y2 |]_j
-- [| [| x1 y1, x1 y2 |], [| x2 y1, x2 y2 |] |]_i_j
[| x1, x2 |]_i .[| y1, y2 |]_i -- [| x1 y1, x2 y2 |]_i
[| x1, x2 |]~i .[| y1, y2 |]_i -- x2 y2 + x1 y1

上記の例のように,Egisonではベクトルは成分を[| |]で囲むことにより表現する. また行列は,ベクトルのベクトルとして表現する. 関数型プログラミングの視点で添字記法は,テンソルの掛け算関数のパラメータとして,テンソルだけでなく添字も追加で渡すことによって一つの関数に複数の役割を持たしていると説明できる.

添字記法を言語機能としてもつプログラミング言語はほかにも存在する. そのなかでEgisonの特徴は,関数型プログラミングとうまく調和するかたちで添字記法を組み込んでいるところである. そのおかげで,高階関数と添字記法を組み合わせて使うことができる. たとえば,\(g_{i_{1}j_{1}} g_{i_{2}j_{2}} ... g_{i_{n}j_{n}}\)のような数式をfoldl関数を使って以下のように記述することができる.

foldl . 1 (map (\x -> g_[i_x]_[j_x]) [1..n])

添字記法はEgisonの重要な言語機能であるため,第12章と第14章で詳しく紹介する.

9.7数式の出力形式

-M”オプションで出力形式を指定することができる. “egison -M latex”を実行すると,LaTeX形式で結果を出力する.

$ egison -M latex
Egison Version 5.0.0
https://www.egison.org
Welcome to Egison Interpreter!
> declare symbol x
> x + sqrt 2
#latex|x + \sqrt{2}|#

egison -M mathematica”を実行すると,Mathematica形式で結果を出力する.

$ egison -M mathematica
Egison Version 5.0.0
https://www.egison.org
Welcome to Egison Interpreter!
> declare symbol x
> x + sqrt 2
#mathematica|x + Sqrt[2]|#

LaTeX出力はJupyter Notebookとも連携している. Jupyter NotebookのEgisonプラグインを利用すれば,図14.2,図14.4のようにインタラクティブに計算結果が数式として表示される.

Jupyter Notebook上でEgisonを実行する様子

この本を別の言語で読む: English, 日本語