⍝ Gauss-Jordan elimination:
try←{⍺←{⍵} ⍝ gauss_jordan vs primitive ⌹.
⎕ct←1e¯13 ⍝ slightly more tolerant comparison.
fuzz0←{↑⍺⍺/⍺ ⍵+⎕ct>|⍺ ⍵} ⍝ 0-fuzzy comparison.
(⍺⌹⍵)≡fuzz0 ⍺ gauss_jordan ⍵
}
⍝ The dfns.dws test scripts are used in Dyalog's overnight QA check to prove
⍝ the integrity of each version of the interpreter on each supported hardware
⍝ platform. In the above test, comparison around 0 has been made slightly more
⍝ tolerant to accommodate small differences in compiler optimisation in differ-
⍝ ent versions.
⍝ monadic ---------------------------
try ↑(1 1 2)(¯1 ¯2 3)(3 ¯7 4)
1
try ↑(4 8 4 0)(1 4 7 2)(1 5 4 ¯3)(1 3 0 ¯2)
1
try ↑(1 2 3)(¯3 1 5)(2 4 ¯1)
1
⍝ dyadic ----------------------------
(3 4⍴⍳12) try ↑(1 1 2)(¯1 ¯2 3)(3 ¯7 4)
1
(3 4⍴⍳12) try ↑(1 2 3)(¯3 1 5)(2 4 ¯1)
1
(4 5⍴⍳12) try ↑(4 8 4 0)(1 4 7 2)(1 5 4 ¯3)(1 3 0 ¯2)
1
hil←{÷1+∘.+⍨(⍳⍵)-⎕io} ⍝ order ⍵ Hilbert matrix.
0 1 0 disp gauss_jordan∘hil¨ 0 to 5
┌─┬─┬─────┬─────────────┬──────────────────────┬───────────────────────────────────┐
│0│1│ 4 ¯6│ 9 ¯36 30│ 16 ¯120 240 ¯140│ 25 ¯300 1050 ¯1400 630│
│ │ │¯6 12│¯36 192 ¯180│¯120 1200 ¯2700 1680│ ¯300 4800 ¯18900 26880 ¯12600│
│ │ │ │ 30 ¯180 180│ 240 ¯2700 6480 ¯4200│ 1050 ¯18900 79380 ¯117600 56700│
│ │ │ │ │¯140 1680 ¯4200 2800│¯1400 26880 ¯117600 179200 ¯88200│
│ │ │ │ │ │ 630 ¯12600 56700 ¯88200 44100│
└─┴─┴─────┴─────────────┴──────────────────────┴───────────────────────────────────┘
⍝ We can modify the function to use rational numbers for perfect precision:
gauss_jordan_rats←{⎕ML ⎕IO←0 ⍝ Rational Gauss-Jordan.
elim←{ ⍝ elimination of row/col ⍺.
mat←⍵ ⍝ name matrix for updating.
v←⍺+{⍵⍳⌈/⍵}{×/*⌿⊃⌽⍵}¨⍺↓mat[;⍺] ⍝ index of pivot row.
mat[⍺ v;]←mat[v ⍺;] ⍝ exchange ⍺th and pth rows.
mat[⍺;]←mat[⍺;]q¨mat[⍺;⍺] ⍝ reduce col diagonal to 1.
mul←mat[;⍺]p¨¯2↑¨⍺≠⍳⊃⍴⍵ ⍝ multiplier.
mat d¨mul∘.p mat[⍺;] ⍝ reduce col off-diagonals to 0.
}
0::⎕SIGNAL ⎕EN ⍝ pass back error to caller.
⍺←{0,⊂⍵ rats 1}¨=/↑⍳⍴⍵ ⍝ id matrix for monadic case.
p←×rats prodquot ⍝ rational product.
q←÷rats prodquot ⍝ .. quotient.
d←+rats diff(-rats) ⍝ .. difference.
(⍴⍺)⍴(0 1×⍴⍵)↓↑elim/(⌽⍳⌊/⍴⍵),⊂⍵,⍺ ⍝ elimination/ ··· 2 1 0 (⍵,⍺)
}
prodquot←{ ⍝ product/quotient.
(s a)(t b)←⍺ ⍵ ⍝ signs and fractions.
((s≠t)>0∊a,b)(a ⍺⍺ b) ⍝ rule-of-signs.
}
diff←{ ⍝ difference.
(s a)(t b)←⍺ ⍵ ⍝ signs and fractions.
s≠t:s(a ⍺⍺ b) ⍝ differing signs:
≥/{×/*⌿⍵}¨a b:s(a ⍵⍵ b) ⍝ ⍺≥⍵: ⍺-⍵
(~s)(b ⍵⍵ a) ⍝ ⍺<⍵: -⍵-⍺
}
hilr←{ ⍝ rational ⍵-Hilbert matrix.
{0 ⍵}∘{1 rats ⍵}¨1+∘.+⍨(⍳⍵)-⎕IO
}
fmt←{ ⍝ ⍺-separated rational format.
↑⍺{ ⍝ ⍺ is separator.
(⍺/'¯'),↑⍺⍺{⍺,⍺⍺,⍵}/{ ⍝ sign, unsigned fraction.
⍕↑×nats/*nats⌿⍵,1 ⍝ using nats for precision.
}¨|(1 0=⊂↑{⍵>0}/↓⍵)/¨⊂⍵ ⍝ numerator/denominator.
}/⍵ ⍝ sign/fraction pair.
}
align←{ ⍝ ⍺-aligned columns.
⍵,⍨¨({(⌈⌿⍵)-[⊃⌽⍳⍴⍴⍵]⍵}⍵⍳¨⍺)/¨' '
}
show←'/'∘align∘('/'∘fmt¨) ⍝ display rational matrix.
show hilr 10 ⍝ 10-Hilbert.
1/1 1/2 1/3 1/4 1/5 1/6 1/7 1/8 1/9 1/10
1/2 1/3 1/4 1/5 1/6 1/7 1/8 1/9 1/10 1/11
1/3 1/4 1/5 1/6 1/7 1/8 1/9 1/10 1/11 1/12
1/4 1/5 1/6 1/7 1/8 1/9 1/10 1/11 1/12 1/13
1/5 1/6 1/7 1/8 1/9 1/10 1/11 1/12 1/13 1/14
1/6 1/7 1/8 1/9 1/10 1/11 1/12 1/13 1/14 1/15
1/7 1/8 1/9 1/10 1/11 1/12 1/13 1/14 1/15 1/16
1/8 1/9 1/10 1/11 1/12 1/13 1/14 1/15 1/16 1/17
1/9 1/10 1/11 1/12 1/13 1/14 1/15 1/16 1/17 1/18
1/10 1/11 1/12 1/13 1/14 1/15 1/16 1/17 1/18 1/19
show gauss_jordan_rats hilr 10 ⍝ inverse of 10-Hilbert.
100/1 ¯4950/1 79200/1 ¯600600/1 2522520/1 ¯6306300/1 9609600/1 ¯8751600/1 4375800/1 ¯923780/1
¯4950/1 326700/1 ¯5880600/1 47567520/1 ¯208107900/1 535134600/1 ¯832431600/1 770140800/1 ¯389883780/1 83140200/1
79200/1 ¯5880600/1 112907520/1 ¯951350400/1 4281076800/1 ¯11237826600/1 17758540800/1 ¯16635041280/1 8506555200/1 ¯1829084400/1
¯600600/1 47567520/1 ¯951350400/1 8245036800/1 ¯37875637800/1 101001700800/1 ¯161602721280/1 152907955200/1 ¯78843164400/1 17071454400/1
2522520/1 ¯208107900/1 4281076800/1 ¯37875637800/1 176752976400/1 ¯477233036280/1 771285715200/1 ¯735869534400/1 382086104400/1 ¯83223340200/1
¯6306300/1 535134600/1 ¯11237826600/1 101001700800/1 ¯477233036280/1 1301544644400/1 ¯2121035716800/1 2037792556800/1 ¯1064382719400/1 233025352560/1
9609600/1 ¯832431600/1 17758540800/1 ¯161602721280/1 771285715200/1 ¯2121035716800/1 3480673996800/1 ¯3363975014400/1 1766086882560/1 ¯388375587600/1
¯8751600/1 770140800/1 ¯16635041280/1 152907955200/1 ¯735869534400/1 2037792556800/1 ¯3363975014400/1 3267861442560/1 ¯1723286307600/1 380449555200/1
4375800/1 ¯389883780/1 8506555200/1 ¯78843164400/1 382086104400/1 ¯1064382719400/1 1766086882560/1 ¯1723286307600/1 912328045200/1 ¯202113826200/1
¯923780/1 83140200/1 ¯1829084400/1 17071454400/1 ¯83223340200/1 233025352560/1 ¯388375587600/1 380449555200/1 ¯202113826200/1 44914183600/1
Back to: code
Back to: Workspaces
Trouble seeing APL font?