2016-06-14

水素原子の表現論(1.5)spectrum generating algebra

水素原子の束縛状態全体は、Lie代数so(4,2)の極小表現になっていることが知られている。このso(4,2)は、dynamical symmetryと呼ばれたり、spectrum generating algebraと呼ばれていたりする。個人的には、ハミルトニアンと可換でないものを"symmetry"と呼ぶのは紛らわしいと思うので、spectrum generating algebraという名前を使用する。spectrum generatingというのも、いまいちな名前とは思う。そして、伝説によれば、散乱状態の全体も、同様にso(4,2)の既約ユニタリー表現になっているらしいのだけど、散乱状態の方は、出典不明で、本当かどうかも、よく分からない。


#表現論的に理解したいこととして、一つには、散乱状態の全体を、何かよい基底を見つけて、so(4,2)の作用を書きたいというのがある。もう一つは、極小表現がso(4)の既約ユニタリ表現の直和に分解できるのと同様、散乱状態の方も、エネルギー一定の状態はso(3,1)の既約ユニタリ表現をなすので、散乱状態全体は、これらの表現の直積分で書けるというようなことが期待される。知る限り、そのような説明を見たことがない


束縛状態に関しては、spectrum generating algebraの具体的な実現は分かっているし、一方で、全ての束縛状態も分かっていて、極小表現は最低ウェイト表現であるから、基底状態(1s軌道解)に対して、最低ウェイト条件と特異ベクトルが消えるという条件の成立を確認すれば、水素原子の束縛状態がso(4,2)の極小表現であるということは証明できる。so(4,2)の極小表現を、so(4)に制限した時の既約分解の仕方などは、表現論の知識のみで分かる


ところで以前、

http://d.hatena.ne.jp/m-a-o/20140130#p1

で『主量子数n+1の電子殻のヒルベルト空間は、so(4)の最高ウェイト(n,0)を持ついくつかの既約表現の直和で書けることが分かったけど、各既約表現が何個出てくるかということは、表現論的に知る方法はなくて、Schrodinger方程式を解く(あるいは、Pauliの時は、まだSchrodinger方程式がなかったのでHeisenberg方程式でやっていたようだけど)しかないように思う(状態のHilbert空間を、どう取るかは物理的な理由によって決まるべきもので、表現論の立場から決める原理はないから)』と書いた


けれどもし、

(1)spectrum generating algebraの生成元が分かって

(2)束縛状態はspectrum generating algebraの何らかの既約ユニタリー表現であるということは分かっているが

(3)にも関わらず、シュレディンガー方程式の解き方は分からない

というような人がいたとしたら、束縛状態の全体について、どのくらいのことが分かるか、考えたい。現在のところ、条件(2)とか、束縛状態の全体が何故so(4,2)の極小表現となる理由を、シュレディンガー方程式を解かずに理解する方法はないと思うけど


とりあえず考えることとして、spectrum generating algebraの生成元の間には、いくつかの非自明な代数的関係式(角運動量とRunge-Lenzベクトルが直交するというような)が成立し、それらの関係式によって、U(so(4,2))の両側イデアル$J$が決定される。この両側イデアルは、極小表現上では、当然0になるはずである。$J$が消えるような、既約ユニタリー表現は、比較的限定されるだろうと期待したい。


一方、極小表現上で消えるU(so(4,2))の元全体を考えることができ、やはり両側イデアルをなす(これにはJosephイデアルという名前が付いている)。まずは、この2つのイデアルに、どれくらいの差があるか、ということを調べたい



Kepler問題に於けるspectrum generating algebraの生成元の具体的な式は、例えば、以下の論文に書いてある

SO(4,2)‐Formulation of the Symmetry Breaking in Relativistic Kepler Problems with or without Magnetic Charges

http://scitation.aip.org/content/aip/journal/jmp/12/5/10.1063/1.1665653


簡単のため、長さをBohr半径で規格化して、プランク定数も1とすると、生成元は、以下の15個

角運動量

L_1 = -y ¥frac{¥partial}{¥partial z} + z ¥frac{¥partial}{¥partial y}

L_2 = -z ¥frac{¥partial}{¥partial x} + x ¥frac{¥partial}{¥partial z}

L_3 = -x ¥frac{¥partial}{¥partial y} + y ¥frac{¥partial}{¥partial x}

scaled Runge-Lenz vector

A_1 = -¥frac{1}{2} x (¥frac{¥partial^2}{¥partial x^2} + ¥frac{¥partial^2}{¥partial y^2} + ¥frac{¥partial^2}{¥partial z^2}) + ¥frac{¥partial}{¥partial x} ¥cdot ( x¥frac{¥partial}{¥partial x} + y ¥frac{¥partial}{¥partial y} + ¥frac{¥partial}{¥partial z}) - ¥frac{1}{2} x

A_2 = -¥frac{1}{2} y (¥frac{¥partial^2}{¥partial x^2} + ¥frac{¥partial^2}{¥partial y^2} + ¥frac{¥partial^2}{¥partial z^2}) + ¥frac{¥partial}{¥partial y} ¥cdot ( x¥frac{¥partial}{¥partial x} + y ¥frac{¥partial}{¥partial y} + ¥frac{¥partial}{¥partial z}) - ¥frac{1}{2} y

A_3 = -¥frac{1}{2} z (¥frac{¥partial^2}{¥partial x^2} + ¥frac{¥partial^2}{¥partial y^2} + ¥frac{¥partial^2}{¥partial z^2}) + ¥frac{¥partial}{¥partial z} ¥cdot ( x¥frac{¥partial}{¥partial x} + y ¥frac{¥partial}{¥partial y} + ¥frac{¥partial}{¥partial z}) - ¥frac{1}{2} z


dual Runge-Lenz vector

B_1 = A_1 + x

B_2 = A_2 + y

B_3 = A_3 + z


other generators

¥Gamma_1 = r ¥frac{¥partial}{¥partial x}

¥Gamma_2 = r ¥frac{¥partial}{¥partial y}

¥Gamma_3 = r ¥frac{¥partial}{¥partial z}

T_1 = -¥frac{r}{2}(¥frac{¥partial^2}{¥partial x^2} + ¥frac{¥partial^2}{¥partial y^2} + ¥frac{¥partial^2}{¥partial z^2}) - ¥frac{r}{2}

T_2 = x ¥frac{¥partial}{¥partial x} + y ¥frac{¥partial}{¥partial y} + z ¥frac{¥partial}{¥partial z} + 1

T_3 = -¥frac{r}{2}(¥frac{¥partial^2}{¥partial x^2} + ¥frac{¥partial^2}{¥partial y^2} + ¥frac{¥partial^2}{¥partial z^2}) + ¥frac{r}{2}


そして、Chevalley生成元を以下のように決める。

e_1 = (A_1 + L_2 + i(L_1 - A_2))/2

f_1 = (A_1 - L_2 + i(L_1 + A_2))/2

h_1 = -(A_3 + i L_3)

e_2 = i(T_1+T_2+B_3-¥Gamma_3)/2

f_2 = i(T_1-T_2+B_3+¥Gamma_3)/2

h_2 = A_3+T_3

e_3 = (A_1 + L_2 -i(L_1 -A_2))/2

f_3 = (A_1 - L_2 -i(L_1+A_2))/2

h_3 = -A_3 + i L_3

末尾に、Risa/Asirでの確認用コードを載せてある。ただし、多項式係数でない微分作用素は扱いづらいので、付加的な変数$r,r'$を使って

p_x = ¥frac{¥partial}{¥partial x} + r’ x ¥frac{¥partial}{¥partial r} - r’^3 x ¥frac{¥partial}{¥partial r’}

etc.

などとして、関係式r’ = r^{-1} , r^2 = x^2+y^2+z^2によるmoduloを取ることで、計算している


長さをBohr半径で規格化しているので、1s軌道解は、規格化定数を無視すれば、単にexp(-r)である。最低ウェイト条件と特異ベクトルが消えることを確認するのは容易い。確認すべき条件は

Maxwell方程式の表現論

http://d.hatena.ne.jp/m-a-o/20150621#p1

で書いた(極小表現は、helicity0のmassless既約表現と同じ)。

e_1 ¥phi = e_3 ¥phi = 0

f_1 ¥phi = f_2 ¥phi = f_3 ¥phi = 0

h_1 ¥phi = (h_2-1)¥phi = h_3 ¥phi = 0


ここまでが準備。それで両側イデアル$J$を計算したいし、多分、機械的に計算していくことは可能だけど、直接計算は面倒そうなので、Josephイデアルの生成元を調べて、それらは勿論極小表現の上では消える(はず)だけど、微分作用素として恒等的に消えるかどうかチェックするという手順を取ろうと思う(続く)


Josephイデアルの構造については、Garfinkleによる結果が知られている

A new construction of the Joseph ideal

https://dspace.mit.edu/handle/1721.1/15620#files-area


その他、本題とは関係ない補足

(1)Josephイデアルのassociated varietyは、極小冪零軌道の閉包になっていて、極小冪零軌道は、Kepler問題の古典的相空間とみなすことができる(散乱状態は無視している)。これには、Kepler多様体という名前がついているらしい。また、表現空間の中で、generalized coherent states全体を考えることができ、Kepler多様体と同型になると思う(要確認)

(2)Coulomb problemは、(一種の拘束条件付きの)4次元調和振動子へ変換できることが知られている(この座標変換にはKustaanheimo-Stiefel変換という名前がついているらしい。この変換は水素原子の経路積分でも使われたらしい(要確認))。数学的にみると、極小表現の4変数多項式環での実現(Segal-Shale-Weilによる?)と同じもの。この実現は

Maxwell方程式の表現論(2)Gupta-Bleuler量子化の表現論的側面

http://d.hatena.ne.jp/m-a-o/20151106#p2

にも書いてある。次数$n$の多項式が(n+1)次元で、量子数n+1の電子軌道と対応する。発見の順番とか、歴史的経緯とかは知らない(この実現自体は、sp(4,R)のmetaplectic表現を、su(2,2)に制限して得られ、metaplectic表現は、物理にinspireされて生まれたと見たことがある気がする)。Kustaanheimo-Stiefel変換はある種のmagicのように私には見えるけど、表現論的には、同じ表現を二通りの方法で実現して、その間のintertwinnerを作っているのだ、と理解できる。そのような理解を得るためには、so(4)対称性だけ見ていたのでは不十分である


(3)普遍展開環$U(so(4,2))$をJosephイデアルで割った代数$hs(4,2)$は、higher spin algebraという名前で、最近(というのは主に21世紀になってから)、物理の人が調べたりしているらしい


/**** common stuffs ****/
load("noro_pd.rr");
load("noro_matrix.rr");

def assert(V){
   if(!V){error(V);}
   return 1;
}


def appdiff_rat(P , F , V){
   N = length(V)/2;
   if(type(P)==9){  /* 分散多項式表現の時 */
      DP = dp_dtop(P, V);
   }else if(type(P)<=2){
      DP = P;
   }else{
      return;
   }
   Terms = p_terms(DP , V , 2);
   for(Ret = 0 , Terms = Terms; Terms!=[] ; Terms=cdr(Terms)){
      T = car(Terms);
      if(type(F)==9){
           Tmp = dp_dtop(F , V);
      }else{
           Tmp = F;
      }
      if(T==1){
         Ret += p_nf(DP , V , V ,0)*Tmp;
      }else{
         C = p_nf(sdiv(DP-p_nf(DP,[T],V,0) , T) , V , V , 0);  /* 単項式Tの係数 */
         for(I = 0 ; I < N ; I++){
            Deg = deg(T , V[N+I]);
            for(J = 0 ; J < Deg ; J++){
               Tmp = red(diff(Tmp , V[I]));
            }
         }
         for(I = 0 ; I < N ; I++){
            Deg = deg(T , V[I]);
            Tmp *= V[I]^Deg;
         }
         Ret += C*Tmp;
      }
   }
   return red(Ret);
}


def appdiff_complex(P , F , V){
    if(type(P)==9){
        DP = dp_dtop(P,V);
    }else{
        DP = P;
    }
    RP = real(DP);
    IP = imag(DP);
    if(type(F)<=2){
        RF = real(F);
        IF = imag(F);
    }else if(type(F)==3){ /*有理式*/
        DN_conj = conj(dn(F));
        if(DN_conj != dn(F)){
           TF = red(DN_conj*nm(F)/(DN_conj*dn(F)));
        }else{
           TF = F;
        }
        RF = red(real(nm(TF))/dn(TF));
        IF = red(imag(nm(TF))/dn(TF));
    }
    Ret = red( appdiff_rat(RP , RF , V) - appdiff_rat(IP,IF,V) );
    Ret += @i*red( appdiff_rat(IP,RF,V) + appdiff_rat(RP,IF,V) );
    return Ret;
}

/* 
微分作用Pを式Fに作用させる 
Vは変数のリスト(dp_dtopなどで使われるのと同様の形式)
*/
def appdiff(P,F,V){
    if(type(P)!=6){
        if(type(F)<=3){
           return appdiff_complex(P,F,V);
        }else if(type(F)==5){   /* vector */
           Ret = newvect(length(F));
           for(I=0;I<length(F);I++){
              Ret[I] = appdiff_complex(P, F[I] , V);
           }
           return Ret;
        }
    }else{ /*行列*/
        Sz = size(P);
        N = Sz[0];
        M = Sz[1];
        if(type(F)==5){
           Ret = newvect(N);
           for (I = 0 ; I < N ; I++){
              for(J = 0 ; J < M ; J++){
                Ret[I] += appdiff_complex(P[I][J] , F[J] , V);
              }
              Ret[I] = red(real(nm(Ret[I]))/dn(Ret[I])) + @i*red(imag(nm(Ret[I]))/dn(Ret[I]));
           }
           return Ret;
        }else if(type(F)==6){
            Sz2 = size(F);
            N2 = Sz2[0];
            M2 = Sz2[1];
            Ret = newmat(N,M2);
            for (I = 0 ; I < N ; I++){
                for(J = 0 ; J < M2 ; J++){
                    for(K = 0 ; K < M ; K++){
                       Ret[I][J] += appdiff_complex(P[I][K] , F[K][J] , V);
                    }
                }
            }
            return Ret;
        }
    }
}


def dp_weyl_matmul(A,B){
    Sz1 = size(A);
    Sz2 = size(B);
    L = Sz1[0];
    M = Sz1[1];
    N = Sz2[0];
    Ret = newmat(L,N);
    for(I = 0; I < L ; I++){
       for(J = 0 ; J < N ; J++){
          for(K = 0 ; K < M ; K++){
               Ret[I][J] += dp_weyl_mul(A[I][K] , B[K][J]);
          }
       }
    }
    return Ret;
}

/*===== end of common stuffs =====*/


/* stuffs for 3D hydrogen atom */
V=[x1,x2,x3,r,ri,d1,d2,d3,dr,dri];
X1=dp_ptod(x1,V);
X2=dp_ptod(x2,V);
X3=dp_ptod(x3,V);
R=dp_ptod(r,V);
RI=dp_ptod(ri,V);
D1=dp_ptod(d1,V);
D2=dp_ptod(d2,V);
D3=dp_ptod(d3,V);
DR = dp_ptod(dr,V);
DRI = dp_ptod(dri,V);
P1 = D1 + dp_weyl_mul(RI*X1 , DR) - dp_weyl_mul(RI*RI*RI*X1 , DRI);
P2 = D2 + dp_weyl_mul(RI*X2 , DR) - dp_weyl_mul(RI*RI*RI*X2 , DRI);
P3 = D3 + dp_weyl_mul(RI*X3 , DR) - dp_weyl_mul(RI*RI*RI*X3 , DRI);
I=gr([r*ri-1,r*r-x1*x1-x2*x2-x3*x3],V,0);
Delta = dp_weyl_mul(P1,P1)+dp_weyl_mul(P2,P2)+dp_weyl_mul(P3,P3); /* Laplacian */


/* angular momentum */
L1 = -dp_weyl_mul(X2,P3) + dp_weyl_mul(X3,P2);
L2 = -dp_weyl_mul(X3,P1) + dp_weyl_mul(X1,P3);
L3 = -dp_weyl_mul(X1,P2) + dp_weyl_mul(X2,P1);

/* check [L1,L2]=L3 , [L2,L3] = L1 , [L3,L1] = L2 */
assert(dp_weyl_mul(L1,L2)-dp_weyl_mul(L2,L1)==L3);
assert(dp_weyl_mul(L2,L3)-dp_weyl_mul(L3,L2)==L1);
assert(dp_weyl_mul(L3,L1)-dp_weyl_mul(L1,L3)==L2);

T1 = (1/2)*(-R*Delta-R);
T2 = X1*P1 + X2*P2 + X3*P3 + 1;
T3 = (1/2)*(-R*Delta+R);
/* check [T1,T2]=T3 , [T2,T3]=-T1 , [T3,T1]=T2 */
/* T1,T2,T3 generates so(2,1) */
assert(p_nf(dp_dtop(dp_weyl_mul(T1,T2)-dp_weyl_mul(T2,T1)-T3,V),I,V,0)==0);
assert(p_nf(dp_dtop(dp_weyl_mul(T2,T3)-dp_weyl_mul(T3,T2)+T1,V),I,V,0)==0);
assert(p_nf(dp_dtop(dp_weyl_mul(T3,T1)-dp_weyl_mul(T1,T3)-T2,V),I,V,0)==0);



