max msp gen~で連立方程式を解く

ガウスジョルダン法を用いてmax/msp gen~で連立方程式を解きます。

ガウス・ジョルダン法

【簡単】線形代数のガウス・ジョルダン法って何?【やり方教えます】

以下のように係数を入力すると解をオーディオで出力してくれます。
2x-4y+6z=6
-x+7y+8z=-3
x+y-2z=2
の三元連立方程式を解いて
x=2.1724…
y=-0.1379…
z=0.01724…
と出力しています。

以下が連立方程式の解を求めるサイトで計算した結果です。
ほぼ同じような解になっています。

以下がgen~の中身です。ここではピボット選択無しのバージョンでガウスジョルダン法を書いています。

以下がgen~の中身です。コピペすると使えます。

Data a(5,5);

poke(a,in1,0,0);poke(a,in2,0,1); poke(a,in3,0,2);  poke(a,in4,0,3);
poke(a,in5,1,0);poke(a,in6,1,1); poke(a,in7,1,2);  poke(a,in8,1,3);
poke(a,in9,2,0);poke(a,in10,2,1);poke(a,in11,2,2);poke(a,in12,2,3);

n=3;q=0;i=0;g=0;j=0;k=0;

for(i=0;i<n;i+=1){
 q=peek(a,i,i);
 for(j=0;j<n+1;j+=1){
   poke(a,(peek(a,i,j)/q),i,j);
 }
 for(k=0;k<n;k+=1){
   g=peek(a,k,i);
   if(k != i){
    for(j=0;j<n+1;j+=1){
      poke(a,peek(a,k,j)-g*peek(a,i,j),k,j);
    }
   }
 } 
}


out1=peek(a,0,0);out2=peek(a,0,1);out3=peek(a,0,2);out4=peek(a,0,3);
out5=peek(a,1,0);out6=peek(a,1,1);out7=peek(a,1,2);out8=peek(a,1,3);
out9=peek(a,2,0);out10=peek(a,2,1);out11=peek(a,2,2);out12=peek(a,2,3);  

以下がガウスジョルダン法　ピボット選択有り　解の判別有り版です

以下がgen~

以下がgen~の中身です。

Data a(5,5);

poke(a,in1,0,0);poke(a,in2,0,1); poke(a,in3,0,2);  poke(a,in4,0,3);
poke(a,in5,1,0);poke(a,in6,1,1); poke(a,in7,1,2);  poke(a,in8,1,3);
poke(a,in9,2,0);poke(a,in10,2,1);poke(a,in11,2,2);poke(a,in12,2,3);

n=3;  //dimension 
m=1; 
epsl=0.00000001;isw=0; i=0;j=0;
nplsm=n+m; nn=0; ip=0; p=0; w=0;
for(nn=0;nn<n;nn+=1){
 p=0;
 for(i=nn;i<n;i+=1){
   if(p<abs(peek(a,i,nn))){
     p= peek(a,i,nn) ;
     p=abs(p);
     ip=i;
   }
 }
 if(p<=epsl){
     isw=1;
 }
 for(j=nn;j<=nplsm;j+=1){
   w=peek(a,nn,j);
   poke(a,peek(a,ip,j),nn,j);
   poke(a,w,ip,j);
 }
 for(j=nn+1;j<nplsm;j+=1){
    poke(a, peek(a,nn,j)/peek(a,nn,nn),nn,j);
 }
 for(i=0;i<n;i+=1 ){
   if(i != nn){
     for(j=nn+1;j<nplsm;j+=1){
       poke(a,peek(a,i,j)-peek(a,i,nn)*peek(a,nn,j),i,j);
     }
   }
 }
}

//out1=peek(a,0,0);out2=peek(a,0,1);out3=peek(a,0,2);out4=peek(a,0,3);
//out5=peek(a,1,0);out6=peek(a,1,1);out7=peek(a,1,2);out8=peek(a,1,3);
//out9=peek(a,2,0);out10=peek(a,2,1);out11=peek(a,2,2);out12=peek(a,2,3);

out1=peek(a,0,3);out2=peek(a,1,3);out3=peek(a,2,3);
out4=isw;