目的
近年、量子コンピュータの開発が実用化に向けて進んでいます。量子コンピュータには大きく分けて量子アニーリング方式と量子ゲート方式とがあります。
ここでは、量子ゲート方式での量子計算を Prolog 上で記述・シミュレートする方法を模索してみます。こんな感じで書けるようにしたいと思います：
?- X ← σx │↓〉.
X = '│↑〉'.
?- X ← hadamard(σx │↑〉), q_write(X), nl.
0.7071067811865475*│↑〉-0.7071067811865475*│↓〉
X = complex_value(0.7071067811865475, 0.0)* ('│↑〉')-complex_value(0.7071067811865475, 0.0)* ('│↓〉').
ブラベクトル、ケットベクトルの表記
SWI-Prolog 上でいろいろ実験した結果、以下のような表現にすれば無難にアトムとして扱えそうだという結論に至りました。
'│↑〉'.
'│↓〉'.
'〈↑│'.
'〈↓│'.
なんとか
∣0⟩
∣1⟩
⟨0∣
⟨1∣
あるいはこれに近い見た目のものにしたかったのでとても残念です。
しぶしぶながらも上記の矢印を使ったごつい表記で頑張ることにします。
便利な場合もありそうなので、述語で分類して定義しておきます：
% 基底ベクトルはブラとケットに分類できる
q_basis(B) :- bra_basis(B).
q_basis(K) :- ket_basis(K).
% ブラ基底
bra_basis(〈↑│).
bra_basis(〈↓│).
% ケット基底
ket_basis(│↑〉).
ket_basis(│↓〉).
ブラベクトルの並び、ケットベクトルの並び、あるいはケットとブラからなる演算子を表現できるように、すべて fy 形式の演算子としても利用可能にしておきます。
:- op(200, fy, │↑〉).
:- op(200, fy, │↓〉).
:- op(200, fy, 〈↑│).
:- op(200, fy, 〈↓│).
これで、なんとなくそれらしい表記が可能になりました。
?- X = │↓〉 〈↑│ + │↑〉 〈↓│.
X = '│↓〉' '〈↑│'+'│↑〉' ('〈↓│').
式を書くときに空白をちゃんと入れてやらないと、例えば │↓〉〈↑│ と空白を入れずに書くと '│↓〉〈↑│ ' というアトムとみなされてしまうのが難点です。
いろいろ難点だらけで気が重くなってきますが、さらに話を進めます。
計算の定義
前に書いた複素数計算と同じ感じで進めていきます。まずは式全体を表現する述語 '←'/2 を用意します：
:- op(900, xfx, '←').
'←'(Result, Exp) :-
q_compute(Exp, Result), !.
q_compute の中身を適宜拡張していけば目的の計算ができるシステムが出来上がるはずです。
準備として、複素数を扱うための述語を定義しておきます：
% 複素数あるいは数値であるか検査する述語
is_complex_number(complex_value(_,_)). % 複素数
is_complex_number(i). % 虚数
is_complex_number(N) :- number(N). % 組み込み型の数値
計算結果が complex_value なんとかとかだと今ひとつ見づらいので、見やすく表示する述語を用意します：
q_write(N) :-
number(N), write(N), !.
q_write(M1*B) :-
is_minus_one(M1),
write('-'), q_write(B), !.
q_write(A*B) :-
q_write_bracket(A), write('*'), q_write_bracket(B), !.
q_write(A+B) :-
q_write(A), write('+'), q_write(B), !.
q_write(A-B) :-
q_write(A), write('-'), q_write_bracket(B), !.
q_write(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), !
)
).
q_write(X) :- write(X), !.
q_write_bracket(A+B) :-
write('('), q_write(A), write('+'), q_write(B), write(')'), !.
q_write_bracket(A-B) :-
write('('), q_write(A), write('-'), q_write(B), write(')'), !.
q_write_bracket(A*B) :-
q_write(A), write('*'), q_write(B), !.
q_write_bracket(X) :-
q_write(X).
複素数に関して、後で必要になるユーティリティ述語も定義してみます。
% 与えられた組み込み型数値を複素数型に変換する述語
to_complex_value(N, complex_value(N, 0)) :- number(N).
to_complex_value(i, complex_value(0, 1)). % i は虚数
to_complex_value(C, C) :- is_complex_number(C). % 複素数はそのまま
% 引数が 0 かどうか判定する述語：
is_zero(0). is_zero(0.0).
is_zero(complex_value(0, 0)). is_zero(complex_value(0, 0.0)).
is_zero(complex_value(0.0, 0)). is_zero(complex_value(0.0, 0.0)).
% 1 かどうか判定する述語：
is_one(1). is_one(1.0).
is_one(complex_value(1, 0)). is_one(complex_value(1.0, 0)).
% -1 かどうか判定する述語：
is_minus_one(-1). is_minus_one(-1.0).
is_minus_one(complex_value(-1, 0)). is_minus_one(complex_value(-1.0, 0)).
% 数値の掛け算：
complex_mult(A, B, complex_value(Rr, Ir)) :-
to_complex_value(A, complex_value(Ra, Ia)),
to_complex_value(B, complex_value(Rb, Ib)),
Rr0 is Ra * Rb - Ia * Ib, Rr = Rr0,
Ir0 is Ra * Ib + Ia * Rb, Ir = Ir0.
% 数値の割り算：
complex_div(A, B, Result) :-
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,
Result = complex_value(Rc, Ic).
complex_conj(
complex_value(R, I),
complex_value(R, Ic)) :-
Im is I * -1, Ic = Im.
% 数値の足し算：
complex_add(A, B, complex_value(Rr, Ir)) :-
to_complex_value(A, complex_value(Ra, Ia)),
to_complex_value(B, complex_value(Rb, Ib)),
Rr0 is Ra + Rb, Rr = Rr0,
Ir0 is Ia + Ib, Ir = Ir0.
% 数値の引き算：
complex_sub(A, B, complex_value(Rr, Ir)) :-
to_complex_value(A, complex_value(Ra, Ia)),
to_complex_value(B, complex_value(Rb, Ib)),
Rr0 is Ra - Rb, Rr = Rr0,
Ir0 is Ia - Ib, Ir = Ir0.
% 平方根の計算：
:- op(180, fx, '√').
complex_sqrt(A, complex_value(Rr, Ir)) :-
to_complex_value(A, complex_value(Ra, Ia)),
( Ia >= 0 ->
Rr0 is sqrt((sqrt(Ra*Ra+Ia*Ia) + Ra) / 2),
Ir0 is sqrt((sqrt(Ra*Ra+Ia*Ia) - Ra) / 2), !
; Rr0 is sqrt((sqrt(Ra*Ra+Ia*Ia) + Ra) / 2),
Ir0 is -sqrt((sqrt(Ra*Ra+Ia*Ia) - Ra) / 2)
), Rr = Rr0, Ir = Ir0.
ここまでで複素数に関する操作の準備ができました。
ここからは記号操作を含めた計算操作を定義していきます。
% 数値（複素数含む）型の値はそのまま返す
q_compute(C, Result) :- is_complex_number(C),
to_complex_value(C, Result).
次に、掛け算を定義してみます。
- 0 との掛け算結果は 0
- 1 との掛け算結果は 元のまま
- 数値計算できるときは数値計算
- 数値計算できないときには元の式をそのまま
という計算を行います。また、数値とブラ、ケット、演算子の積が出てくる場合には必ず数値を前に出すことにします:
q_compute(Exp, Result) :- Exp = A * B, !,
A0 ← A, B0 ← B, % 引数を解決しておく
( is_zero(A0) -> Result = 0 % 0 := 0 * X
; is_zero(B0) -> Result = 0 % 0 := X * 0
; is_one(A0) -> Result = B0 % X := 1 * X
; is_one(B0) -> Result = A0 % X := X * 1
; complex_mult(A0, B0, Result) % 数値として計算できれば計算
; B0 = B1 + B2 -> Result ← A0 * B1 + A0 * B2 % A0 * B1 + A0 * B2 ← A0 * (B1+ B2)
; B0 = B1 - B2 -> Result ← A0 * B1 - A0 * B2 % A0 * B1 - A0 * B2 ← A0 * (B1- B2)
; B0 = B1 * B2, is_complex_number(B1)
-> Result ← (B1 * A0) * B2, ! % 数値を優先して計算する
; B0 = B1 / B2, is_complex_number(B2) -> Result ← (A0 / B2) * B1
; A0 = A1 * A2, is_complex_number(A1)
-> Result ← A1 * (A2 * B0), !
; Result = A0 * B0 ).
q_compute(Exp, Result) :- Exp = A / B, !,
A0 ← A, B0 ← B, % 引数を解決しておく
( is_zero(A0) -> Result = 0 % 0 / X は X がなんであれ 0
; complex_div(A0, B0, Result), ! % 数値として計算できれば計算
; Result = A / B ).
0 との足し算は元のままにします：
q_compute(Exp, Result) :- Exp = A + B, !,
A0 ← A, B0 ← B, % 引数を解決しておく
( is_zero(A0) -> Result = B0 % 0 + X は X がなんであれ X
; is_zero(B0) -> Result = A0 % X + 0 は X がなんであれ X
; complex_add(A0, B0, Result), ! % 数値として計算できれば計算
; Result = A0 + B0 ).
引き算についても同様です：
q_compute(Exp, Result) :- Exp = A - B, !,
A0 ← A, B0 ← B, % 引数を解決しておく
( is_zero(A0) -> Result = -1 * B0 % 0 - X は X がなんであれ -X
; is_zero(B0) -> Result = A0 % X - 0 は X がなんであれ X
; complex_sub(A0, B0, Result), ! % 数値として計算できれば計算
; Result = A0 - B0 ).
アダマール演算で √2 を多用するので、計算できるようにしておきます。
q_compute(√ A, Result) :-
A0 ← A, % 引数を解決しておく
( is_complex_number(A0) -> complex_sqrt(A0, Result), !
; Result = √ A ).
量子ビット操作
さて、いよいよ量子ビット操作に関する定義です。
基底については数値と同じく到達したら計算終了です。
q_compute(Q, Result) :- q_basis(Q), Result = Q.
ようやく量子計算らしい定義ができるようになってきました：
q_compute( 〈↑│ │↑〉, 1).
q_compute( 〈↓│ │↓〉, 1).
q_compute( 〈↓│ │↑〉, 0).
q_compute( 〈↑│ │↓〉, 0).
q_compute( │↑〉 〈↑│, │↑〉 〈↑│).
q_compute( │↓〉 〈↑│, │↓〉 〈↑│).
q_compute( │↑〉 〈↓│, │↑〉 〈↓│).
q_compute( │↓〉 〈↓│, │↓〉 〈↓│).
数値に基底を作用させるものについては、数値×基底の掛け算に直します。
q_compute(│↑〉 X, Result) :-
X0 ← X, % 引数を解決しておく
( is_zero(X0) -> Result = 0 % 0 * X は X がなんであれ 0
; is_complex_number(X0) -> Result ← X0 * │↑〉, !
; X0 = A + B -> Result ← │↑〉 A + │↑〉 B, ! % │↑〉A+│↑〉B := │↑〉 (A+B)
; X0 = A - B -> Result ← │↑〉 A + │↑〉 B, ! % │↑〉A-│↑〉B := │↑〉 (A-B)
; X0 = C * X1, is_complex_number(C) -> Result ← C * │↑〉 X1, !
; X0 = X1 / C, is_complex_number(C) -> Result ← (1 / C) * │↑〉 X1, !
; Result ← │↑〉 X0 ).
q_compute(│↓〉 X, Result) :-
X0 ← X,
( is_zero(X0) -> Result = 0 % 0 * X は X がなんであれ 0
; is_complex_number(X0) -> Result ← X0 * │↓〉, !
; X0 = A + B -> Result ← │↓〉 A + │↓〉 B, ! % │↓〉A+│↓〉B := │↓〉 (A+B)
; X0 = A - B -> Result ← │↓〉 A + │↓〉 B, ! % │↓〉A-│↓〉B := │↓〉 (A-B)
; X0 = C * X1, is_complex_number(C) -> Result ← C * │↓〉 X1, !
; X0 = X1 / C, is_complex_number(C) -> Result ← (1 / C) * │↓〉 X1, !
; Result ← │↓〉 X0 ).
q_compute(〈↑│ X, Result) :-
X0 ← X,
( is_zero(X0) -> Result = 0 % 0 * X は X がなんであれ 0
; is_complex_number(X0) -> Result ← X0 * 〈↑│, !
; X0 = A + B -> Result ← 〈↑│ A + 〈↑│ B, ! % 〈↑│A+〈↑│B := 〈↑│ (A+B)
; X0 = A - B -> Result ← 〈↑│ A + 〈↑│ B, ! % 〈↑│A+〈↑│B := 〈↑│ (A-B)
; X0 = C * X1, is_complex_number(C) -> Result ← C * 〈↑│ X1, !
; X0 = X1 / C, is_complex_number(C) -> Result ← (1/C) * 〈↑│ X1, !
; Result ← 〈↑│ X0 ).
q_compute(〈↓│ X, Result) :-
X0 ← X,
( is_zero(X0) -> Result = 0 % 0 * X は X がなんであれ 0
; is_complex_number(X0) -> Result ← X0 * 〈↓│, !
; X0 = A + B -> Result ← 〈↓│ A + 〈↓│ B, ! % 〈↓│A+〈↓│B := 〈↓│ (A+B)
; X0 = A - B -> Result ← 〈↓│ A + 〈↓│ B, ! % 〈↓│A-〈↓│B := 〈↓│ (A-B)
; X0 = C * X1, is_complex_number(C) -> Result ← C * 〈↓│ X1, !
; X0 = X1 / C, is_complex_number(C) -> Result ← (1/C) * 〈↓│ X1, !
; Result ← 〈↓│ X0 ).
基本的な演算については以上です。
単一ビット量子ゲートの定義
ここまであれば、パウリ演算子 σx, σy, σz を以下のように定義できます：
:- op(200, fy, σx).
:- op(200, fy, σy).
:- op(200, fy, σz).
q_compute(σx Q, Result) :-
Result ← │↓〉 〈↑│ Q + │↑〉 〈↓│ Q.
q_compute(σy Q, Result) :-
Result ← i * │↓〉 〈↑│ Q - i * │↑〉 〈↓│ Q.
q_compute(σz Q, Result) :-
Result ← │↑〉 〈↑│ Q - │↓〉 〈↓│ Q.
アダマールゲートも同様に書けます：
% アダマールゲート。本当は H とかで書きたい。
q_compute(hadamard(Q), Result) :-
Result ← (1 / √2) * (
│↑〉 〈↑│ Q + │↑〉 〈↓│ Q
+ │↓〉 〈↑│ Q - │↓〉 〈↓│ Q).
これで単一量子ビットについてパウリゲートを適用できるようになりました。
?- X ← σx │↑〉, q_write(X), nl.
│↓〉
X = ('│↓〉').
?- X ← σx │↓〉, q_write(X), nl.
│↑〉
X = ('│↑〉').
?- X ← σy │↑〉, q_write(X), nl.
i*│↓〉
X = complex_value(0, 1)* ('│↓〉').
?- X ← σy │↓〉, q_write(X), nl.
-i*│↑〉
X = -1* (complex_value(0, 1)* ('│↑〉')).
?- X ← σz │↑〉, q_write(X), nl.
│↑〉
X = ('│↑〉').
?- X ← σz │↓〉, q_write(X), nl.
-│↓〉
X = -1* ('│↓〉').
?- X ← hadamard(│↑〉), q_write(X), nl.
0.7071067811865475*│↑〉+0.7071067811865475*│↓〉
X = complex_value(0.7071067811865475, 0.0)* ('│↑〉')+complex_value(0.7071067811865475, 0.0)* ('│↓〉').
?- X ← hadamard(│↓〉), q_write(X), nl.
0.7071067811865475*│↑〉-0.7071067811865475*│↓〉
X = complex_value(0.7071067811865475, 0.0)* ('│↑〉')-complex_value(0.7071067811865475, 0.0)* ('│↓〉').
感想
どうにか単一量子ビットについて計算操作を定義できました。テンソル積を '⊗' とかで定義してやれば複数量子ビットに拡張することもできそうです。
この記事は Coq で証明量子プログラミングを行う QWIRE からインスパイヤされました。
QWIRE のソースは見ているだけで楽しいので、同じようなことを Prolog でやろうと試みました。しかし、QWIRE で使われている記号が Prolog で簡単には扱えなかったために、やや不格好な感じがします。また、空白が抜けると連続したアトムとして解釈されてしまうので空白に気を使う必要があるだとか、記号の入力はやっぱり面倒（コピーペーストで頑張るのが良さそう）とか、見た目をそれらしくするのにはそれなりに苦労が必要なようです。