/* scaled Runge-Lenz vector */
A1 = -(1/2)*X1*Delta + dp_weyl_mul(P1,X1*P1+X2*P2+X3*P3) - X1/2;
A2 = -(1/2)*X2*Delta + dp_weyl_mul(P2,X1*P1+X2*P2+X3*P3) - X2/2;
A3 = -(1/2)*X3*Delta + dp_weyl_mul(P3,X1*P1+X2*P2+X3*P3) - X3/2;

/* dual Runge-Lenz vector */
B1 = -(1/2)*X1*Delta + dp_weyl_mul(P1,X1*P1+X2*P2+X3*P3) + X1/2;
B2 = -(1/2)*X2*Delta + dp_weyl_mul(P2,X1*P1+X2*P2+X3*P3) + X2/2;
B3 = -(1/2)*X3*Delta + dp_weyl_mul(P3,X1*P1+X2*P2+X3*P3) + X3/2;
   

assert(dp_weyl_mul(A1,A2) - dp_weyl_mul(A2,A1) == -L3);
assert(dp_weyl_mul(A2,A3) - dp_weyl_mul(A3,A2) == -L1);
assert(dp_weyl_mul(A3,A1) - dp_weyl_mul(A1,A3) == -L2);
assert(dp_weyl_mul(A1,L2) - dp_weyl_mul(L2,A1)==A3);
assert(dp_weyl_mul(@i*A1+L1,@i*A2+L2)-dp_weyl_mul(@i*A2+L2,@i*A1+L1)==2*(@i*A3+L3));


assert(dp_weyl_mul(B1,B2)-dp_weyl_mul(B2,B1)==L3);
assert(dp_weyl_mul(B2,B3)-dp_weyl_mul(B3,B2)==L1);
assert(dp_weyl_mul(B3,B1)-dp_weyl_mul(B1,B3)==L2);
assert(dp_weyl_mul(B1,L2) - dp_weyl_mul(L2,B1)==B3);
assert(dp_weyl_mul(B2,L3) - dp_weyl_mul(L3,B2)==B1);
assert(dp_weyl_mul(B3,L1) - dp_weyl_mul(L1,B3)==B2);
assert(dp_weyl_mul(B1+L1,B2+L2)-dp_weyl_mul(B2+L2,B1+L1)==2*(B3+L3));


assert(dp_weyl_mul(A1,B1)-dp_weyl_mul(B1,A1)==T2);
assert(dp_weyl_mul(A2,B2)-dp_weyl_mul(B2,A2)==T2);
assert(dp_weyl_mul(A3,B3)-dp_weyl_mul(B3,A3)==T2);
assert(dp_weyl_mul(T2,A1)-dp_weyl_mul(A1,T2)+B1==0);
assert(dp_weyl_mul(T2,A2)-dp_weyl_mul(A2,T2)+B2==0);
assert(dp_weyl_mul(T2,A3)-dp_weyl_mul(A3,T2)+B3==0);
assert(dp_weyl_mul(T2,B3)-dp_weyl_mul(B3,T2)+A3==0);


G1 = R*P1;
G2 = R*P2;
G3 = R*P3;
assert(p_nf(dp_dtop(dp_weyl_mul(G1,T1)-dp_weyl_mul(T1,G1)-A1,V),I,V,0)==0);
assert(p_nf(dp_dtop(dp_weyl_mul(G2,T1)-dp_weyl_mul(T1,G2)-A2,V),I,V,0)==0);
assert(p_nf(dp_dtop(dp_weyl_mul(G3,T1)-dp_weyl_mul(T1,G3)-A3,V),I,V,0)==0);
assert(p_nf(dp_dtop(dp_weyl_mul(G1,T3)-dp_weyl_mul(T3,G1)-B1,V),I,V,0)==0);
assert(p_nf(dp_dtop(dp_weyl_mul(T1,A2)-dp_weyl_mul(A2,T1)-G2,V),I,V,0)==0);
assert(p_nf(dp_dtop(dp_weyl_mul(T1,A3)-dp_weyl_mul(A3,T1)-G3,V),I,V,0)==0);
assert(p_nf(dp_dtop(dp_weyl_mul(T1,B1)-dp_weyl_mul(B1,T1),V),I,V,0)==0);
assert(p_nf(dp_dtop(dp_weyl_mul(T1,B2)-dp_weyl_mul(B2,T1),V),I,V,0)==0);
assert(p_nf(dp_dtop(dp_weyl_mul(T1,B3)-dp_weyl_mul(B3,T1),V),I,V,0)==0);
assert(p_nf(dp_dtop(dp_weyl_mul(T3,A1)-dp_weyl_mul(A1,T3),V),I,V,0)==0);
assert(p_nf(dp_dtop(dp_weyl_mul(T3,A2)-dp_weyl_mul(A2,T3),V),I,V,0)==0);
assert(p_nf(dp_dtop(dp_weyl_mul(T3,B1)-dp_weyl_mul(B1,T3)+G1,V),I,V,0)==0)
assert(p_nf(dp_dtop(dp_weyl_mul(T3,B2)-dp_weyl_mul(B2,T3)+G2,V),I,V,0)==0);
assert(p_nf(dp_dtop(dp_weyl_mul(T3,B3)-dp_weyl_mul(B3,T3)+G3,V),I,V,0)==0);
assert(p_nf(dp_dtop(dp_weyl_mul(T3,G3)-dp_weyl_mul(G3,T3)+B3,V),I,V,0)==0);
assert(p_nf(dp_dtop(dp_weyl_mul(A3,G3)-dp_weyl_mul(G3,A3)+T1,V),I,V,0)==0);


assert(p_nf(dp_dtop(dp_weyl_mul(A3,T3)-dp_weyl_mul(T3,A3),V),I,V,0)==0);
assert(p_nf(dp_dtop(dp_weyl_mul(L3,T3)-dp_weyl_mul(T3,L3),V),I,V,0)==0);
assert(p_nf(dp_dtop(dp_weyl_mul(L3,A3)-dp_weyl_mul(A3,L3),V),I,V,0)==0);



/*
$T_1 (e^{-r}) = (1-r)e^{-r}$
$T_2 (e^{-r}) = (1-r)e^{-r}$
$T_3 (e^{-r}) = e^{-r}$

$A_1 (e^{-r}) = 0$
$A_2 (e^{-r}) = 0$
$A_3 (e^{-r}) = 0$

$L_1 (e^{-r}) = 0$
$L_2 (e^{-r}) = 0$
$L_3 (e^{-r}) = 0$

$\phi$: lowest weight vector for minrep of so(4,2)

$f_1(\phi) = f_2(\phi) = f_3(\phi) = 0$
$h_1 \phi = (h_2-1)\phi = h_3 \phi = 0$
$e_1 \phi = e_3 \phi = 0$

$h_4 = h_1+h_2 , h_5 = h_2+h_3 , h_6 = h_1+h_2+h_3$
*/

/* Chevalley generators */
E1 = (A1+L2+@i*L1-@i*A2)/2;
F1 = (A1-L2 +@i*L1+@i*A2)/2;
H1 = -(A3+@i*L3);

E2 = @i*(T1+T2+B3-G3)/2;
F2 = @i*(T1-T2+B3+G3)/2;
H2 = A3+T3;

E3 = -(@i*(@i*A1+L1)-(@i*A2+L2))/2;
F3 = (A1-L2-@i*L1-@i*A2)/2;
H3 = -A3+@i*L3;

assert(dp_weyl_mul(E1,F1)-dp_weyl_mul(F1,E1)==H1);
assert(dp_weyl_mul(H1,E1)-dp_weyl_mul(E1,H1)==2*E1);
assert(dp_weyl_mul(H1,F1)-dp_weyl_mul(F1,H1)==-2*F1);

assert(p_nf(@i*dp_dtop(dp_weyl_mul(H2,E2)-dp_weyl_mul(E2,H2)-2*E2,V),I,V,0)==0);
assert(p_nf(dp_dtop(dp_weyl_mul(E2,F2)-dp_weyl_mul(F2,E2)-H2,V),I,V,0)==0);
assert(p_nf(@i*dp_dtop(dp_weyl_mul(H2,F2)-dp_weyl_mul(F2,H2)+2*F2,V),I,V,0)==0);

assert(dp_weyl_mul(E3,F3)-dp_weyl_mul(F3,E3)==H3);
assert(dp_weyl_mul(H3,E3)-dp_weyl_mul(E3,H3)==2*E3);
assert(dp_weyl_mul(H3,F3)-dp_weyl_mul(F3,H3)==-2*F3);

assert(dp_weyl_mul(H1,E3)-dp_weyl_mul(E3,H1)==0);
assert(dp_weyl_mul(H3,E1)-dp_weyl_mul(E1,H3)==0);
assert(dp_weyl_mul(E1,F3)-dp_weyl_mul(F3,E1)==0);
assert(dp_weyl_mul(E3,F1)-dp_weyl_mul(F1,E3)==0);

/* [H2 , E1] = -E1 */
assert(p_nf(dp_dtop(dp_weyl_mul(A3+T3,A1+L2)-dp_weyl_mul(A1+L2,A3+T3)+(A1+L2),V),I,V,0)==0);
assert(p_nf(dp_dtop(dp_weyl_mul(A3+T3,L1-A2)-dp_weyl_mul(L1-A2,A3+T3)+(L1-A2),V),I,V,0)==0);

/* [H2 , F1] = F1 */
assert(p_nf(dp_dtop(dp_weyl_mul(A3+T3,L1+A2)-dp_weyl_mul(L1+A2,A3+T3)-(L1+A2),V),I,V,0)==0);
assert(p_nf(dp_dtop(dp_weyl_mul(A3+T3,A1-L2)-dp_weyl_mul(A1-L2,A3+T3)-(A1-L2),V),I,V,0)==0);

/* [H1,E2] = -E2 */
assert(p_nf(dp_dtop(@i*(dp_weyl_mul(H1,E2)-dp_weyl_mul(E2,H1)+E2),V),I,V,0)==0);

/* [H1,F2] = F2 */
assert(p_nf(dp_dtop(@i*(dp_weyl_mul(H1,F2)-dp_weyl_mul(F2,H1)-F2),V),I,V,0)==0);



/* some vanishing conditions */ 
assert(dp_weyl_mul(L1,A1)+dp_weyl_mul(L2,A2)+dp_weyl_mul(L3,A3)==0);

/* A1*A1+A2*A2+A3*A3-L1*L1-L2*L2-L3*L3==T3*T3-1 */
assert(p_nf(dp_dtop(dp_weyl_mul(A1,A1)+dp_weyl_mul(A2,A2)+dp_weyl_mul(A3,A3)-dp_weyl_mul(L1,L1)-dp_weyl_mul(L2,L2)-dp_weyl_mul(L3,L3)-dp_weyl_mul(T3,T3)+1,V),I,V,0)==0);

/* B1*B1+B2*B2+B3*B3+L1*L1+L2*L2+L3*L3==T1*T1+1 */
assert(p_nf(dp_dtop(dp_weyl_mul(B1,B1)+dp_weyl_mul(B2,B2)+dp_weyl_mul(B3,B3)+dp_weyl_mul(L1,L1)+dp_weyl_mul(L2,L2)+dp_weyl_mul(L3,L3)-dp_weyl_mul(T1,T1)-1,V),I,V,0)==0);


三角形の合同類のなす空間と有限変形理論

三角形の合同条件と不変式論

http://d.hatena.ne.jp/m-a-o/20130104#p1

みたいなことを昔書いた。ここで書いた考え方だと、三角形が退化した場合(三点が一直線上にある場合とか、一点に集中している場合)の扱いが面倒で、そのへんのことは、詳しく書いてない。が、最近、三角形の合同類のなす空間は、もっと単純な記述ができることに気付いた。


答えとしては、退化してない三角形の合同類の空間は$GL_{+}(2,\mathbf{R})/SO(2)$になる。実際、三角形の一点を原点に持ってきて、残りの二点の座標が、(a,b),(c,d)とすると、非退化な三角形であれば、ad-bc≠0が成立するから、そのような2点の配置の空間は$GL(2,\mathbf{R})$そのものである。一点を原点に置いてるので、これで並進の自由度は消えていて、残りは回転(と鏡映)の自由度で、商空間を考えれば、退化してない三角形の合同類の空間を与える。一般の次元でも、$GL_{+}(n,\mathbf{R})$(行列式が正のn行n列の行列全体)が$n$-単体の合同類の空間に推移的に作用して、固定部分群が$SO(n)$になるので、非退化な$n$-単体の合同類の空間は、$GL_{+}(n,\mathbf{R})/SO(n)$になる。この空間は、正定値対称行列全体と同一視できる($GL_{+}(n,\mathbf{R})$の元が、正定値対称行列と直交行列の積に一意に分解できるというのが、極分解定理)


#$SL(2,\mathbf{R})/SO(2)$は有名なポアンカレの上半平面になる。$GL_{+}(1,\mathbf{R})$は定数倍であるから、単なる拡大・縮小であり、$SO(n)$と合わせて相似な三角形全体を通る。$GL_{+}(2,\mathbf{R})$は$GL_{+}(1,\mathbf{R})$と$SL(2,\mathbf{R})/SO(2)$の直積であるから、$SL(2,\mathbf{R})/SO(2)$は三角形の相似類の空間である。$SL(n,\mathbf{R})/SO(n)$も対称空間の代表的な例である。数学的には、合同類より相似類の方が自然なのかもしれない


並進の自由度をad hocに消しているけど、$n$-単体の一様並進も入れると、$n$-単体の合同類の空間に、アフィン変換群が推移的に作用し、その場合、固定部分群は、合同変換群$Iso(n,\mathbf{R})$となる。回転と鏡映を考えると$GL(n,\mathbf{R})$が作用し、固定部分群は$O(n)$になる。


これを考えたのは、有限要素法の大変形解析のことを考えていた時。有限要素法では物体の変形がアフィン変換(並進の自由度を消せば一般線形変換)で線形補間すると十分よい近似になる微小領域を使う。この微小領域の初期状態(平衡状態)として、三角形や四面体を取ると、線形変換やアフィン変換によって、任意の三角形や四面体が得られる。そして、微小領域の形状のみを考えるなら、合同変換の自由度(剛体並進・剛体回転)は邪魔である。これは、$n$-単体の合同類の空間を考えた時と同じ考え方になっている。結局、(局所的な)有限歪みは、$GL(n,\mathbf{R})/O(n)$で分類され、超弾性体などの場合、弾性エネルギー/ひずみエネルギーは、この空間上の関数とみなすことができる。次元は、n^2-n(n-1)/2=n(n+1)/2で、n=2なら3、n=3なら6で、これは三角形や四面体の辺の数と等しい


※)三角形や四面体の定ひずみ要素を考える場合、平衡状態(要素の歪みエネルギーが最小となる配置)に於ける節点の座標と、変形後の節点の座標から、前者を後者へ変換するアフィン変換が一意に定まり、これの線形変換部分が変形勾配テンソル(テンソルと名前が付いてるけど、単なる行列)となる。変形勾配テンソルは直接、極分解するか、あるいは、転置行列との積を取ることにより、剛体回転成分が除かれ、有限歪みを得ることができる。応力と歪みの構成則があれば、応力成分と歪み成分の積を足し合わせて、弾性エネルギー密度を得られる。こうして得られた弾性エネルギーは、変形後の節点座標の関数と見なすことができるので、変形後の各節点の座標成分で微分すると、節点力が定まる


連続体が非圧縮性条件を満たす場合、局所的には体積が保存するので、$GL_{+}(n,\mathbf{R})$の代わりに$SL(n,\mathbf{R})$を考えることになる。従って、非圧縮性有限歪みの分類空間は、対称空間$SL(n,\mathbf{R})/SO(n)$で記述される(という事実が、何かの役に立つわけではない)

日本語/英語の組み合わせ範疇文法パーサーを作った

toyccg

https://github.com/vertexoperator/toyccg


英語については、CCG(Combinatory Categorial Grammar)のparserは公開されているものが、いくつか存在する(疑問文とか命令文は、あんまり対応してなかったりする。toyccgも疑問文や命令文は、あんまり対応できてない。やればできるとは思う)けど、日本語については、知る限りではない(公開せずに実装している人はいると思うけど)ので。原理的に新しいことは何もないけど、実用的に使うのに、必要そうな機能は、一通り実装してあると思う


自然言語用のパーサーとしては

・pythonで1500〜2000行程度((自然)言語非依存のコア部分は、1000行くらい。外部ライブラリは使ってないし他(プログラミング)言語への移植は楽と思う)で比較的小さい

・日本語と英語同時対応

・機械学習フリー

などの特徴がある。機械学習を使用してないのは、そもそも、日本語の学習用データが存在しない(多分)(英語用では、CCGBankがあって、普通これを使うと思うけど、どうやって入手するのか謎なので、やっぱり使えない)のもあるし、機械学習頼りだと、手で修正するということがやりづらくなるという面もある。手でルール作るとか、時代の流れに逆行している気もするけど、実用上は悪いことばかりでもない(CCGBankはTreeBankからの自動生成っぽいけど、学習用データを自分で作る場合、やっぱり手間がかかるわけだし)

