⍝ 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 compilers.

    ⍝ 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.

    gauss_jordan∘hil¨ 0 to 5
┌┬─┬─────┬─────────────┬──────────────────────┬───────────────────────────────────┐
││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←{⎕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 5                                ⍝ 5-Hilbert.
 1/1  1/2  1/3  1/4  1/5 
 1/2  1/3  1/4  1/5  1/6 
 1/3  1/4  1/5  1/6  1/7 
 1/4  1/5  1/6  1/7  1/8 
 1/5  1/6  1/7  1/8  1/9 

    ⍕show gauss_jordan_rats hilr 5              ⍝ inverse of 5-Hilbert.
    25/1    ¯300/1     1050/1    ¯1400/1     630/1 
  ¯300/1    4800/1   ¯18900/1    26880/1  ¯12600/1 
  1050/1  ¯18900/1    79380/1  ¯117600/1   56700/1 
 ¯1400/1   26880/1  ¯117600/1   179200/1  ¯88200/1 
   630/1  ¯12600/1    56700/1   ¯88200/1   44100/1 

⍝∇ gauss_jordan nats rats factors to

Back to: code

Back to: Workspaces