LoginSignup
2
0

More than 3 years have passed since last update.

Prolog で複素数の計算

Last updated at Posted at 2019-06-12

目的

最近、github に公開されているとあるソースコードを眺めていると、数式の表現にテンソル積の記号 ⊗ や根号記号 √ が使われていて、驚くと共にその見易さに感銘を受けました。幸か不幸か日本人はマルチバイト文字にまみれた日常を送っていますが、プログラミングにマルチバイト文字を多用するという人、あるいはそうしたことが許容されている組織は少ないのではないでしょうか。

構文ルールを拡張可能なプログラミング言語、かつマルチバイト文字を識別子に使用可能な言語処理系を使えば、数式を手書きするときに近い形でプログラミングが可能になりそうです。

そういう可能性のあるプログラミング言語として、今回 Prolog に目を向けてみました。練習として、複素数計算を記述できるように目指してみます。

| ?- X  (4 + 2 * i) / (2 + i), printc(X), nl.
2.0

前提条件あるいは注意事項

ここでは、Unicode を使った項の表現が可能な Prolog 処理系の利用を前提としています。当方は SWI-Prolog および UTF-8 拡張を行った GNU prolog で動作確認しています。

複素数の内部表現方法

複素数は項 complex_value/2 で表現することにします。まずは、与えられた値が複素数かどうか判定する述語を書いてみます。

is_complex_term(complex_value(_, _)).

complex_value/2 の値について加減乗除を行う述語を準備します。いずれも、第一引数と第二引数を用いて演算を行い、その結果を第三引数に与える形式です。

% 加算
complex_add(
        complex_value(Ra, Ia),
        complex_value(Rb, Ib),
        complex_value(Rc, Ic)) :-
    Rab is Ra + Rb,
    Iab is Ia + Ib,
    Rc = Rab, Ic = Iab.
% 減算
complex_sub(
        complex_value(Ra, Ia),
        complex_value(Rb, Ib),
        complex_value(Rc, Ic)) :-
    Rab is Ra - Rb,
    Iab is Ia - Ib,
    Rc = Rab, Ic = Iab.
% 乗算
complex_mult(
        complex_value(Ra, Ia),
        complex_value(Rb, Ib),
        complex_value(Rc, Ic)) :-
    Rab is Ra * Rb - Ia * Ib,
    Iab is Ra * Ib + Ia * Rb,
    Rc = Rab, Ic = Iab.
% 除算
complex_div(A, B, C) :-
    complex_conj(B, Bc),
    complex_mult(A, Bc, complex_value(Rab, Iab)),
    complex_mult(B, Bc, complex_value(B2, _)),
    Rc is Rab / B2,
    Ic is Iab / B2,
    C = complex_value(Rc, Ic).

Prolog で特徴的なのは、数値計算を行うときに使う is/2 の左辺が未定義の変数でなければならないことです。計算とユニファイを独立して書いている部分はそういった事情によります。さて、上に書いた加減乗除のうち、除算は共役複素数を用いた計算を行っています。共役複素数を求める述語 complex_conj/2 を以下のように定義します。

% 共役複素数
complex_conj(
        complex_value(R, I),
        complex_value(R, Ic)) :-
    Im is I * -1, Ic = Im.

ここまであれば、複素数の計算が可能です。おつかれさまでした。

?- complex_div(complex_value(1, 0), complex_value(2, 3), Result).
Result = complex_value(0.15384615384615385, -0.23076923076923078).

これをベースとして、なるべく人に優しい自然な記法で記述できるように頑張ってみます。

見栄えを良くする

毎回 complex_なんとか、だとか complex_value だとか書く苦痛なので、以下のように書けるようにします。

?- compute(1 ÷ c(2, 3), Result), printc(Result).
0.15384615384615385-0.23076923076923078i
Result = complex_value(0.15384615384615385, -0.23076923076923078).

× や ÷ など、以下の演算子を使えるようにします。