#CFGを確率化した確率文脈自由文法(PCFG)と同様、CCGを確率化した、確率的組み合わせ範疇文法(PCCG)のような方法も存在する


nltkのCCGパーサーとの違いは

・デフォルトで辞書を用意している

・normal form parsingをサポートしている

・nltkのパーサーで日本語を解析する場合、単語単位への分かち書きは別に行う必要があるけど、toyccgは、分かち書きと構文解析を同時にやる

というあたり。分かち書き対応したおかげで、英語でも、New YorkとかWhite Houseみたいなのを一単語のようにして扱うことができるようになった(やろうと思えば、"This is a pen"のような文単位で統語範疇を辞書に登録することもできる)。あとまぁ、長単位・短単位みたいなのを気にする必要も、あまりない


#"数学的構造"みたいなのを辞書登録すると、数学+的+構造で分解するのと、"数学的構造"で単語扱いする結果が、ダブって出現することになる。こういう単語が多い文では、計算時間の増大にも繋がるし、本質的に同じ解析結果が沢山出現することにもなる。toyccgで、このパターンで計算量が増えるケースは結構ある


#CCGでは、長単位・短単位以外に、品詞についても、考える必要がない(これをメリットと捉えるかは人によるかもしれないけど)。名詞・動詞は分かりやすいけど、実際に作業していて、日本語文法の知識が殆どない私には、助詞・助動詞の区別は迷うことがあるし、例えば、もう死語のような気がするけど

「なう」の品詞ってなんですか?

http://togetter.com/li/258854

という話でも、CCGだと、"S\N"みたいな統語範疇を設定してやれば"渋谷なう"みたいな文章を解析できる。これに類似した統語範疇(つまり、"S[*]\N[*]")を持つ単語は、今のところ他にはない。動詞や形容詞は統語範疇として"S\NP"を持つので、近いと言えるかもしれない。"ご飯食べてるなう"だと、"S[now]\S"とするとか。統語素性[now]を付けているのは、"ご飯食べてるなうなう"みたいな文章を是としないため。



統語素性の設計について。日本語の統語素性は、ベースになるようなものがないので、勝手に作ったものが沢山ある(一応、「日本語文法の形式理論」という本が一冊出ていて、それは参考にした)。英語の場合は、CCGBankがデファクトスタンダードになっているので、なるべくそれに従っているけれど、改変してある部分もある(例えば、toyccgは前置詞句PPを使わない)。これらは、辞書の問題なので、辞書を差しかえれば、実装の方を変えることなく、他の統語範疇・統語素性設計に対応できる(はず)

#CCGBankは、結構謎のunary rulesがあって、本来のCCGからは、はみ出しているので、これが気に入らない。toyccgも、このへんの問題はクリアできてなくて、適切に統語範疇・統語素性を設計すれば解決するのか、それともどうしようもないのかは、不明。英語ではCCGBankが標準のようになっているけど、本当に、これでいいのかはよく検討するべきことと個人的には思う


unary rulesとして、どういうものがあるかは、

https://github.com/mikelewis0/easyccg/blob/e42d58e08eb2a86593d52f730c5afe222e939781/training/unaryRules

などを見ると分かる



機械学習ベースじゃないことのデメリットとして、構文解析できないケースが、増えてしまう(実用上は、何でもいいから、とりあえず結果を出してくれた方が嬉しい時もあれば、ダメな時には、あからさまに失敗してくれた方がいい時もあると思うので、一概にデメリットとも言い切れない。toyccgの場合、日本語として、文法が崩壊している文章は、解析されない可能性が高い)。一応、それなりには頑張っていて

>>> import toyccg.japanese as jpn
>>> jpn.run(u"すももももももももの内。")
test run : sentence=すももももももももの内。
すもも	N
も	(NP[sbj]/NP[nom-enum])\N
もも	N
も	NP[nom-enum]\N
もも	N
の	((S[nom]\NP[sbj])/(S[nom]\NP[sbj]))\N
内	S[nom]\NP[sbj]
。	ROOT\S[nom]

>>> jpn.run(u"象は鼻が長い。")
test run : sentence=象は鼻が長い。
象	N
は	NP[sbj]\N
鼻	N
が	NP[ga-acc]\N
長い	(S\NP[sbj])\NP[ga-acc]
。	ROOT\S

くらいの解析はできる。sentences.ja.txtには、解析できる文章として、どれくらいのものがあるか、100以上の例がある。一応、ニュースや論文に出てくるような、比較的固めの文章なら、そこそこ解析できるように頑張っている(解析できるとは言っていない)。twitterや話し言葉で出てくるような文章となると、なかなか難しい。これらの文章では、格助詞の省略がよく起きて、難しい原因の一つとなる。その他、未知語が含まれている場合(これは、ニュースや論文でも、よくある。一応、文字種に基づく未知語推定を行えるしているが)や、倒置法とかは全く対応されてない


助詞が省略されていても、

>>> jpn.run(u"お腹すいた。")
test run : sentence=お腹すいた。

>>> jpn.parser.lexicon.setdefault(u"お腹",[]).append( "NP[sbj]" )
>>> jpn.run(u"お腹すいた。")
test run : sentence=お腹すいた。
お腹	NP[sbj]
すい	IV[euph]
た	(S\NP[sbj])\IV[euph]
。	ROOT\S

みたいに、適切な統語範疇を辞書登録してやれば、対応できないわけではないけど。安易にこういうことをして、間違った構文解析結果が増えてしまうと、それはそれで困る。


#既存部分への影響を、頭を使わずに、確実に抑えたいなら、NP[x]とか適当な新しい素性を導入して、動詞や形容詞のような単語にも、対応する統語範疇("S\NP[x]"など)を追加していくなどの手もある。スマートな解決策とは言えないけど、実用的な解析を行うには、そういうことも必要かもしれない。日本語は、話し言葉レベルでは、「私、昨日、学校行った」みたいな格助詞が全部落ちたような文章も可能であるから、そこまで行くと、現時点では、どれが主語か正しく判定するのは難しい



jpn.runは、形態素解析的な結果しか出さないけど、内部的には、構文木が作られている

>>> import toyccg.japanese as jpn
>>> r = jpn.parser.parse(u"あらゆる現実を、全て自分の方へ捻じ曲げたのだ。")
>>> t = r.next()
>>> print(t.show())
(LApp (LApp (SkipCommaJP (LApp (RApp [あらゆる:N/N] [現実:N]) [を:NP[obj]\N]) [、:COMMA]) (RBx [全て:S[null]/S[null]] (RBx (LApp (RApp (LApp [自分:N] [の:(N/N[base])\N]) [方:N[base]]) [へ:(S[null]/S[null])\N]) (LB (LApp [捻じ曲げ:TV[euph]] [た:(S[null]\NP[obj])\TV[euph]]) [のだ:S[null]\S[null]])))) [。:ROOT\S[null]])

あんまり見やすい出力ではないけど。多くの場合、構文解析結果は一意でないので、結果はジェネレータで、一つずつ返ってくる。

#デフォルトでは、ソースコード中で、shortcutと書いてある"近似"が有効になっていて、例えば、"かわいい私の妹"は、二種類の係り受け構造("かわいい"が"私"に係る場合と、"妹"に係る場合)が存在するが、どちらの場合も統語範疇は"N"となる(これは、本質的に構造が違うので、normal form parsingでは対処できない)ので、それ以降の解析には影響しない。こういうのが沢山あると、長文では、解析結果が沢山出てきて、計算時間が膨大になる。こういう状況を回避するため、先に出てきた結果だけを取って、それ以降のは使わないという処理をしている。最終的な構文木が一つしか必要ない場合は、shortcutをTrueにしておく方が、計算が速く済む(大体、ニュースとか論文とかの文章だと、これがないと一文の計算が全然終わらないことが多い)けど、本当に、全ての結果が必要な場合は、shortcutをFalseにする必要がある


分かち書きしたいだけなら、

>>> print [c.token for c in t.leaves()]
[u'\u3042\u3089\u3086\u308b', u'\u73fe\u5b9f', u'\u3092', u'\u3001', u'\u5168\u3066', u'\u81ea\u5206', u'\u306e', u'\u65b9', u'\u3078', u'\u637b\u3058\u66f2\u3052', u'\u305f', u'\u306e', u'\u3060', u'\u3002']

とかできる



感想。色々と実験するプロトタイプ的なつもりだったのだけど、思ったよりも、よく解析できるな、という印象(自画自賛)。


課題

・複合名詞対応

・格助詞落ち対応


今後。もう少し改善しようかと思ってたけど、飽きてしまったので、当分いじることはなさげ(完)

2016-01-03

非可換Groebner基底と表現論(1)概要メモ

Groebner基底は、多項式環のイデアルの"よい"生成元のことで、Groebner基底を計算するBuchbergerアルゴリズムなどの存在によって、計算機代数で広く使われている。Groebner基底を非可換多項式環に拡張しようというのは、誰でも考えるようなことではあると思うけど、1986年に、Moraという人が、非可換多項式環の両側イデアルに対する非可換Groebner基底を計算するように、Buchbergerアルゴリズムを拡張したらしい(これは、停止するとは限らない半アルゴリズム)。


ところで、非可換代数では、両側イデアルと左イデアル/右イデアルの区別が生じる。Lie環の表現論の文脈では、普遍展開環の定義関係式は、非可換多項式環の両側イデアルを与え、普遍展開環の左イデアルによる商加群を考えると、Lie環の表現が得られる(抽象論としては、unitalな非可換結合代数$R$の単純加群は、ある極大左イデアル$I$によって、$R/I$の形で書ける)。そういうわけで、表現論の主役は、むしろ左イデアルであるから、両側イデアルの非可換Groebner基底が分かっても、特別嬉しいことはない。


などと思っていたところ

GROEBNER-SHIRSHOV BASES FOR REPRESENTATION THEORY

http://www.math.uconn.edu/~khlee/Papers/GS-pair.pdf

という論文を読んだ。もしかしたら、これ以前に同じことを考えた人はいるかもしれないけど、確認はしてない(求情報)。タイトルにあるGroebner-Shirshov basisというのは、非可換Groebner基底と同じものだ(と思う)けど、上記論文では、Groebner-Shirshov pairという両側イデアルの"よい"生成元と左イデアルの"よい"生成元の組を定義している。論文では、Lie環sl(3)の有限次元既約表現について、Groebner-Shirshov pairを与えている(全ての有限次元既約表現に対する一般形を知るため、証明を行っているけど、Buchbergerのアルゴリズムの拡張として、Groebner-Shirshov pairを計算する半アルゴリズムも与えている)


#その後の論文で、Lie環sl(n)の有限次元既約表現についても具体形を与えている

Grobner-Shirshov Bases for Irreducible sl_{n + 1}-Modules

http://www.sciencedirect.com/science/article/pii/S0021869300983819

http://citeseerx.ist.psu.edu/viewdoc/summary?doi=10.1.1.514.1152



Groebner-Shirshov pairが分かると、何ができるのかというと、イデアルメンバーシップ問題を解いたり、イデアルが等しいかどうかの判定を行ったりできるはずだけど、それは表現論的には嬉しい場面がなさそうな気がする。応用の一つとして、標準単項式の計算がある。これは、Groebner基底の時にもあったものであるけど、標準単項式は表現の基底を与える(Groebner"基底"とかGroebner-Shirshov"基底"は、生成元であるけど、こっちは表現空間の基底を指す)。ある代数の表現論で、具体的な計算を実行しようと思うと、表現の"よい"基底を取るということは、重要な問題となる。この種の問題に対する答えとして、Gelfand-Tsetlin基底や結晶基底があるけども、その存在や構成は、代数や表現の性質に強く依存している。標準単項式は、より一般的な代数に対して定義可能なものとなる


#Groebner-Shirshov pairの前に、左イデアルの生成元が分かっている必要はあるわけだけど、例えば、複素半単純Lie環の有限次元既約表現の場合、これは既約最高ウェイト表現であり、特異ベクトルがどこに出るかも全部分かるので、左イデアルの生成元は分かる。既約表現の構成の仕方によっては、左イデアルの生成元を調べること自体が難しい場合もあるとは思うし、最高ウェイト表現でない場合、一般にどれくらいのことが知られているのかよく分からない(Casimir作用素=普遍展開環の中心の元は、既約表現に定数倍で作用するので、それによって、左イデアルの生成元がいくつか得られる場合もある。それ以外に、どれくらいの付加的な条件があるのかは、分からない)



とりあえず、証明とかは抜きで、ざっと流れを把握する。可換Groebner基底の時と殆ど同じではあるけど、detailは、非可換性のために多少複雑になっている


まず単項式順序。非可換の時は、length-lexicographic order(長いのでdeglex順序と呼ぶことにする)というもの以外使われてるのを見ない(これしか適切な単項式順序がないのかどうかは知らない。要検討)。ちゃんとした定義は上論文参照。2変数x,yでx<yの場合、1 < x < y < x^2 < xy < yx < y^2 < x^3 < x^2y < xyx < xy^2 < yx^2 < ....などとなる。それで、可換の時と同様、leading monomialとかleading coefficientが定まる。leading coefficientが1の時、monicであるという


normal formもしくは剰余。可換の時と同様、剰余によって、normal formを計算/定義する。これが一意に定まるのが、非可換Groebner基底やGroebner-Shirshov pairで重要な性質(定義としてもいい)。$S$と$T$をそれぞれ、monicな非可換多項式の有限集合として($S$を両側イデアルの生成元、$T$を左イデアルの生成元と思う)、非可換多項式$p$が

p = ¥sum ¥alpha_i (a_i s_i b_i) + r

と書けるとする。但し

(1)¥alpha_iは定数(複素数とか有理数)

(2)a_i,b_iは任意の単項式、s_i ¥in Sで、a_i LT(s_i) b_iLT(p)を超えない($LT(f)$は$f$のleading monomial)

(3)$r$の任意の単項式は、任意の$S$の元$s$について$a LT(s) b$と書けないとする

このような表示は、剰余の手続きによって得ることができる($a_i$,$b_i$は単項式なので、和は有限和であれば、異なる$i,j$について$s_i=s_j$となってもよい)。この時$r$を$p$の$S$に関するnormal formとする。$S$に関するnormal formが任意の非可換多項式$p$について一意に定まる時、$S$の元は、非可換Groebner基底をなすという。


同様に$(S,T)$に関するnormal formを定める。Groebner-Shirshov pairの場合、両側剰余と左剰余が混ざることに注意する必要がある。非可換多項式$p$が

p = ¥sum ¥alpha_i (a_i s_i b_i) + ¥sum ¥beta_j c_j t_j + r

と書けるとする。但し

(1)¥alpha_i,¥beta_jは定数

(2)a_i,b_i,c_jは任意の単項式、s_i ¥in S , t_j ¥in T

(3)$r$の任意の単項式は、任意の$S$の元$s$についてa LT(s) bと書けず、任意の$T$の元$t$についてc LT(t)とも書けないとする

この時$r$を$p$の$(S,T)$に関するnormal formとする。$(S,T)$に関するnormal formが任意の非可換多項式$p$について一意に定まる時、$(S,T)$はGroebner-Shirshov pairであるという。


ちゃんと書くと式は何か複雑であるけど、$S$は両側イデアルの生成元であるので、両側剰余を取っている。$T$は($S$が生成する両側イデアルによる剰余代数の)左イデアルの生成元となるものなので、左側からの掛け算しか許していない


非可換Groebner基底もGroebner-Shirshov pairも、Buchbergerアルゴリズムを拡張した完備化手続きによって計算できる。注意することとして、非可換Groebner基底の場合、可換の場合のS-polynomialに相当するものが一つとは限らないという事情がある。例えば、f = xy - x , g = yx - yでx<yとするdeglex順序が入っているとすると、fとgのleading monominalはそれぞれxyとyxになる。可換の時、これらが消えるようにLCMを選んでいたけども、fx - xg = -x^2 + xyで(-x^2+xy=-x^2+(f+x)なので-x^2+xにreductionされることに注意)、一方、yf - gy = -yx + y^2でもある(-yx+y^2=-(g+y)+y^2なのでy^2-yにreductionされる)。こうして、次は4つの多項式{f,g,x^2-x,y^2-y}の各ペアに対して、同様の計算を繰り返す(この例では、これで終わりとなる)。ついでに停止性も保証されない


$(S,T)$のGroebner-Shirshov pairを計算する場合は、まず$S$の非可換Groebner基底を計算して、その後、$T$の方を完備化していくことになる。$T$の方の完備化は、少し込み入っていて、2ステップの繰り返しとなる。まず、$T$同士の元で、right-justified composition(これも、S-polynomialの計算の一般化)というものを計算して簡約、新しい元の追加を行う。次に、完備化した$S$の元とupdateした$T$の元の対で、$T$に元を追加していく。これを繰り返す(詳細は論文参照。言葉だけで説明するのが辛くなってきた)。


アルゴリズムの正しさは、非可換Groebner基底/Groebner-Shirshv pairでも、剰余の一意性と"closed under composition"という性質(可換の場合の、全てのS-polynomialを簡約すると0になるという性質に対応)が同値であることを示せば自明となる


