program ioi94day2prb1invert(input, output);
{ Tom Verhoeff, Eindhoven University of Technology }
{ Pre-processor to invert a system of linear algebraic equations }

const
  Test = true ;

procedure Compute ;
  var
    A: array[1..9, 1..9] of integer; { matrix of coefficients }
    B: array [1..9, 1..9] of integer; { right-hand side matrix }
    { The system of linear algebraic equations to solve: }
    {   (Sum j: j in [1..9]: A[i,j] * T[j,k]) = B[i,k]  for i in [1..9] }

  procedure WriteMatrix ; { for testing }
    var i, j: integer;
    begin
    for i := 1 to 9 do begin
      for j := 1 to 9 do write(A[i, j]:2) ;
      write('  | ') ;
      for j := 1 to 9 do write(B[i, j]:2) ;
      writeln
      end { for i } ;
    writeln
    end { WriteMatrix } ;

  procedure SetUp ;
    var f: text; i, j: integer;
    begin
    assign(f, 'matrix.dat') ;
    reset(f) ;
    for i := 1 to 9 do begin
      for j := 1 to 9 do read(f, A[i, j]) ;
      readln(f) ;
      end { for i } ;
    { set up identity matrix in in B[1..9] }
    for i := 1 to 9 do begin
      for j := 1 to 9 do
        B[i, j] := ord(i = j)
      end { for i } ;
    if Test then begin
      writeln('Initial set up is:') ;
      WriteMatrix
      end { if }
    end { SetUp } ;

  procedure Solve ; { by Gauss-Jordan elimination }
    var h, i, j, k: integer;
    begin
    { transform A into the unit matrix and B into the inverse matrix }
    for i := 1 to 9 do begin { process column i }
      { find pivot by bounded linear search for first 1 or 3 in column i }
      k := i ; h := 10 ;
      while k <> h do
        if A[k, i] in [1, 3] then h := k else k := succ(k) ;
      if k = 10 then begin
        writeln('Pivot not found in step ', i:1) ;
        halt
        end { if } ;
      if Test then writeln('Pivot for column ', i:1, ' found in row ', k:1) ;
      { swap rows i and k }
      for j := i to 9 do begin
        h := A[i, j] ; A[i, j] := A[k, j] ; A[k, j] := h
        end { for j } ;
      for j := 1 to 9 do begin
        h := B[i, j] ; B[i, j] := B[k, j] ; B[k, j] := h
        end { for j } ;
      { normalize row i, yielding A[i, i] = 1 }
      h := A[i, i] ; { h * A[i, i] = 1 (mod 4) }
      for j := i to 9 do
        A[i, j] := (h * A[i, j]) mod 4 ;
      for j := 1 to 9 do
        B[i, j] := (h * B[i, j]) mod 4 ;
      { sweep column i to zeroes in rows other than i }
      for k := 1 to 9 do begin { take care of row k }
        if k <> i then begin
          h := 3*A[k, i] ; { -1 = 3 (mod 4) }
          for j := i to 9 do
            A[k, j] := (A[k, j] + h*A[i, j]) mod 4 ;
          for j := 1 to 9 do
            B[k, j] := (B[k, j] + h*B[i, j]) mod 4 ;
          end { if }
        end { for j } ;
      if Test then begin
        writeln('Situation after step ', i:1) ;
        WriteMatrix
        end { if }
      end { for i }
    end { Solve } ;

  procedure WriteInverse ;
    var f: text; i, j: integer;
    begin
    assign(f, 'inverse.dat') ;
    rewrite(f) ;
    for i := 1 to 9 do begin
      for j := 1 to 9 do
        write(f, B[i, j]:2) ;
      writeln(f)
      end { for i } ;
    close(f)
    end { WriteInverse } ;

  begin { Compute }
  SetUp ;
  Solve ;
  WriteInverse
  end { Compute } ;

begin
writeln('Pre-processor to invert matrix') ;
Compute
end.