% 共役複素数
:- op(150, xf, '^*').
% 乗算
:- op(400, yfx, '×').
% 除算
:- op(400, yfx, '÷').
% 加算
:- op(500, yfx, '+').

複素数を与えるときに complex_value(~) と書く代わりに c(~) で済ませることにします。1 なら c(1, 0)、i なら c(0, 1) という具合です。
c(~) で与えた複素数に対して加減乗除 (+, -, ×, ÷) したものを引数として与えると、それを複素数の計算式として評価してくれる compute/2 を定義してみたのが以下です。

% 複素数はそのまま
compute(Exp, Result) :- is_complex_term(Exp), !,
    Result = Exp.
% c(〜) は複素数として解釈する
compute(Exp, Result) :- Exp = c(R, I), !,
    Result = complex_value(R, I).
% 実数は複素数にする
compute(Exp, Result) :- number(Exp), !,
    Result = complex_value(Exp, 0).
% 負数
compute(Exp, Result) :- Exp = - E, !,
    compute(E, E0), compute(c(-1, 0) * E0, Result).
% 共役複素数
compute(Exp, Result) :- Exp = C ^*, !,
    compute(C, C0), complex_conj(C0, Result).
% 乗算
compute(Exp, Result) :- ( Exp = A × B ; Exp = A * B ), !,
    compute(A, A0),
    compute(B, B0),
    complex_mult(A0, B0, Result).
% 除算
compute(Exp, Result) :- ( Exp = A ÷ B ; Exp = A / B ), !,
    compute(A, A0),
    compute(B, B0),
    complex_div(A0, B0, Result).
% 加算
compute(Exp, Result) :- ( Exp = A  B ; Exp = A + B ), !,
    compute(A, A0),
    compute(B, B0),
    complex_add(A0, B0, Result).
% 減算
compute(Exp, Result) :- Exp = A - B, !,
    compute(A, A0),
    compute(B, B0),
    complex_sub(A0, B0, Result).

と書いたところで、やっぱり虚数は i として書きたいよね、と思ったので以下も追加します。

% i は虚数
compute(Exp, Result) :- Exp = i, !,
    Result = complex_value(0, 1).

さらに、自然に書けるように '←'/2 を複素数計算の演算子として定義することもできます。

:- op(900, xfx, '←').
'←'(Result, Exp) :-
    compute(Exp, Result).

計算結果を見やすくするために、複素数を表示する printc/1 を定義してみます。

printc(complex_value(R, I)) :-
    (
     % 虚部がなければ実部をそのまま表示
     (I =:= 0 ; I =:= 0.0) -> write(R), !
     ;
     (R =:= 0 ; R =:= 0.0) -> (
         % 虚部のみ、かつ虚部が 1 なら i のみ表示
         (I =:= 1 ; I =:= 1.0) -> write(i), !
         ;
         % 虚部のみ、かつ虚部が -1 なら -i のみ表示
         (I =:= -1 ; I =:= -1.0) -> write('-i'), !
         ;
         write(I), write(i), !
     )
     ;
     % 0 でない実部を表示
     write(R),
     (
      % 虚部が負なら、マイナス記号があるので実部に続けてそのまま表示できる
      (I < 0 ; I < 0.0) -> write(I), write(i), !
      ;
      % 虚部が正なら、実部との間にプラス記号を表示する
      write('+'), write(I), write(i) !
     )
    ).
printc(N) :-
    number(N), write(N), !.

こうして書いたものを complex.pl としてファイルにし、Prolog (Unicode が使えるものを使ってください)のプロンプト上で "[complex]." と投入すると complex.pl を読み込んでくれて、定義した compute/2 や printc/1 が利用可能になるはずです。

?- [complex].
true.

| ?- X  (2 + 3 * i) * (4 + 5 * i), printc(X), nl.   
-7+22i

X = complex_value(-7,22)

yes
| ?- X  (-7 + 22 * i) / (4 + 5 * i), printc(X), nl. 
2.0+3.0i

X = complex_value(2.0,3.0)

yes
2
0
0

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
2
0