例。Lie環$sl(3)$の最高ウェイト¥Lambda_1+¥Lambda_2の既約表現を考える。$sl(3)$の下三角部分しか重要でないので、Chevalley生成元$f_1,f_2$のみを考える。単項式順序は$f_1 < f_2$で定まるdeglex順序とする。これらには2つのSerre関係式がある

p = f_2^2f_1 - 2 f_2f_1f_2 + f_1f_2^2

q = f_2f_1^2 - 2 f_1f_2f_1 + f_1^2f_2

が両側イデアルの方の生成元。唯一の"S-polynomial"p f_1 - f_2 qが存在するけども、$\{p,q\}$による剰余を計算すると0になるので、$S=\{p,q\}$は非可換Groebner基底をなしている。


一方、片側イデアルの方は、2つの非可換多項式

T=¥{t_1 = f_1^2 , t_2 = f_2^2¥}

で生成される。$T$の完備化の第一段階であるright-justified compositionは、今回することがない。次に、$S$と$T$の間の完備化計算を行うと、次の3つのS-polynomialが出る($S$と$T$による剰余も計算している)

 t_3 = q - f_2 t_1 = -2 f_1 f_2 f_1 + f_1^2 f_2

 t_4 = p f_1 - f_2^2 t_1 = -2 f_2 f_1 f_2 f_1 + 2 f_1 f_2 f _1 f_2 - f_1^2 t_2 ¥equiv -2 f_2 f_1 f_2 f_1 + 2 f_1 f_2 f _1 f_2

 t_5 = q f_1 - f_2 f_1 t_1 = -3 f_1 f_2 t_1/2 - f_1 q/2 + f_1^3 f_2/2 ¥equiv f_1^3 f_2/2

これら3つの多項式+元の$T$の2つの多項式について、"right-justified completion"を行う。のだけど、その前に、2つ目の多項式は、

 f_2 t_3 - q f_2 + f_1^2 t_2

なので、除いてしまっても問題ない。ということを繰り返すのだけど、手計算では、予想より大分辛かったので以下略


というわけで、何はともあれ、多くのことが可換Groebner基底の時とパラレルに進む


TODO:あとで、Groebner-Shirshov pairを計算するコードを追加する

2015-11-06

Combinatory Categorical Grammar

というものを最近調べた。前身として、Categorical Grammar(範疇文法)というものがあったらしい。


歴史的には、範疇文法は文脈自由文法(CFG)より古いっぽいが、範疇文法は、文脈自由文法と同等の計算能力しかないことが知られているらしい。基本的には、

S -> NP VP

のような生成規則があったら、VP=S\NPないしNP=S/VPとして、前者であれば、VPは、"左に型NPの引数を受け取って、型Sの値を返す関数"と見なし、後者なら、NPは、"右に型VPの引数を受け取って型Sの値を返す関数"と見なす。一般的なプログラム言語で、普通の中置きされた(int型同士の)"+"演算子であれば、(int\int)/intのような型になる感じだと理解すると、分かりやすい気がする。


CFGと比較した場合

・CFGと比べるとラムダ計算で意味付けする方法が明確なので、意味論に繋げやすい

・各単語ごとに、"型"が付いて、統語規則の方は変更しない語彙化文法というものになっている(実用上のメリットとして、CFGで文法ルールを作って、生成規則を増やす場合、影響範囲を予測するのが難しいけど、語彙化文法なら、追加した語彙が出てこない文には影響しない)

というメリットがある。CFGから範疇文法への変換は、一回標準形にする必要があり、例えば、二重他動詞DTVは、CFGでは生成規則VP->DTV NP NPを持つと思われるけど、直接対応する範疇文法は作れない(範疇文法なら、DTV=(VP/NP)/NPという形にするのが普通)。


NPやVPやSは、型とは呼ばずに、(統語)範疇と呼ぶのが普通ぽい。時制や人称のような情報は、統語素性として扱って、単数形名詞をN[sg]、複数形名詞をN[pl]などと書いたりする(場合もある)。この場合、実用上は、N[sg]とN[pl]は、特に何の関係もない、別の素性という扱いになる。現実的に、文法的に許される文章を過不足なく記述しようと思うと、かなり多くの統語素性は必要になると思う(時制や人称がないと、"He meet ..."とか"him looks ..."みたいな文章も合法となったりする)。



組み合わせ範疇文法(CCG)。範疇文法には関数適用しかないけど、ラムダ計算で意味付けされることを考えると、関数合成を許すのは自然に思える。関数適用が左と右の2つあるので、合成のバリエーションも多い。関数合成があると、

He conjectured and might prove completeness.

のような文章で、

conjectured |- (S\NP)

might |- (S\NP)/(S\NP)

prove |- (S\NP)

から、

might prove |- S\NP

とできて、conjecturedとmight provedが同じ型になって対をなすようにできる。andは、forall x,(x\x)/xという多相関数みたいな型として実現されるとする。他に、

(X/Y)/Z Y/Z => X/Z

(X/Y)\Z (Y\Z) => X\Z

Y\Z (X\Y)\Z => X\Z

Y/Z (X\Y)/Z => X/Z

のような代入規則も、CCGでは使われる


ここまでは、いわば"二項演算"であるけども、CCGには、もう一つ、type raisingというのもあって、これは、一項演算みたいなもの。インフォーマルには、S->NP VPをCG/CCGに変換するとき、VP=S\NPとしても、NP=S/VPとしてもよいという任意性があり(意味論の自然さから、普通は、NPを基本範疇にして、VP=S\NPとする)、VP=S\NPと、NP=S/VPを組み合わせると、NP=S/(S\NP)と"解ける"。一般に、XとY=T\Xがあれば、X=T/Yでもあり、X=T/(T\X)と解ける。これがtype raisingの由来と思う

#左右の区別を忘れると、Type raisingは、イメージ的には、型Aから、多相型forall X,(A->X)->Xを作り出すような操作と見なせる。一方、多相型ラムダ計算では、最小不動点のChurch Encodingとして、forall X,(F(X) -> X) -> Xというのがあって、その特殊な場合と思える。

#日本語(トルコ語とかもそうらしい)などは、語順の自由度が英語より高い(と一般的に言われる)が、それらを系統的に扱うために、scramblingというものも追加する場合があるらしい


type raisingの恩恵として、関係代名詞の扱いが自然になるという説明が、よく書いてある。例えば、the company which Microsoft bought

だと、(boughtは他動詞なので)

Microsoft |- NP

bought |- (S\NP)/NP

だけど、type raisingによって

Microsoft |- S/(S\NP)

となって、関数合成によって

Microsoft bought |- S/NP

となって、

which |- (N\N)/(S\NP)

が適用できる形になる



意味が扱えるという話。例として、"0"と"1"と"+"だけで構成される算術式を、範疇文法(<=>文脈自由文法)で考える(個人的な感覚としては、CCGは、数式やプログラミング言語のような形式言語を扱うには、自由度が高すぎて、制御しづらいという気がする)。この場合、"0","1","+"の範疇として

"0" |- S : 0

"0" |- S\S : (\x -> 10*x)

"1" |- S : 1

"1" |- S\S : (\x -> 10*x + 1)

"+" |- (S\S)/S : (\x y -> x+y)

があるとする("001"のような表記も許容する)。':'の後ろにあるのが、意味。範疇文法なので、右適用と左適用だけ考えればよく、右適用を>で、左適用を<で表すと、"110+1"の導出木は、例えば

(< (< (< "1"|S "1"|(S\S)) "0"|(S\S)) (> "+"|( (S\S)/S ) "1"|S))

みたいになる。対応する意味を

(< (< (< 1 (\x -> 10*x + 1)) (\x -> 10*x)) (> (\x y -> x+y) 1))

のように置き換えて、簡約化していくと

(< (< (< 1 (\x -> 10*x + 1)) (\x -> 10*x)) (> (\x y -> x+y) 1))

=> (< (< 11 (\x -> 10*x)) (> (\x y -> x+y) 1))

=> (< 110 (> (\x y -> x+y) 1))

=> (< 110 (\y -> 1+y))

=> 111

という風に計算できる。というわけで、"110+1"の意味は、今の場合、111となる。


どういう意味を付けるかは、CGやCCGとは独立の話で、JSONパーサーなんかだと、対応するJSONオブジェクトが意味ということになるし、英語の意味を表すのに、日本語を使ってもいい。自然言語形式意味論としては、discourse representation theory(DRT)というものが比較的よく使われているっぽい(あまり調べていない)。一階述語論理との違いがよく分からないけど、東ロボくんの数学部分は、問題を計算機代数システムで処理できる形に変換するのに、CCG+DRTを使っているらしい(と書いてある)

cf.深い言語理解と数式処理の接合による入試数学問題解答システム

https://kaigi.org/jsai/webprogram/2013/pdf/622.pdf



CCGはCFGより計算能力が高いことが知られている(文脈依存文法よりは弱いらしい)けど、CCGまでいくと、主に関数合成を許しているせいで、導出木は、一般に非常に沢山ある。が、意味の方で見ると同じになる擬似的曖昧性(spurious ambiguity。導出木は違うけど、意味は同じと見なせる)と、そもそも、異なる意味に対応する、本質的な曖昧性が存在する

(本質的曖昧性の例1)有名な例

I saw a girl with a telescope

で、withが動詞句"saw a girl"にかかるか、名詞句"a girl"にかかるか、で、withの範疇は、( (S\NP)\(S\NP) )/NPと(NP\NP)/NPになって、意味的に異なる

(本質的曖昧性の例2)

The fish twisted and turned on the bent hook.

で、"on"の範疇が、( (S\NP)\(S\NP) )/NPであるとすると、"on the bent hook"は、turnedのみにかかると見ることもできるし、"twisted and turned"にかかると見ることもできる(あくまで、文法的には)。というわけで、全ての範疇が決まっても、、意味の異なる導出木が複数存在する場合がある



慣用表現の扱い。特に体系的な扱いは存在しない。英語には、色々な慣用表現がある。特に、句動詞(phrasal verb)というものが多い。look at ...みたいなの。こういう時には、at ...がlookにかかるというより、look atが他動詞であるかのように見たい。こういう時は、"at ..."の範疇をPP(前置詞句)とする。atの範疇はPP/NPとなる。lookの範疇は(S\NP)/PPとなる。look atの範疇は合成によって(S\NP)/NPとなるので、目的は果たせられる。ただ、どこからが句動詞とするのか、明確な線引きはない気がする。take ... apartみたいなパターンを慣用表現として捉えたければ、apartの範疇を何か新しいXPとして、takeの範疇を(VP/XP)/NPとすれば、一応、そういう扱いもできる


別の慣用表現の例。"It is easy to show that."とか"There is ...."みたいな文で、ItやThereに、もはや意味はない。これらは、CCGBankなどでは、それぞれNP[expl]とか、NP[thr]という範疇を持ち、isは((S\NP[expl])/VP[to])/VP[adj]とか(S\NP[thr])/NPみたいな特殊な範疇となる。


実際の構文解析は、例えばCYK法で行うことができる。英語のCCG parserを手っ取り早く試すには、

EasyCCG

http://homepages.inf.ed.ac.uk/s1049478/easyccg.html

とかがいい気がする(多分、これはProbabilistic CCG parserで、CCGBankで学習してるのだと思う。命令文は、実質対応できてないっぽい?)(CCGBankはPenn Treebankから機械的に生成されたものだと思う)(何故か、自然言語parser業界では、Javaが使われることが多く、これもJavaで書かれている)。python nltkにもCCG parserはあるけど、lexiconを自分で定義する必要がある。


日本語のCCG parserは、自分で作るしかない気がする

Maxwell方程式の表現論(2)Gupta-Bleuler量子化の表現論的側面

Photons and gravitons in conformal field theory

http://link.springer.com/chapter/10.1007%2F3540171630_73

という論文の2ページ目の冒頭を見たら、次のようなことが書いてあった。



$V$を、Poincare代数のhelicity 0のmassless既約ユニタリ表現$V_0$と、Poincare代数の自然な4次元表現(並進成分は自明に働くので、実質的に$so(3,1)$のベクトル表現)のテンソル積表現空間とする。$V$は、ユニタリー表現ではないし、既約でもないが、直既約表現になっているらしい。$V$の重要な部分表現として、

$V_L$:Lorenz gauge条件を満たす$V$の部分空間

$V_g$:A_{¥mu}=¥partial_{¥mu} ¥Lambdaと書けるような$V$の元全体

があって、$V_g$は$V_L$の部分空間でもある。そして、以下の3つの命題が成立する

(1)商表現$V/V_L$は、Poincare代数のhelicity 0のmassless既約ユニタリ表現と同値

(2)商表現$V_L/V_g$は、Poincare代数のhelicity +1/-1の2つのmassless既約ユニタリ表現の直和表現と同値

(3)$V_g$は、Poincare代数のhelicity 0のmassless既約ユニタリ表現と同値


Gupta-Bleuler量子化の時は、非負ノルムを持つ状態の部分空間を0ノルム空間で割って、正ノルムを持つ状態を取り出していた。同じことを、1粒子状態に対してやっているのが、上の操作と考えられ、$(V,V_L,V_g)$を指して、Gupta-Bleuler tripletと呼んでいる文献もある。$V$には、平面波のような、ノルムが定まらない状態は含まれておらず、高々可算個の基底を持つ、"小さな"空間になっている


#重力波は、Klein-Gordon方程式を満たす、10成分の場である。Poincare代数の自然表現の2次の対称テンソル積表現は10次元表現で、これとhelicity 0のmassless既約ユニタリ表現のテンソル積によって、新しい表現を作ることができる。この表現から出発して、上記と類似の方法で、helicity +2/-2のmassless既約ユニタリ表現を得ることができるっぽい(要確認)


Poincare代数のhelicity +1/-1の2つのmassless既約ユニタリ表現の直和表現は、Maxwell方程式の解空間と同型になっているというのが、twistor理論の教えだった(多分、Maxwell方程式のノルム有限な解は、これ(の完備化したもの)で尽きると思うのだけど、証明とかは知らない)。この解空間は、$V_L/V_g$と同型であるが、一方、$V_L$の元は、電磁ポテンシャルなので、標準的な電磁気学の知識から、対応する電磁場を直接作ることができ、これはMaxwell方程式を満たす。真空中の電磁場は、波動方程式を満たすので、helicity 0のmassless既約ユニタリ表現と、Poincare代数の6次元表現(自然表現の2次の反対称テンソル積表現)のテンソル積空間$V_2$に含まれる解を考えるのは自然。そうすると、Poincare代数の作用と可換な写像

d_1 : V_L ¥to V_2

があって、その像$Im(d_1)$は、$V_2$の中で、Maxwell方程式を満たす元全体に等しいと期待される。そして

Ker(d_1) ¥simeq V_g

であれば、$Im(d_1)$は、(2)の表現空間$V_L/V_g$と同型になるので、辻褄が合う。


$V_0$は、共形代数$su(2,2)$の最低ウェイト表現に一意に拡張できるので、具体的な基底を取って、Poincare代数の作用を明示的に書くことができる。同様に、$V_L$や$V_g$の基底も具体的に記述することができるはず。ここでは、$V_0$を、以下のように、4変数多項式環の部分空間として実現する

V_0 = ¥{f ¥in ¥mathbf{C}¥[z_1,z_2,z_3,z_4¥] | (z_1 ¥partial_1 + z_2 ¥partial_2 - z_3 ¥partial_3 - z_4 ¥partial_4)f = 0¥}

$su(2,2)$の複素化sl(4,¥mathbf{C})のChevalley生成元の作用は、

e_1 = z_1 ¥partial_2

f_1 = z_2 ¥partial_1

h_1 = z_1 ¥partial_1 - z_2 ¥partial_2

e_2 = - z_2 z_3

f_2 = ¥partial_2 ¥partial_3

h_2 = 1 + z_2 ¥partial_2 + z_3 ¥partial_3

e_3 = - z_4 ¥partial_3

f_3 = - z_3 ¥partial_4

h_3 = - z_3 ¥partial_3 + z_4 ¥partial_4

のようになる。1が最低ウェイトベクトルとなり、これに上記の元を作用させていくことで、$V_0$の基底が得られる。別に、この実現を使うことに特に意味はないけど、コンピュータに計算させるときに、ちょっと楽。変数の数はたまたま4つだけど、これ自体は、通常の時空の座標とは関係ない


#$V_0$に含まれる同次多項式の次数は、全て偶数で、次数2n(n=0,1,2,...)の同次多項式が張る部分空間の次元は、$(n+1)^2$になる。これは、$V_0$の定義によって、z1,z2の部分の次数と、z3,z4の部分の次数が等しくないといけないことから、容易に分かる。

e_4 = ¥[e_1,e_2¥] = -z_1 z_3

e_5 = ¥[e_2,e_3¥] = -z_2 z_4

e_6 = ¥[e_1,e_5¥] = -z_1 z_4

P_0 = -(e_2+e_6) = z_2 z_3 + z_1 z_4

P_1 = -¥sqrt{-1}(e_4 - e_5) = ¥sqrt{-1}(z_1z_3 - z_2 z_4)

P_2 = -(e_4+e_5) = z_1 z_3 + z_2 z_4

P_3 = e_6 - e_2 = z_2 z_3 - z_1 z_4

となる。P_0^2 - P_1^2 - P_2^2 - P_3^2は恒等的に0なので、質量0のKlein-Gordon方程式が満たされる。空間回転$M_1,M_2,M_3$及びLorentzブースト$L_1,L_2,L_3$は、

M_1 = -(z_1 ¥partial_2 - z_2 ¥partial_1 - z_3 ¥partial_4 + z_4 ¥partial_3)/2

M_2 = ¥sqrt{-1}(z_1 ¥partial_2 - z_3¥partial_4 - z_4 ¥partial_3 + z_2 ¥partial_1)/2

M_3 = -¥sqrt{-1}(z_1 ¥partial_1 - z_2 ¥partial_2 + z_3 ¥partial_3 - z_4 ¥partial_4)/2

L_1 = -¥sqrt{-1}(z_1 ¥partial_2 + z_3 ¥partial_4 - z_2 ¥partial_1 - z_4 ¥partial_3)/2

L_2 = (-z_4¥partial_3 - z_1 ¥partial_2 - z_2 ¥partial_1 - z_3 ¥partial_4)/2

L_3 = (z_1 ¥partial_1 - z_2 ¥partial_2 - z_3 ¥partial_3 + z_4 ¥partial_4)/2



以降は、普通の電磁気学でやる計算と同じになる。

V = V_0 ¥otimes ¥mathbf{C}^4として、

$w_0 = (1,0,0,0)$

$w_1 = (0,1,0,0)$

$w_2 = (0,0,1,0)$

$w_3 = (0,0,0,1)$

を$\mathbf{C}^4$の基底とする

$d_0 : V_0 \to V$を

d_0(f) = (e_2+e_6)(f) ¥otimes w_0 - ¥sqrt{-1} (e_4-e_5)(f) ¥otimes w_1 - (e_4+e_5)(f) ¥otimes w_2 -(e_2-e_6)(f) ¥otimes w_3 = (-P_0(f) , P_1(f) , P_2(f) , P_3(f))

で定義し、

$V_g = Im(d_0)$とする。$d_0$は、Poincare代数の作用と可換なはずで、従って、$V_g$は$V_0$と同値な表現空間となる。$V_g$の最低ウェイトベクトルは

(z_2 z_3+z_1z_4 , ¥sqrt{-1}(z_1z_3 - z_2 z_4) , z_1 z_3 + z_2 z_4 , z_2 z_3 - z_1 z_4)

になる。第一成分がスカラーポテンシャル、残りがベクトルポテンシャルに相当する


$d_0$とPoincare代数の作用の可換性は、例えば

(1 ¥otimes M_1 + M_1 ¥otimes 1) d_0 = (-M_1 P_0 , M_1 P_1 , M_1 P_2-P_3 , M_1 P_3+P_2) = (-P_0 M_1, P_1 M_1 , P_2 M_1 , P_3 M_1)

などで、この場合、

¥[M_1,P_2¥] = P_3

¥[M_1,P_3¥] = -P_2

より従う。他の生成元についても、同様に計算すればいい


次に、d_1 : V_0 ¥otimes ¥mathbf{C}^4 ¥to V_2は、通常の外微分と同様、A = (A_0,A_1,A_2,A_3) ¥in V_0 ¥otimes ¥mathbf{C}^4に対して

d_1(A) = d_0(A_0) ¥wedge w_0 + d_0(A_1) ¥wedge w_1 + d_0(A_2) ¥wedge w_2 + d_0(A_3) ¥wedge w_3

で定義する。6成分出てくるけど、当然、電場と磁場に対応する


Mikowski計量に対する、Hodgeのスター作用素* : ¥wedge^2 ¥mathbf{C}^4 ¥to ¥wedge^2 ¥mathbf{C}^4

*(w_0 ¥wedge w_1) = -w_2 ¥wedge w_3

*(w_0 ¥wedge w_2) = -w_3 ¥wedge w_1

*(w_0 ¥wedge w_3) = -w_1 ¥wedge w_2

*(w_1 ¥wedge w_2) = w_0 ¥wedge w_3

*(w_3 ¥wedge w_1) = w_0 ¥wedge w_2

*(w_2 ¥wedge w_3) = w_0 ¥wedge w_1

で定義すると、$*^2=-1$なので、¥wedge^2 ¥mathbf{C}^4を$*$の固有値+i/-iの空間に固有分解できる。$*$は、Poincare代数の作用と可換なので、$V_2$も2つの表現空間の直和に分解する。この分解は、$Im(d_1)$に制限すると、Maxwell方程式の解空間が、正helicity/負helicityの解に分解する(それぞれ、自己双対・反自己双対Maxwell方程式の解空間)のと対応している。


自己双対・反自己双対Maxwell方程式の解空間は、共形代数の$su(2,2)$の最低ウェイト表現空間となっているので、最低ウェイトベクトルを計算するのが、構造を知る上で、手っ取り早い。$V_2$には、$su(2,2)$自体は作用しないが、Poincare代数の作用のみでも、最低ウェイトベクトル条件のいくつかは書くことができ、かつ、最低ウェイトベクトルを決定するにはそれで十分である。このへんは、表現論の知識が必要となる。最低ウェイト条件の詳細は

Maxwell方程式の表現論

http://d.hatena.ne.jp/m-a-o/20150621#p1

で書いたものと同じ。


具体的に、自己双対Maxwell方程式の"最低ウェイト解"を計算してみる。まず、$*$について、固有値+iの空間の基底として、

w_0 ¥wedge w_1 + ¥sqrt{-1} w_2 ¥wedge w_3

w_0 ¥wedge w_2 + ¥sqrt{-1} w_3 ¥wedge w_1

w_0 ¥wedge w_3 + ¥sqrt{-1} w_1 ¥wedge w_2

を取ることができる。この基底に対する、成分を(¥psi_1,¥psi_2,¥psi_3) ¥in V_0 ¥otimes ¥mathbf{C}^3とする。$V_0$の定義より

(z_1 ¥partial_1 + z_2 ¥partial_2 - z_3 ¥partial_3 - z_4 ¥partial_4)¥psi_i = 0

が分かる。$su(2,2)$作用のChevalley生成元のいくつかは、Poincare代数の元のみで書けて、例えば

h_1 = L_3 - ¥sqrt{-1} M_3 = z_1 ¥partial_1 - z_2 ¥partial_2

の場合、

h_1(¥psi_i) + 2 ¥psi_i = 0

が最低ウェイトベクトルの満たす条件となる。また、

f_1 = (M_1 - ¥sqrt{-1}M_2)/2 - ¥sqrt{-1}(L_1 - ¥sqrt{-1}L_2)/2 = z_2 ¥partial_1

に対して、

f_1(¥psi_i) = 0

なども条件としてある。ここまでの条件で、¥psi_1,¥psi_2,¥psi_3

z_2^2 z_3 z_4 , z_2^2z_3^2, z_2^2z_4^2

の3つの元の一次結合で書けることが分かる。他に、$h_2,h_3,f_3$などに関する条件も使うと、定数倍を除いて

¥psi_1 = -¥frac{¥sqrt{-1}}{2}(z_2^2 z_4^2 + z_2^2 z_3^2)

¥psi_2 = ¥frac{1}{2}(z_2^2 z_4^2 - z_2^2 z_3^2)

¥psi_3 = z_2^2 z_3 z_4

という解を得る。これらの解を与える電磁ポテンシャルとして、例えば

(¥phi,A_x,A_y,A_z) = (-z_2 z_4/2 , ¥sqrt{-1}z_2 z_3/2 , z_2 z_3/2 , -z_2 z_4/2)

などが取れる。反自己双対Maxwell方程式の最低ウェイト解も同様に計算できる。


光子の状態空間は、光子の1粒子状態の対称テンソル積空間として得られる。物理の教科書では、波数ごとに、(非可算無限個の)生成消滅演算子があって、それを真空に作用させて...ということになるが、平面波は、Maxwell方程式のノルム有限な解ではなく、$P_0,P_1,P_2,P_3$の連続スペクトルに対する一般化固有状態として得られるはず。そして、($V_0$の各基底に対応して)高々可算無限個の生成消滅演算子があって、それを真空に作用させていくことで、多光子状態が得られるという形になる。通常、物理の教科書に書いてある、量子電気力学の生成消滅演算子による光子の状態空間と、表現論的("群論的")な光子の一粒子状態空間の説明が、等価だというのは、長い間納得できなかったけど、数学的に正当化できそうではある



電磁気学をやるには、電子の状態空間も知りたい。Poincare代数のmassless既約表現については、共形代数の表現に拡張することで、その構造が数学的によく分かる(通常の群の誘導表現による構成では、詳しい構造を調べるのは結構難しい)のだけど、massive表現については、その手は使えない。が、de Sitter代数$so(4,1)$の主系列表現は、$so(4,1)$からPoincare代数へのcontraction(半径→∞の極限を取る)を行うことによって、Poincare代数のmassive既約ユニタリ表現2つの直和に分解するらしい。

GENERALIZED EIGENVECTORS AND GROUP REPRESENTATIONS -- THE CONNECTION BETWEEN REPRESENTATIONS OF SO(4, 1) AND THE POINCARE GROUPt

http://link.springer.com/chapter/10.1007/978-94-010-2669-7_8

$so(4,1)$の主系列表現は、最高ウェイトも最低ウェイトもないので、Lie環の表現論としてはまだ難しい部類ではある(というのは、現時点では、系統的な手法が使えず、個別に対処する必要があるという意味で)けども、Gelfand-Zetlin basisなど、具体的な基底を構成することができる(はず←やったことない)。これによって、Poincare代数のmassive既約表現の具体的な記述を手にすることができるはず(続く)


#contractionによる構成自体は、Euclid群の既約ユニタリ表現を作る場合にも、同じような方法が取られるので、さほど問題なく出来るはずではある。masslessの方の基底は分かっても、massiveの方がどうにもならなくて、ふと、Eudlic群なら、誰かやってる人いそうだ、と思って調べたら、q-Euclid代数でやってる人がいて、Poincare代数でも、同じことを考えた人が1960年代にいたという...


#http://www.kyoritsu-pub.co.jp/90thmath/math03.html

によると『リー群のユニタリ表現論』という本が刊行予定らしく、"非コンパクトリー群の例として一般化ローレンツ群SO(n,1)と擬ユニタリ群SU(n,1)をとり,それらの(非ユニタリ)主系列表現の構造を,上のGZ基底を拡張した基底を作って詳細に調べることにより"とあり、欲しい情報を書いてくれてそうな感じであるが、いつ出版されるか不明なので、自分で計算する必要がありそう。ぐぬぬ


以下は、計算の一部をRisa/Asirで確かめた時のコード。文章と式が合ってなかったら、文章の方が間違ってると思う

/* Risa/Asir code */

def assert(V){
   if(!V){error(V);}
   return 1;
}


V=[z1,z2,z3,z4,d1,d2,d3,d4];
Z1=dp_ptod(z1,V);
Z2=dp_ptod(z2,V);
Z3=dp_ptod(z3,V);
Z4=dp_ptod(z4,V);
D1=dp_ptod(d1,V);
D2=dp_ptod(d2,V);
D3=dp_ptod(d3,V);
D4=dp_ptod(d4,V);

P0 = dp_ptod(z2*z3+z1*z4 , V);
P1 = dp_ptod(@i*(z1*z3-z2*z4) , V);
P2 = dp_ptod(z1*z3 + z2*z4 , V);
P3 = dp_ptod(z2*z3 - z1*z4 , V);
M1 = -(dp_weyl_mul(Z1,D2) - dp_weyl_mul(Z2,D1) - dp_weyl_mul(Z3,D4) + dp_weyl_mul(Z4,D3))/2;
M2 = @i*(dp_weyl_mul(Z1,D2) - dp_weyl_mul(Z3,D4) - dp_weyl_mul(Z4,D3) + dp_weyl_mul(Z2,D1))/2;
M3 = -@i*(dp_weyl_mul(Z1,D1) - dp_weyl_mul(Z2,D2) + dp_weyl_mul(Z3,D3) - dp_weyl_mul(Z4,D4))/2;

L1 = -@i*(dp_weyl_mul(Z1,D2) + dp_weyl_mul(Z3,D4) - dp_weyl_mul(Z2,D1) - dp_weyl_mul(Z4,D3))/2;
L2 = (-dp_weyl_mul(Z4,D3) - dp_weyl_mul(Z1,D2) - dp_weyl_mul(Z2,D1) - dp_weyl_mul(Z3,D4))/2;
L3 = (dp_weyl_mul(Z1,D1) - dp_weyl_mul(Z2,D2) - dp_weyl_mul(Z3,D3) + dp_weyl_mul(Z4,D4))/2;


RM1 = newmat(4,4,[[0,0,0,0],[0,0,0,0],[0,0,0,-1],[0,0,1,0]]);
RM2 = newmat(4,4,[[0,0,0,0],[0,0,0,1],[0,0,0,0],[0,-1,0,0]]);
RM3 = newmat(4,4,[[0,0,0,0],[0,0,-1,0],[0,1,0,0],[0,0,0,0]]);
RL1 = newmat(4,4,[[0,1,0,0],[1,0,0,0],[0,0,0,0],[0,0,0,0]]);
RL2 = newmat(4,4,[[0,0,1,0],[0,0,0,0],[1,0,0,0],[0,0,0,0]]);
RL3 = newmat(4,4,[[0,0,0,1],[0,0,0,0],[0,0,0,0],[1,0,0,0]]);

assert( dp_weyl_mul(P0,M1)==dp_weyl_mul(M1,P0) );
assert( dp_weyl_mul(P1,M1)==dp_weyl_mul(M1,P1) );
assert( dp_weyl_mul(L1,P2)==dp_weyl_mul(P2,L1) );
assert( dp_weyl_mul(L1,P3)==dp_weyl_mul(P3,L1) );
assert( dp_weyl_mul(P0,L1)-dp_weyl_mul(L1,P0)==P1 );
assert( dp_weyl_mul(P1,L1)-dp_weyl_mul(L1,P1)==P0 );
assert( dp_weyl_mul(M1,P2)-dp_weyl_mul(P2,M1)==P3 );
assert( dp_weyl_mul(M1,P3)-dp_weyl_mul(P3,M1)==-P2 );

assert( dp_weyl_mul(M1,M2)-dp_weyl_mul(M2,M1)==M3 );
assert( dp_weyl_mul(M2,M3)-dp_weyl_mul(M3,M2)==M1 );
assert( dp_weyl_mul(M3,M1)-dp_weyl_mul(M1,M3)==M2 );
assert( dp_weyl_mul(L1,L2)-dp_weyl_mul(L2,L1)==-M3 );
assert( dp_weyl_mul(L2,L3)-dp_weyl_mul(L3,L2)==-M1 );
assert( dp_weyl_mul(L3,L1)-dp_weyl_mul(L1,L3)==-M2 );
assert( dp_weyl_mul(M1,L2)-dp_weyl_mul(L2,M1)==L3 );
assert( dp_weyl_mul(L1,M2)-dp_weyl_mul(M2,L1)==L3 );
assert( dp_weyl_mul(M2,L3)-dp_weyl_mul(L3,M2)==L1 );
assert( dp_weyl_mul(L2,M3)-dp_weyl_mul(M3,L2)==L1 );
assert( dp_weyl_mul(L3,M1)-dp_weyl_mul(M1,L3)==L2 );
assert( dp_weyl_mul(M3,L1)-dp_weyl_mul(L1,M3)==L2 );
assert( RM1*RM2-RM2*RM1==RM3 );
assert( RM2*RM3-RM3*RM2==RM1 );
assert( RM3*RM1-RM1*RM3==RM2 );
assert( RL1*RL2-RL2*RL1==-RM3 );
assert( RL2*RL3-RL3*RL2==-RM1 );
assert( RL3*RL1-RL1*RL3==-RM2 );
assert( RM1*RL2-RL2*RM1==RL3 );
assert( RL1*RM2-RM2*RL1==RL3 );
assert( RM2*RL3-RL3*RM2==RL1 );
assert( RL2*RM3-RM3*RL2==RL1 );
assert( RM3*RL1-RL1*RM3==RL2 );
assert( RL3*RM1-RM1*RL3==RL2 );

/* ---------------------------------- */
P0 = z2*z3+z1*z4;
P1 = @i*(z1*z3-z2*z4);
P2 = z1*z3 + z2*z4;
P3 = z2*z3 - z1*z4;

/* lowest weight solution of SD-Maxwell equation */
Ex = -@i*(z2^2*z4^2 + z2^2*z3^2)/2;
Ey = (z2^2*z4^2 - z2^2*z3^2)/2;
Ez = z2^2*z3*z4;
Bx = -@i*Ex;
By = -@i*Ey;
Bz = -@i*Ez;

/* half-part of Maxwell equation */
assert( P1*Ex + P2*Ey + P3*Ez==0 );  /* div E = 0 */
assert( P2*Ez-P3*Ey==-P0*Bx );
assert( P3*Ex-P1*Ez==-P0*By );
assert( P1*Ey-P2*Ex==-P0*Bz );

/* lowest weight potential of SD-Maxwell equation */
A0 = -z2*z4/2;
A1 = @i*z2*z3/2;
A2 = z2*z3/2;
A3 = -z2*z4/2;

assert( -P1*A0 - P0*A1 == Ex );
assert( -P2*A0 - P0*A2 == Ey );
assert( -P3*A0 - P0*A3 == Ez );
assert( P2*A3-P3*A2 == Bx );
assert( P3*A1-P1*A3 == By );
assert( P1*A2-P2*A1 == Bz );

2015-06-21

Maxwell方程式の表現論

Maxwell方程式の表現論(0)

http://d.hatena.ne.jp/m-a-o/20150125#p2

からの続き(本編)。


Poincare代数の離散ヘリシティを持つmassless既約ユニタリ表現は、共形代数$su(2,2)$のladder representationに一意に拡張できて、twistor理論では、Dolbeaultコホモロジーを使って、この表現空間の実現が与えられた。一方で、共形変換はMaxwell方程式の解空間にも自然に作用していて、DolbeaultコホモロジーからMaxwell方程式の(適当なクラスの)解空間(正確には、自己双対Maxwell方程式と反自己双対Maxwell方程式それぞれの解空間)への表現の同値を与えるintertwinning operatorは、Penrose変換によって具体的に実現される。と言いたいのだけど、Penrose変換の像がどうなっているか、明示的に書いてある文献を見つけることができなかった。


というわけで、真空中の自己双対Maxwell方程式(SD-Maxwell方程式)と反自己双対Maxwell方程式(ASD-Maxwell方程式)の解の中で、ladder representationと同値な表現を与えるような基底を見つけることを考える。twistor空間上の"関数"は、物理的な解釈がしづらいものなので、Minkowski空間上の場に対して、ladder representationが実現できてしまえば、Penrose変換とかコホモロジーとか難しいものは忘れてしまってもよい



方針。$su(2,2)$のladder representationは、lowest weight representationであって、この既約ユニタリー表現を、純表現論的に記述することは(原理的には)難しくない。Verma加群の商としての直接的な記述は、例えば

THE MASSLESS REPRESENTATIONS OF THE CONFORMAL QUANTUM ALGEBRA

http://iopscience.iop.org/0305-4470/27/14/012

などで見ることができる。表現の構造は、これで分かったとして、SD/ASD-Maxwell方程式の解の中で、lowest weight vectorに対応するものを発見して、それに対して、無限小共形変換を、どんどん作用させていけば、ASD/SD-Maxwell方程式の解空間の基底が得られる(はず)。これがladder representationと同値であることを確認するには、lowest weight vectorが満たす最低ウェイト条件の確認と、singular vectorがちゃんと消えることを言えばいい。


lowest weight vectorは適当なPenrose変換の像の中にあると信じて探すのだが、

On the density of twistor elementary states.

http://projecteuclid.org/euclid.pjm/1102637077

に、elementary statesとは$K$-finite vector($K$は$SU(2,2)$の極大コンパクト部分群)と同じものであると書いてある。なので、elementary statesの像の中で探せばよい。勿論、最低ウェイト条件もsingular vectorが消えるという条件も線形微分方程式になるので、直接解いてもいいけど


#これらの基底は、いずれも"ハミルトニアン"(ここでは、Hamiltonianは、Poincare代数の時間並進演算子のこととする)の固有状態ではない(多分、エネルギー固有状態は、一つも存在しない。これは、量子力学に於いて、平面波は、ラプラシアンの"固有状態"であるが、二乗可積分でない、というのと似た事態だと思う)。ところで、Weinbergの教科書「場の量子論」1巻の2.5節の冒頭を見ると、"4元エネルギー・運動量テンソルは互いに可換だから、物理的状態ベクトルをこの4元ベクトルの固有値を使って表すことは自然だ."と書いてある。これは、連続スペクトルと固有値を同じように扱おうということだと思うけど、Weinbergの教科書に限らず、物理の教科書は、大体こういう議論で始まる気がする(Wignerの論文自体で、こういう議論がされていたっぽい。cf: arXiv:0809.4942)。


準備。まず、spinor場への共形変換の作用を具体的に知る必要がある。これは

Finite-component field representations of the conformal group

http://www.sciencedirect.com/science/article/pii/0003491669902784

で調べられている。論文では、spinor場とは限らない、一般の多成分場を扱っていて、共形変換の多成分場への表現は、適当な1-cocycleでtwistした誘導表現で実現されると考え、このcocycleが並進に対しては自明である時、スケール変換・Lorentz変換・特殊共形変換のなす群(の被覆群)の有限次元表現とcocycleが対応する。という具合であるが、今考えるのはspinor場であるので、特殊共形変換の作用は自明であるとして、結局、$Spin(3,1)$の有限次元表現から、cocycleを作り、誘導表現を得ることになる。結局、これは多くの場の理論の教科書に書いてあることと同じで、massless spinor場の場合、$Spin(3,1)$の有限次元表現として、$D(j,0)$あるいは$D(0,j)$の形のものを取る($j$は非負の整数or半整数で、それぞれhelicityが+jあるいは-jの場に対応する)。こうして、spinor場への共形変換の作用が分かる。$Spin(3,1)$の有限次元表現と$so(3,1)$の有限次元表現は、1:1に対応するので、以下、両者を混同する


#Poincare代数の作用は、Klein-Gordon作用素と可換(Klein-Gordon作用素はPoincare代数のCasimir元)であるけども、特殊共形変換やスケール変換は、Klein-Gordon作用素と可換ではない(けど、Klein-Gordon方程式の解を別の解を移す)ので、同じ"対称性"といっても、Poincare代数の場合と、共形代数の場合では、少し意味が違う(対称性という言葉は、濫用されすぎな気もする)


練習。"spin 1"(helicity +1/-1)の場合を考える前に、masslessなスカラー粒子の場合を眺める。この場合、spinor場の満たす方程式は、質量0のKlein-Gordon方程式。Penrose変換を計算する必要があるのだけど、幸いなことに

(Pre-)Hilbert spaces in twistor quantization

http://arxiv.org/abs/1210.0349

で必要な計算をしてくれているので、カンニングすることにする。結論を先に書くと、

¥phi = ¥frac{1}{x_0^2 - x_1^2 - x_2^2 - x_3^2}

が最低ウェイトベクトルを与える。


まず、最低ウェイト条件+特異ベクトルが消える条件をChevalley生成元を使って書く。共形代数は、複素化すると、sl(4)ないしso(6)と同型なLie代数となり、これはrankが3である。最低ウェイト条件は、まず

f_1 ¥phi = f_2 ¥phi = f_3 ¥phi =0

で、Cartan部分環の作用は、非負の半整数(or整数)の組$(j_1,j_2)$と共形次元$d = j_1+j_2+1$を使って

(h_1 + 2 j_1)¥phi = (h_2 -d - j_1 - j_2)¥phi = (h_3 + 2 j_2)¥phi= 0

と書ける(ことが分かっている)。スカラー粒子の場合は、j1=j2=0で、

h_1 ¥phi = (h_2 - 1)¥phi = h_3 ¥phi = 0

となる。また、2つの特異ベクトルが存在し、それは、$e_1$と$e_3$をそれぞれ一回ずつ作用させて得られる。これが0でないといけないので

e_1 ¥phi = e_3 ¥phi = 0

という2つの条件も確認する必要がある。


#[宿題]ユニタリー性を満たすためには、j1=0かj2=0でないといけないが、ユニタリー性を諦めれば、j1とj2が共に0でないことも可能である。そして、そのような表現は、$so(3,1)$の有限次元既約表現$D(j1,j2)$に付随するspinor場と関連すると予想するのは自然に思える。例えば、ベクトルポテンシャル(あるいは、一般のベクトル場)は、表現$D(1/2,1/2)$に付随するspinor場である(とされている)。この場合、Gupta-Bleuler形式を思い出して、最低ウェイトVerma表現から、物理でよくやるように、ゴーストを取り除いて商表現を作ると、2つのladder表現の直和になるだろうか(特異ベクトルを除くのも、ゴーストを除くのも、ユニタリーになるように商を取るだけなので、作業としては同じようなもんである)。このVerma表現を純表現論的に記述することは、多分そんなに難しくない。一方で、電磁ポテンシャルの中で、(j1,j2)=(1/2,1/2)の最低ウェイトベクトル条件を満たすものを見出すことも、多分できる。この"最低ウェイトベクトル場"に共形代数を作用させて得られる表現は、Verma表現と同値なのか、それとも、その何らかの商表現となっているのだろうか。


次に、Chevalley生成元と、共形代数の"幾何学的な生成元"の関係を知らないといけない。天下り的であるが

e_1 = ¥frac{1}{2}(M_1 + i M_2 + i L_1 + L_2)

e_2 = ¥frac{1}{2}(-P_{0} - P_{3})

e_3 = -¥frac{1}{2}(M_1 - i M_2 -i L_1 - L_2)

h_1 = L_3 - i M_3

h_2 = -(D + L_3)

h_3 = L_3 + i M_3

f_1 = ¥frac{1}{2}(-M_1 + i M_2 - i L_1 - L_2)

f_2 = ¥frac{1}{2}(K_{0} - K_{3})

f_3 = ¥frac{1}{2}(M_1 + i M_2 - i L_1 + L_2)

が答え。$D$,$P_a$,$K_a$,$M_i$,$L_i$はそれぞれ、スケール変換、並進、特殊共形変換、空間回転、Lorentz boostのgenerator。そして、これらの生成元のspinor場への作用は、既に分かっているので、最低ウェイト条件+特異ベクトルの消滅を確認するのは、ただの計算練習となる。


(注意)スケール変換と特殊共形変換の生成子には、cocycleを作る時に使ったスケール変換の一次元表現から来る"余分な"項δ(論文"Finite-component field representations of the conformal group"参照)が付いている。δ=0とすると、普通の特殊共形変換の無限小生成子が得られるけど、スカラー粒子の場合でも、δは0でない(0としてしまうと、lowest weight vectorは特殊共形変換で不変でなくなり、困ったことになってしまう。当初スカラーの場合はδ=0でいいと思いこんでいて、半日くらい悩んでしまった。。)


本題。スカラーの場合が分かれば、"spin 1"の場合も、殆ど同じようなものである。spin 1の場合は、正ヘリシティ・負ヘリシティのそれぞれの場は、3成分を持つので、計算は面倒になる。面倒になるけど、計算自体は高校生でもできる程度のものである。末尾に、Risa/Asirによる計算を書いたので、そちらを参考(Phi2とPsi2が、ヘリシティ+1/-1の最低ウェイト状態)。空間回転とLorentz boostの作用には、$so(3)$の3次元表現に対応する項が加わるが、この表現の形を見る(AsirコードのT1,T2,T3を参照)と、twistor理論で扱われている、標準的な基底は、あんまりよくないという気がする(基底の取り方は、任意であって、良い悪いは完全に、人間の主観であるが、Riemann?Silberstein vectorのようなものを使って書く方がよいのかもしれない)。


あと、Minkowski空間上の場を考えると、共形変換群の作用は、純幾何学的に書けるので、Perelomovの意味でのgeneralized coherent states(lowest weight vectorへ$SU(2,2)$を作用させて得られる軌道の元)の具体的表示を得るのが容易になる(純表現論的に計算するのは、しんどい)。最低ウェイト条件から、特殊共形変換は、最低ウェイトベクトルを不変に保つし、並進は、純幾何学的な並進。

#このgeneralized coherent statesは、普通の量子光学で出てくるコヒーレント状態とは、全く別物で無関係。通常のコヒーレント状態は、様々な光子数状態の重ね合わせになっているけども、SU(2,2)-generalized coherent statesの方は、光子数は一個の状態(と考えるべきなのだと思う)。マクロなMaxwell方程式と考えてしまうと、物質がないのに電磁場が存在するというのは、奇妙な気がするし、その場合エネルギーや運動量は場の関数であるので、勝手に規格化することも物理的に変である



得られた最低ウェイト状態は、式としては、twistor理論の説明をしている文献で、時々見かけるものではある。それでも、こんな形で、ladder representationが実現されているというのは、ちょっと意外だった



/**** common stuffs ****/
load("noro_pd.rr");
load("noro_matrix.rr");

def assert(V){
   if(!V){error(V);}
   return 1;
}


def appdiff_rat(P , F , V){
   N = length(V)/2;
   if(type(P)==9){  /* 分散多項式表現の時 */
      DP = dp_dtop(P, V);
   }else if(type(P)<=2){
      DP = P;
   }else{
      return;
   }
   Terms = p_terms(DP , V , 2);
   for(Ret = 0 , Terms = Terms; Terms!=[] ; Terms=cdr(Terms)){
      T = car(Terms);
      if(type(F)==9){
           Tmp = dp_dtop(F , V);
      }else{
           Tmp = F;
      }
      if(T==1){
         Ret += p_nf(DP , V , V ,0)*Tmp;
      }else{
         C = p_nf(sdiv(DP-p_nf(DP,[T],V,0) , T) , V , V , 0);  /* 単項式Tの係数 */
         for(I = 0 ; I < N ; I++){
            Deg = deg(T , V[N+I]);
            for(J = 0 ; J < Deg ; J++){
               Tmp = red(diff(Tmp , V[I]));
            }
         }
         for(I = 0 ; I < N ; I++){
            Deg = deg(T , V[I]);
            Tmp *= V[I]^Deg;
         }
         Ret += C*Tmp;
      }
   }
   return red(Ret);
}


def appdiff_complex(P , F , V){
    if(type(P)==9){
        DP = dp_dtop(P,V);
    }else{
        DP = P;
    }
    RP = real(DP);
    IP = imag(DP);
    if(type(F)<=2){
        RF = real(F);
        IF = imag(F);
    }else if(type(F)==3){ /*有理式*/
        DN_conj = conj(dn(F));
        if(DN_conj != dn(F)){
           TF = red(DN_conj*nm(F)/(DN_conj*dn(F)));
        }else{
           TF = F;
        }
        RF = red(real(nm(TF))/dn(TF));
        IF = red(imag(nm(TF))/dn(TF));
    }
    Ret = red( appdiff_rat(RP , RF , V) - appdiff_rat(IP,IF,V) );
    Ret += @i*red( appdiff_rat(IP,RF,V) + appdiff_rat(RP,IF,V) );
    return Ret;
}

/* 
微分作用Pを式Fに作用させる 
Vは変数のリスト(dp_dtopなどで使われるのと同様の形式)
*/
def appdiff(P,F,V){
    if(type(P)!=6){
        if(type(F)<=3){
           return appdiff_complex(P,F,V);
        }else if(type(F)==5){   /* vector */
           Ret = newvect(length(F));
           for(I=0;I<length(F);I++){
              Ret[I] = appdiff_complex(P, F[I] , V);
           }
           return Ret;
        }
    }else{ /*行列*/
        Sz = size(P);
        N = Sz[0];
        M = Sz[1];
        if(type(F)==5){
           Ret = newvect(N);
           for (I = 0 ; I < N ; I++){
              for(J = 0 ; J < M ; J++){
                Ret[I] += appdiff_complex(P[I][J] , F[J] , V);
              }
              Ret[I] = red(real(nm(Ret[I]))/dn(Ret[I])) + @i*red(imag(nm(Ret[I]))/dn(Ret[I]));
           }
           return Ret;
        }else if(type(F)==6){
            Sz2 = size(F);
            N2 = Sz2[0];
            M2 = Sz2[1];
            Ret = newmat(N,M2);
            for (I = 0 ; I < N ; I++){
                for(J = 0 ; J < M2 ; J++){
                    for(K = 0 ; K < M ; K++){
                       Ret[I][J] += appdiff_complex(P[I][K] , F[K][J] , V);
                    }
                }
            }
            return Ret;
        }
    }
}


def dp_weyl_matmul(A,B){
    Sz1 = size(A);
    Sz2 = size(B);
    L = Sz1[0];
    M = Sz1[1];
    N = Sz2[0];
    Ret = newmat(L,N);
    for(I = 0; I < L ; I++){
       for(J = 0 ; J < N ; J++){
          for(K = 0 ; K < M ; K++){
               Ret[I][J] += dp_weyl_mul(A[I][K] , B[K][J]);
          }
       }
    }
    return Ret;
}


V=[x0,x1,x2,x3,d0,d1,d2,d3];
X0=dp_ptod(x0,V);
X1=dp_ptod(x1,V);
X2=dp_ptod(x2,V);
X3=dp_ptod(x3,V);
D0=dp_ptod(d0,V);
D1=dp_ptod(d1,V);
D2=dp_ptod(d2,V);
D3=dp_ptod(d3,V);
C2 = D0*D0 - D1*D1 - D2*D2 - D3*D3;  /* Klein-Gordon operator */

/**** The end of common stuffs ****/


/**** scalar field stuffs ****/
Phi0 = 1/(x0^2-x1^2-x2^2-x3^2);    /* lowest weight state of helicity-0 massless particle */

/* infinitesimal conformal transformation for scalar field */
D = dp_weyl_mul(X0,D0) + dp_weyl_mul(X1,D1) + dp_weyl_mul(X2,D2) + dp_weyl_mul(X3,D3) + dp_ptod(1,V);
P0=D0;
P1=D1;
P2=D2;
P3=D3;
K0 = 2*dp_weyl_mul(X0,D) - dp_weyl_mul(X0*X0-X1*X1-X2*X2-X3*X3,D0);  /* special conformal transformation */
K1 = -2*dp_weyl_mul(X1,D) - dp_weyl_mul(X0*X0-X1*X1-X2*X2-X3*X3,D1); /* special conformal transformation */
K2 = -2*dp_weyl_mul(X2,D) - dp_weyl_mul(X0*X0-X1*X1-X2*X2-X3*X3,D2); /* special conformal transformation */
K3 = -2*dp_weyl_mul(X3,D) - dp_weyl_mul(X0*X0-X1*X1-X2*X2-X3*X3,D3); /* special conformal transformation */
M1 = dp_weyl_mul(X2,D3) - dp_weyl_mul(X3,D2);
M2 = dp_weyl_mul(X3,D1) - dp_weyl_mul(X1,D3);
M3 = dp_weyl_mul(X1,D2) - dp_weyl_mul(X2,D1);
L1 = dp_weyl_mul(X0,D1) + dp_weyl_mul(X1,D0);  /* Lorentz boost */
L2 = dp_weyl_mul(X0,D2) + dp_weyl_mul(X2,D0);  /* Lorentz boost */
L3 = dp_weyl_mul(X0,D3) + dp_weyl_mul(X3,D0);  /* Lorentz boost */


/* Chevalley generators for sl(4) and so(6) */
E1 = (M1 + @i*M2)/2 + @i*(L1 + @i*L2)/2;
F1 = -(M1 - @i*M2)/2 - @i*(L1 - @i*L2)/2;
H1 = L3 - @i*M3;
E2 = (-P0 - P3)/2;
F2 = (K0 - K3)/2;
H2 = -(D + L3);
E3 = -(M1 - @i*M2)/2 + @i*(L1 - @i*L2)/2;
F3 = (M1 + @i*M2)/2 - @i*(L1 + @i*L2)/2;
H3 = L3 + @i*M3;

/* and other Cartan-Weyl Basis */
E4 = dp_weyl_mul(E1,E2) - dp_weyl_mul(E2,E1);   /* @i*(P1 + @i*P2)/2 */
E5 = dp_weyl_mul(E2,E3) - dp_weyl_mul(E3,E2);   /* -@i*(P1 - @i*P2)/2 */
E6 = dp_weyl_mul(E1,E5) - dp_weyl_mul(E5,E1);   /* (-P0+P3)/2 */
F4 = dp_weyl_mul(F1,F2) - dp_weyl_mul(F2,F1);   /* @i*(K1 - @i*K2)/2 */
F5 = dp_weyl_mul(F2,F3) - dp_weyl_mul(F3,F2);   /* -@i*(K1+@i*K2)/2 */
F6 = dp_weyl_mul(F1,F5) - dp_weyl_mul(F5,F1);   /* (K0+K3)/2 */
H4 = H1+H2;
H5 = H2+H3;
H6 = H1+H2+H3;

assert(appdiff(C2 , Phi0 , V)==0); /* Klein-Gordon equation */
assert(appdiff(F1,Phi0,V)==0);     /* lowest weight condition */
assert(appdiff(F2,Phi0,V)==0);     /* lowest weight condition */
assert(appdiff(F3,Phi0,V)==0);     /* lowest weight condition */
assert(appdiff(H1,Phi0,V)==0);     /* lowest weight condition */
assert(appdiff(H2,Phi0,V)==Phi0);  /* lowest weight condition */
assert(appdiff(H3,Phi0,V)==0);      /* lowest weight condition */
assert(appdiff(E1,Phi0,V)==0);     /* first singular vector */
assert(appdiff(E3,Phi0,V)==0);     /* second singular vector */


/**** The end of scalar field stuffs ****/


/**** "spin 1/2 fields" ****/
/* Pauli matrix */
S0 = newmat(2,2,[[1,0],[0,1]]);
S1 = newmat(2,2,[[0,1],[1,0]]);
S2 = newmat(2,2,[[0,-@i],[@i,0]]);
S3 = newmat(2,2,[[1,0],[0,-1]]);

/**** helicity +1/2 ****/
Phi = newvect(2,[(x1-@i*x2)/(x0^2-x1^2-x2^2-x3^2)^2 , (x0-x3)/(x0^2-x1^2-x2^2-x3^2)^2]);

DS = dp_weyl_mul(X0,D0) + dp_weyl_mul(X1,D1) + dp_weyl_mul(X2,D2) + dp_weyl_mul(X3,D3);
D = (DS + dp_ptod(3/2,V))*S0;
P0 = D0*S0;
P1 = D1*S0;
P2 = D2*S0;
P3 = D3*S0;
M1 = S0*(dp_weyl_mul(X2,D3) - dp_weyl_mul(X3,D2)) + dp_ptod(@i/2,V)*S1;
M2 = S0*(dp_weyl_mul(X3,D1) - dp_weyl_mul(X1,D3)) + dp_ptod(@i/2,V)*S2;
M3 = S0*(dp_weyl_mul(X1,D2) - dp_weyl_mul(X2,D1)) + dp_ptod(@i/2,V)*S3;
L1 = S0*(dp_weyl_mul(X0,D1) + dp_weyl_mul(X1,D0)) + dp_ptod(-1/2,V)*S1;
L2 = S0*(dp_weyl_mul(X0,D2) + dp_weyl_mul(X2,D0)) + dp_ptod(-1/2,V)*S2;
L3 = S0*(dp_weyl_mul(X0,D3) + dp_weyl_mul(X3,D0)) + dp_ptod(-1/2,V)*S3;
K0 = (2*dp_weyl_mul(X0,DS) - dp_weyl_mul(X0*X0-X1*X1-X2*X2-X3*X3,D0))*S0 + 3*X0*S0 - X1*S1 - X2*S2 - X3*S3; 
K1 = (-2*dp_weyl_mul(X1,DS) - dp_weyl_mul(X0*X0-X1*X1-X2*X2-X3*X3,D1))*S0 - 3*X1*S0 + X0*S1 - @i*X2*S3 + @i*X3*S2;
K2 = (-2*dp_weyl_mul(X2,DS) - dp_weyl_mul(X0*X0-X1*X1-X2*X2-X3*X3,D2))*S0 - 3*X2*S0 + X0*S2 + @i*X1*S3 - @i*X3*S1;
K3 = (-2*dp_weyl_mul(X3,DS) - dp_weyl_mul(X0*X0-X1*X1-X2*X2-X3*X3,D3))*S0 - 3*X3*S0 - @i*X1*S2 + @i*X2*S1 + X0*S3;

E1 = (M1 + @i*M2)/2 + @i*(L1 + @i*L2)/2;
F1 = -(M1 - @i*M2)/2 - @i*(L1 - @i*L2)/2;
H1 = L3 - @i*M3;
E2 = (-P0 - P3)/2;
F2 = (K0 - K3)/2;
H2 = -(D + L3);
E3 = -(M1 - @i*M2)/2 + @i*(L1 - @i*L2)/2;
F3 = (M1 + @i*M2)/2 - @i*(L1 + @i*L2)/2;
H3 = L3 + @i*M3;
E4 = dp_weyl_matmul(E1,E2) - dp_weyl_matmul(E2,E1);
E5 = dp_weyl_matmul(E2,E3) - dp_weyl_matmul(E3,E2);
E6 = dp_weyl_matmul(E1,E5) - dp_weyl_matmul(E5,E1);

/* parts of sl(4) relation */
assert( dp_weyl_matmul(E1,F1) - dp_weyl_matmul(F1,E1)==H1 );
assert( dp_weyl_matmul(E2,F2) - dp_weyl_matmul(F2,E2)==H2 );
assert( dp_weyl_matmul(E3,F3) - dp_weyl_matmul(F3,E3)==H3 );
assert( dp_weyl_matmul(E2,F1)==dp_weyl_matmul(F1,E2) );
assert( dp_weyl_matmul(H1,E1) - dp_weyl_matmul(E1,H1)==2*E1 );
assert( dp_weyl_matmul(H2,E2) - dp_weyl_matmul(E2,H2)==2*E2 );


assert( appdiff(D0*S0+D1*S1+D2*S2+D3*S3,Phi,V)==0*Phi );   /* Weyl equation */
assert( appdiff(F1,Phi,V)==0*Phi );
assert( appdiff(F2,Phi,V)==0*Phi );
assert( appdiff(F3,Phi,V)==0*Phi );
assert( appdiff(H1,Phi,V)+Phi==0*Phi );          /* j1==1/2 */
assert( appdiff(H2,Phi,V)-2*Phi==0*Phi );        /* d+j1+j2==2 */
assert( appdiff(H3,Phi,V)==0*Phi );              /* j2==0 */
/* singular vector vanishing */
assert( appdiff(E1,appdiff(E1,Phi,V),V)==0*Phi );
assert( appdiff(E3,Phi,V)==0*Phi ); 
assert( appdiff(E4,Phi,V)-appdiff(E2,appdiff(E1,Phi,V),V)==0*Phi );
assert( appdiff(E6 , appdiff(E2,Phi,V) , V) - appdiff(E4 , appdiff(E5,Phi,V) , V)==0*Phi);


/**** helicity -1/2 ****/
Psi = newvect(2, [-x0+x3 , x1+@i*x2])/(x0^2-x1^2-x2^2-x3^2)^2;
DS = dp_weyl_mul(X0,D0) + dp_weyl_mul(X1,D1) + dp_weyl_mul(X2,D2) + dp_weyl_mul(X3,D3);
D = (DS + dp_ptod(3/2,V))*S0;
P0 = D0*S0;
P1 = D1*S0;
P2 = D2*S0;
P3 = D3*S0;
M1 = S0*(dp_weyl_mul(X2,D3) - dp_weyl_mul(X3,D2)) + dp_ptod(@i/2,V)*S1;
M2 = S0*(dp_weyl_mul(X3,D1) - dp_weyl_mul(X1,D3)) + dp_ptod(@i/2,V)*S2;
M3 = S0*(dp_weyl_mul(X1,D2) - dp_weyl_mul(X2,D1)) + dp_ptod(@i/2,V)*S3;
L1 = S0*(dp_weyl_mul(X0,D1) + dp_weyl_mul(X1,D0)) + dp_ptod(1/2,V)*S1;
L2 = S0*(dp_weyl_mul(X0,D2) + dp_weyl_mul(X2,D0)) + dp_ptod(1/2,V)*S2;
L3 = S0*(dp_weyl_mul(X0,D3) + dp_weyl_mul(X3,D0)) + dp_ptod(1/2,V)*S3;
K0 = (2*dp_weyl_mul(X0,DS) - dp_weyl_mul(X0*X0-X1*X1-X2*X2-X3*X3,D0))*S0 + 3*X0*S0 + X1*S1 + X2*S2 + X3*S3; 
K1 = (-2*dp_weyl_mul(X1,DS) - dp_weyl_mul(X0*X0-X1*X1-X2*X2-X3*X3,D1))*S0 - 3*X1*S0 - X0*S1 - @i*X2*S3 + @i*X3*S2;
K2 = (-2*dp_weyl_mul(X2,DS) - dp_weyl_mul(X0*X0-X1*X1-X2*X2-X3*X3,D2))*S0 - 3*X2*S0 - X0*S2 + @i*X1*S3 - @i*X3*S1;
K3 = (-2*dp_weyl_mul(X3,DS) - dp_weyl_mul(X0*X0-X1*X1-X2*X2-X3*X3,D3))*S0 - 3*X3*S0 - @i*X1*S2 + @i*X2*S1 - X0*S3;

E1 = (M1 + @i*M2)/2 + @i*(L1 + @i*L2)/2;
F1 = -(M1 - @i*M2)/2 - @i*(L1 - @i*L2)/2;
H1 = L3 - @i*M3;
E2 = (-P0 - P3)/2;
F2 = (K0 - K3)/2;
H2 = -(D + L3);
E3 = -(M1 - @i*M2)/2 + @i*(L1 - @i*L2)/2;
F3 = (M1 + @i*M2)/2 - @i*(L1 + @i*L2)/2;
H3 = L3 + @i*M3;
E4 = dp_weyl_matmul(E1,E2) - dp_weyl_matmul(E2,E1);
E5 = dp_weyl_matmul(E2,E3) - dp_weyl_matmul(E3,E2);
E6 = dp_weyl_matmul(E1,E5) - dp_weyl_matmul(E5,E1);

/* parts of sl(4) relation */
assert( dp_weyl_matmul(E1,F1) - dp_weyl_matmul(F1,E1)==H1 );
assert( dp_weyl_matmul(E2,F2) - dp_weyl_matmul(F2,E2)==H2 );
assert( dp_weyl_matmul(E3,F3) - dp_weyl_matmul(F3,E3)==H3 );
assert( dp_weyl_matmul(E2,F1)==dp_weyl_matmul(F1,E2) );
assert( dp_weyl_matmul(H1,E1) - dp_weyl_matmul(E1,H1)==2*E1 );
assert( dp_weyl_matmul(H2,E2) - dp_weyl_matmul(E2,H2)==2*E2 );


assert( appdiff(-D0*S0+D1*S1+D2*S2+D3*S3 , Psi , V)==0*Psi );   /* Weyl equation */
assert( appdiff(H1,Psi,V)==0*Psi );        /* j1==0 */
assert( appdiff(H2,Psi,V)-2*Psi==0*Psi );  /* d+j1+j2==2 */
assert( appdiff(H3,Psi,V)+Psi==0*Psi );    /* j2==1/2 */
assert( appdiff(F1,Psi,V)==0*Psi );
assert( appdiff(F2,Psi,V)==0*Psi );
assert( appdiff(F3,Psi,V)==0*Psi );
assert( appdiff(E1,Psi,V)==0*Psi );
assert( appdiff(E3,appdiff(E3,Psi,V),V)==0*Psi );
assert( appdiff(E5,Psi,V)+appdiff(E2,appdiff(E3,Psi,V),V)==0*Psi );
assert( appdiff(E6 , appdiff(E2,Psi,V) , V) - appdiff(E4 , appdiff(E5,Psi,V) , V)==0*Psi);


/**** "spin 1" fields ****/
T0 = dp_ptod(1,V)*newmat(3,3,[[1,0,0],[0,1,0],[0,0,1]]);
T1 = dp_ptod(1,V)*@i*newmat(3,3,[[0,-2,0],[-1,0,-1],[0,-2,0]])/2;
T2 = dp_ptod(1,V)*@i*newmat(3,3,[[0,2*@i,0],[-@i,0,@i],[0,-2*@i,0]])/2;
T3 = dp_ptod(1,V)*@i*newmat(3,3,[[-2,0,0],[0,0,0],[0,0,2]])/2;

/* (T1,T2,T3) is a representaion of so(3) */
assert(T1*T2-T2*T1==T3);
assert(T2*T3-T3*T2==T1);
assert(T3*T1-T1*T3==T2);


/**** helicity +1 ****/
Phi2 = newvect(3,[(x1-@i*x2)^2 , (x0-x3)*(x1-@i*x2) , (x0-x3)^2])/(x0^2-x1^2-x2^2-x3^2)^3;
DS = dp_weyl_mul(X0,D0) + dp_weyl_mul(X1,D1) + dp_weyl_mul(X2,D2) + dp_weyl_mul(X3,D3);

D = (DS + 2)*T0;
P0 = D0*T0;
P1 = D1*T0;
P2 = D2*T0;
P3 = D3*T0;
L1 = (dp_weyl_mul(X0,D1) + dp_weyl_mul(X1,D0))*T0 - @i*T1;
L2 = (dp_weyl_mul(X0,D2) + dp_weyl_mul(X2,D0))*T0 - @i*T2;
L3 = (dp_weyl_mul(X0,D3) + dp_weyl_mul(X3,D0))*T0 - @i*T3;
M1 = (dp_weyl_mul(X2,D3) - dp_weyl_mul(X3,D2))*T0 - T1;
M2 = (dp_weyl_mul(X3,D1) - dp_weyl_mul(X1,D3))*T0 - T2;
M3 = (dp_weyl_mul(X1,D2) - dp_weyl_mul(X2,D1))*T0 - T3;
K0 = (2*dp_weyl_mul(X0,DS) - dp_weyl_mul(X0*X0-X1*X1-X2*X2-X3*X3,D0))*T0 + 4*X0*T0 - 2*@i*(X1*T1+X2*T2+X3*T3);
K1 = (-2*dp_weyl_mul(X1,DS) - dp_weyl_mul(X0*X0-X1*X1-X2*X2-X3*X3,D1))*T0 - 4*X1*T0 + 2*@i*X0*T1 + 2*X2*T3 - 2*X3*T2;
K2 = (-2*dp_weyl_mul(X2,DS) - dp_weyl_mul(X0*X0-X1*X1-X2*X2-X3*X3,D2))*T0 - 4*X2*T0 + 2*@i*X0*T2 - 2*X1*T3 + 2*X3*T1;
K3 = (-2*dp_weyl_mul(X3,DS) - dp_weyl_mul(X0*X0-X1*X1-X2*X2-X3*X3,D3))*T0 - 4*X3*T0 + 2*X1*T2 - 2*X2*T1 + 2*@i*X0*T3;

E1 = (M1 + @i*M2)/2 + @i*(L1 + @i*L2)/2;
F1 = -(M1 - @i*M2)/2 - @i*(L1 - @i*L2)/2;
H1 = L3 - @i*M3;
E2 = (-P0 - P3)/2;
F2 = (K0 - K3)/2;
H2 = -(D + L3);
E3 = -(M1 - @i*M2)/2 + @i*(L1 - @i*L2)/2;
F3 = (M1 + @i*M2)/2 - @i*(L1 + @i*L2)/2;
H3 = L3 + @i*M3;
E4 = dp_weyl_matmul(E1,E2) - dp_weyl_matmul(E2,E1);
E5 = dp_weyl_matmul(E2,E3) - dp_weyl_matmul(E3,E2);
E6 = dp_weyl_matmul(E1,E5) - dp_weyl_matmul(E5,E1);

/* parts of sl(4) relation */
assert( dp_weyl_matmul(E1,F1) - dp_weyl_matmul(F1,E1)==H1 );
assert( dp_weyl_matmul(E2,F2) - dp_weyl_matmul(F2,E2)==H2 );
assert( dp_weyl_matmul(E3,F3) - dp_weyl_matmul(F3,E3)==H3 );
assert( dp_weyl_matmul(E2,F1)==dp_weyl_matmul(F1,E2) );
assert( dp_weyl_matmul(H1,E1) - dp_weyl_matmul(E1,H1)==2*E1 );
assert( dp_weyl_matmul(H2,E2) - dp_weyl_matmul(E2,H2)==2*E2 );


/* SD-Maxwell equation */
assert( appdiff(D0*S0+D1*S1+D2*S2+D3*S3 , newvect(2,[Phi2[0],Phi2[1]]) , V)==newvect(2,[0,0]) );
assert( appdiff(D0*S0+D1*S1+D2*S2+D3*S3 , newvect(2,[Phi2[1],Phi2[2]]) , V)==newvect(2,[0,0]) );

assert( appdiff(H1,Phi2,V)+2*Phi2==0*Phi2 );
assert( appdiff(H2,Phi2,V)-3*Phi2==0*Phi2 );
assert( appdiff(H3,Phi2,V)==0*Phi2 );
assert( appdiff(F1,Phi2,V)==0*Phi2 );
assert( appdiff(F2,Phi2,V)==0*Phi2 );
assert( appdiff(F3,Phi2,V)==0*Phi2 );
assert( appdiff(E1,appdiff(E1,appdiff(E1,Phi2,V),V),V)==0*Phi2 );
assert( appdiff(E3,Phi2,V)==0*Phi2 );
assert( 2*appdiff(E4,Phi2,V)-appdiff(E2,appdiff(E1,Phi2,V),V)==0*Phi2 );
assert( appdiff(E6 , appdiff(E2,Phi2,V) , V) - appdiff(E4 , appdiff(E5,Phi2,V) , V)==0*Phi2);


/**** helicity -1 ****/
Psi2 = newvect(3,[(-x0+x3)^2 , (x1+@i*x2)*(-x0+x3) , (x1+@i*x2)^2])/(x0^2-x1^2-x2^2-x3^2)^3;
DS = dp_weyl_mul(X0,D0) + dp_weyl_mul(X1,D1) + dp_weyl_mul(X2,D2) + dp_weyl_mul(X3,D3);

D = (DS + 2)*T0;
P0 = D0*T0;
P1 = D1*T0;
P2 = D2*T0;
P3 = D3*T0;
L1 = (dp_weyl_mul(X0,D1) + dp_weyl_mul(X1,D0))*T0 + @i*T1;
L2 = (dp_weyl_mul(X0,D2) + dp_weyl_mul(X2,D0))*T0 + @i*T2;
L3 = (dp_weyl_mul(X0,D3) + dp_weyl_mul(X3,D0))*T0 + @i*T3;
M1 = (dp_weyl_mul(X2,D3) - dp_weyl_mul(X3,D2))*T0 - T1;
M2 = (dp_weyl_mul(X3,D1) - dp_weyl_mul(X1,D3))*T0 - T2;
M3 = (dp_weyl_mul(X1,D2) - dp_weyl_mul(X2,D1))*T0 - T3;
K0 = (2*dp_weyl_mul(X0,DS) - dp_weyl_mul(X0*X0-X1*X1-X2*X2-X3*X3,D0))*T0 + 4*X0*T0 + 2*@i*(X1*T1+X2*T2+X3*T3);
K1 = (-2*dp_weyl_mul(X1,DS) - dp_weyl_mul(X0*X0-X1*X1-X2*X2-X3*X3,D1))*T0 - 4*X1*T0 - 2*@i*X0*T1 + 2*X2*T3 - 2*X3*T2;
K2 = (-2*dp_weyl_mul(X2,DS) - dp_weyl_mul(X0*X0-X1*X1-X2*X2-X3*X3,D2))*T0 - 4*X2*T0 - 2*@i*X0*T2 - 2*X1*T3 + 2*X3*T1;
K3 = (-2*dp_weyl_mul(X3,DS) - dp_weyl_mul(X0*X0-X1*X1-X2*X2-X3*X3,D3))*T0 - 4*X3*T0 + 2*X1*T2 - 2*X2*T1 - 2*@i*X0*T3;


E1 = (M1 + @i*M2)/2 + @i*(L1 + @i*L2)/2;
F1 = -(M1 - @i*M2)/2 - @i*(L1 - @i*L2)/2;
H1 = L3 - @i*M3;
E2 = (-P0 - P3)/2;
F2 = (K0 - K3)/2;
H2 = -(D + L3);
E3 = -(M1 - @i*M2)/2 + @i*(L1 - @i*L2)/2;
F3 = (M1 + @i*M2)/2 - @i*(L1 + @i*L2)/2;
H3 = L3 + @i*M3;
E4 = dp_weyl_matmul(E1,E2) - dp_weyl_matmul(E2,E1);
E5 = dp_weyl_matmul(E2,E3) - dp_weyl_matmul(E3,E2);
E6 = dp_weyl_matmul(E1,E5) - dp_weyl_matmul(E5,E1);

/* parts of sl(4) relation */
assert( dp_weyl_matmul(E1,F1) - dp_weyl_matmul(F1,E1)==H1 );
assert( dp_weyl_matmul(E2,F2) - dp_weyl_matmul(F2,E2)==H2 );
assert( dp_weyl_matmul(E3,F3) - dp_weyl_matmul(F3,E3)==H3 );
assert( dp_weyl_matmul(E2,F1)==dp_weyl_matmul(F1,E2) );
assert( dp_weyl_matmul(H1,E1) - dp_weyl_matmul(E1,H1)==2*E1 );
assert( dp_weyl_matmul(H2,E2) - dp_weyl_matmul(E2,H2)==2*E2 );

/* ASD Maxwell equation */
assert( appdiff(-D0*S0+D1*S1+D2*S2+D3*S3 , newvect(2,[Psi2[0],Psi2[1]]) , V)==newvect(2,[0,0]) );
assert( appdiff(-D0*S0+D1*S1+D2*S2+D3*S3 , newvect(2,[Psi2[1],Psi2[2]]) , V)==newvect(2,[0,0]) );

assert( appdiff(H1,Psi2,V)==0*Psi2 );
assert( appdiff(H2,Psi2,V)-3*Psi2==0*Psi2 );
assert( appdiff(H3,Psi2,V)+2*Psi2==0*Psi2 );
assert( appdiff(F1,Psi2,V)==0*Psi2 );
assert( appdiff(F2,Psi2,V)==0*Psi2 );
assert( appdiff(F3,Psi2,V)==0*Psi2 );
assert( appdiff(E1,Psi2,V)==0*Psi2 );
assert( appdiff(E3,appdiff(E3,appdiff(E3,Psi2,V),V),V)==0*Psi2 );
assert( 2*appdiff(E5,Psi2,V)+appdiff(E2,appdiff(E3,Psi2,V),V)==0*Psi2 );
assert( appdiff(E6 , appdiff(E2,Psi2,V) , V) - appdiff(E4 , appdiff(E5,Psi2,V) , V)==0*Psi2);


2015-04-14

[]2度目のCubical

Thierry Coquandが、cubical interpreterでHigher Inductive Typeが使えるようになったとアナウンスしていた。理論的な面で進歩があったわけではないようだけど、cubical interpreterの実装が少し進歩したらしい。

https://github.com/mortberg/cubicaltt


基本的には、以前の実装を踏襲しているけど、primitive扱いだったreflとかmapOnPathが

refl (A : U) (a : A) : Id A a a = <i> a

みたいな形で定義できるようになっている。前も、内部的には、こういう形で扱われていた(実装上はTermとValueの区別が、あんまりなくなった?)のだけど、明示的に書けるようになったので、cubical type theoryを理解するのに便利になった。この変更のおかげで、実装も以前より、すっきりしたように見える。


<i> aは直感的には、i=0とi=1にaを対応させるintervalのことと思う(一次元"cube")。

Variation on Cubical sets

http://www.cse.chalmers.se/~coquand/comp.pdf

の13節に、operational semanticsとかがある。一般に、p : Id A x yは、i=0にx、i=1にyを対応させる"関数のようなもの"と解釈されて、新しいcubical interpreterでは、p@0とp@1で、それぞれx,yを取り出せる。直感的には、iは区間[0,1]のことであるが、p@0.5とかはできない


Syntaxとかも少し変更されているようなので、

はじめてのCubical

http://d.hatena.ne.jp/m-a-o/20131227#p1

の後半の例を書き直してみた(Syntaxは、また変わるかもしれないが)

module firsttest where

-- definitions from examples/prelude.ctt
-- -- "IdP P a b" is a built-in notion where P : Id U A B, a : A and b : B
-- "<a>" abstracts the name/color/dimension a.
Id (A : U) (a0 a1 : A) : U = IdP (<a> A) a0 a1

refl (A : U) (a : A) : Id A a a = <i> a

mapOnPath (A B : U) (f : A -> B) (x y : A)
          (p : Id A x y) : Id B (f x) (f y) = <a> f (p @ a)

funExt (A : U) (B : A -> U) (f g : (x : A) -> B x)
       (p : (x : A) -> Id (B x) (f x) (g x)) :
       Id ((y : A) -> B y) f g = <i> \(a : A) -> (p a) @ i

subst (A : U) (P : A -> U) (a b : A) (p : Id A a b) (e : P a) : P b =
  transport (mapOnPath A U P a b p) e


-- end of quote


ext (A B:U) (f g:A->B) (p : (x:A) -> Id B (f x) (g x)) : Id (A->B) f g = funExt A (\(a:A) -> B) f g p

data Nat = zero | succ (n:Nat)

add : Nat -> Nat -> Nat = split
  zero -> \(y:Nat) -> y
  succ n -> \(y:Nat) -> succ (add n y)

add0 : Nat -> Nat = \(n:Nat) -> add zero n

id (A:U) (a:A) : A = a

add0 : Nat -> Nat = \(n:Nat) -> add zero n

add0_eq_id : Id (Nat->Nat) add0 (id Nat) = 
   ext Nat Nat add0 (id Nat) p
      where
         p (n:Nat) : (Id Nat (add0 n) (id Nat n)) = refl Nat n

test : Nat = subst (Nat->Nat) f add0 (id Nat) add0_eq_id  zero
  where
    f (x:Nat -> Nat) : U = Nat 

前半部分は、prelude.cttから取ってきたので、import preludeで済む(transportとIdPはprimitive)。ファイルをloadして、testを評価すると、これが、ちゃんとzeroであることが分かる


以前は、内部的に扱われていたKan compositionも明示的に書けるようになった(Variation on Cubical setsの4節に説明がある)。例えば

module kantest where

Id (A : U) (a0 a1 : A) : U = IdP (<a> A) a0 a1
kan (A : U) (a b c d : A) (p : Id A a b) (q : Id A a c)
                          (r : Id A b d) : Id A c d =
  <i> comp A (p @ i) [ (i = 0) -> q , (i = 1) -> r ]

これは

A model of type theory in cubical sets

http://www.cse.chalmers.se/~coquand/mod1.pdf

の操作$X^{+}$に対応している。<i> comp A (p @ i) [ (i = 0) -> q , (i = 1) -> r ]は、i=0の時

(<i> comp A (p @ i) [ (i = 0) -> q , (i = 1) -> r ])@0 = comp A a [(i=0)->q]

で、qの型はq : Id A a cであるけど、compは点a=q@0をq : Id A a cをc=q@1に"繋ぐ"ような働きをして、

(<i> comp A (p @ i) [ (i = 0) -> q , (i = 1) -> r ])@0 = c

となる。同様に、i=1の時は

(<i> comp A (p @ i) [ (i = 0) -> q , (i = 1) -> r ])@1 = d

で、<i> comp A (p @ i) [ (i = 0) -> q , (i = 1) -> r ]自体の型は、Id A c dとなる。これを使えば、Idの推移律を示すことができる。対称律も、本来は、Kan filling/Kan compositionを使って作るのだけど、もっと単純に

eqsym (A:U) (x y:A) (p:Id A x y) : Id A y x= <i> p @ -i

と書ける。-iは、0を1に、1を0にする



compはintervalだけじゃなく、squareとか、任意のn-cubeも作れる。

module sqtest where
Id (A:U) (x y:A) : U = IdP (<i> A) x y
prop (A : U) : U = (a b : A) -> IdP (<i> A) a b
--         u
--    a0 -----> a1
--    |         |
-- r0 |         | r1
--    |         |
--    V         V
--    b0 -----> b1
--         v

Square (A : U) (a0 a1 : A) (u : Id A a0 a1)
               (b0 b1 : A) (v : Id A b0 b1)
               (r0 : Id A a0 b0) (r1 : Id A a1 b1) : U
  = IdP (<i> (IdP (<j> A) (u @ i) (v @ i))) r0 r1

propSquare  (A:U) (H:prop A) (a a0 a1 b0 b1:A) (r0 : Id A a0 b0) (r1 : Id A a1 b1) (u : Id A a0 a1) (v : Id A b0 b1) 
  : Square A a0 a1 u b0 b1 v r0 r1
  = <i j> comp A a [(i=0)-> H a (r0@j) ,(i=1) -> H a (r1@j) , (j=0) -> H a (u@i) , (j=1) -> H a (v@i)]

とか。<i j> comp A a [(i=0)-> H a (r0@j) ,(i=1) -> H a (r1@j) , (j=0) -> H a (u@i) , (j=1) -> H a (v@i)]は、頂点(i,j)=(0,0)に於いては、

comp A a [(i=0) -> H a (r0@0)]

と同じで、これはa0である。または

comp A a [(j=0) -> H a (u@0)]

でもあり、やっぱり、a0である。これは、当然ながら、ちゃんとcompatibleでないといけない。[]の中は空でもよくて(0-cube相当?)、comp A a [ ]は、aになる。


glueingもあるけど略


Higher Inductive Typeの作り方は,circle.cttとかsusp.cttとかinterval.cttとかにある

data S1 = base
        | loop_ @ base ~ base

loop : Id S1 base base = <i> loop_{S1} @ i

みたいにやる。なんでだかよく分からないけど、コンストラクタloop_の型は、Id S1 base baseでなく、上のような形で作る必要がある。それ以外は、普通



Quotients of sets

http://ncatlab.org/nlab/show/higher+inductive+type#QuotientsOfSets

を真似て、HITを使ってquotient typeを作ろうと思ったのだけど、書いてある定義は出来ない気がする

module quotient where
import prelude

rel (A:U) : U = A -> A -> U

prop (A:U) : U = (x y:A) -> Id A x y

propRel (A:U) (R:rel A) : U = (x y : A) -> prop (R x y)

{-
data quotient (A:U) (R:rel A) (p:propRel A R) 
  = proj (x:A) 
  | relate (x y:A) (r:R x y) @ proj x ~ proj y
  | contr1 (x y:quotient A R p) (u v : Id (quotient A R p) x y) @ u ~ v
-}

data quotient (A:U) (R:rel A) 
  = proj (x:A) 
  | relate (x y:A) @ proj x ~ proj y

よく考えてみると、contr1は、商をhSetにするために必要なのだと思うけど、そもそも、AがhSetでない場合、商もhSetであるとは限らないで、contr1が定義出来てしまうと、矛盾してしまう気がするけど、どうなんだろう。

quotient type in cubical type theory

http://d.hatena.ne.jp/m-a-o/20140825#p1

で定義したquotと同等な定義は上の形でいいんじゃないかという気がする。推測では、これだけで、(RがpropRelかつequivalenceの時)、商の普遍性が言え("inductive"なんだから普遍性は満たすだろうと思うけど)て、HITを使わないquotient typeと等しいことが言えると思う(けど、証明は書いていない。そのうち書く)


あと、こういう↓感じで、assocを定義しようとすると無理っぽい。2回addが出てきているので、だめっぽい

module testdat where

-- type constructor "assoc" is invalid
data TestAlg (A:U)
  = empty
  | lift (x:A)
  | add (x:TestAlg A) (y:TestAlg A)
  | comm (x y:TestAlg A) @ add x y ~ add y x
--  | assoc (x y z:TestAlg A) @ add (add x y) z ~ add x (add y z)


homotopy groups of spheres

http://ncatlab.org/homotopytypetheory/show/homotopy+groups+of+spheres

Guillaume has proved that there exists an $n$ such that $\pi_4(S^3)$ is $\mathbb{Z}_n$. Given a computational interpretation, we could run this proof and check that $n$ is 2.

という一文があるのだけど、cubical interpreterは、univalence axiomのcomputational interpretationを与えていると(予想)され、かつHITを使えるようになった(Suspensionがあるので、S^3は定義できる)ので、試金石として試すことができるはず