MODULE LR2x2 EXPORTS LR2x2, LR2x2Extras;
(*
  Hand-optimized implementation of the "LR2x2" and "LR2x2Extra" interface. 

  Created by J. Stolfi with contributions by R. L. W. Liesenfeld. *)
  
IMPORT LR2Rep, LR2x2Rep, LRMxN, Wr, Thread, Math;

(* LR2x2 *)

TYPE RowT = ARRAY [0..N-1] OF ElemT;

PROCEDURE Null(): T =
  BEGIN
    RETURN T{RowT{0.0d0, ..}, ..}
  END Null;
  
PROCEDURE Identity(): T =
  BEGIN
    RETURN T{
      RowT{1.0d0, 0.0d0},
      RowT{0.0d0, 1.0d0}
    }
  END Identity;
    
PROCEDURE ToReal(READONLY m: T): RealT =
  VAR rm: RealT;
  BEGIN
    rm[0,0] := FLOAT(m[0,0], REAL);
    rm[0,1] := FLOAT(m[0,1], REAL);
    rm[1,0] := FLOAT(m[1,0], REAL);
    rm[1,1] := FLOAT(m[1,1], REAL);
    RETURN rm
  END ToReal;
  
PROCEDURE ToLongReal(READONLY m: T): LongRealT =
  BEGIN
    RETURN m
  END ToLongReal;
    
PROCEDURE FromReal(READONLY c: RealT): T =
  VAR m: T;
  BEGIN
    m[0,0] := FLOAT(c[0,0], ElemT);
    m[0,1] := FLOAT(c[0,1], ElemT);
    m[1,0] := FLOAT(c[1,0], ElemT);
    m[1,1] := FLOAT(c[1,1], ElemT);
    RETURN m
  END FromReal;

PROCEDURE FromLongReal(READONLY c: LongRealT): T =
  BEGIN
    RETURN c
  END FromLongReal;

PROCEDURE Get(READONLY m: T; i, j: [0..N-1]): LONGREAL =
  BEGIN
    RETURN m[i, j]
  END Get;
    
PROCEDURE Set(VAR m: T; i, j: [0..N-1]; x: LONGREAL) =
  BEGIN
    m[i, j] := x
  END Set;

PROCEDURE MapRow(READONLY a: LR2Rep.T; READONLY m: T): LR2Rep.T =
  BEGIN
    WITH
      a0 = a[0],
      a1 = a[1],
      
      m00 = m[0,0], m01 = m[0,1],
      m10 = m[1,0], m11 = m[1,1]
    DO
      RETURN LR2Rep.T{a0*m00 + a1*m10, a0*m01 + a1*m11}
    END;
  END MapRow;

PROCEDURE MapCol(READONLY m: T; READONLY a: LR2Rep.T): LR2Rep.T =
  BEGIN
    WITH
      m00 = m[0,0], m01 = m[0,1],
      m10 = m[1,0], m11 = m[1,1],
      
      a0 = a[0],
      a1 = a[1]
    DO
      RETURN LR2Rep.T{m00*a0 + m01*a1, m10*a0 + m11*a1}
    END;
  END MapCol;

PROCEDURE Mul(READONLY m, n: T): T =
  BEGIN
    WITH
      m00 = m[0,0], m01 = m[0,1],
      m10 = m[1,0], m11 = m[1,1],
      
      n00 = n[0,0], n01 = n[0,1],
      n10 = n[1,0], n11 = n[1,1]
    DO
      RETURN T{
        LR2Rep.T{m00*n00 + m01*n10, m00*n01 + m01*n11},
        LR2Rep.T{m10*n00 + m11*n10, m10*n01 + m11*n11}
      }
    END
  END Mul;
  
PROCEDURE Adj(READONLY m: T): T =
  BEGIN
    WITH
      m00 = m[0,0], m01 = m[0,1],
      m10 = m[1,0], m11 = m[1,1]
    DO
      RETURN T{
        LR2Rep.T{+ m11, - m01},
        LR2Rep.T{- m10, + m00}
      }
    END
  END Adj;

PROCEDURE Inv(READONLY m: T): T =
  BEGIN
    WITH
      m00 = m[0,0], m01 = m[0,1],
      m10 = m[1,0], m11 = m[1,1],
      
      d = + m00 * m11 - m01 * m10
    DO
      RETURN T{
        LR2Rep.T{+ m11/d, - m01/d},
        LR2Rep.T{- m10/d, + m00/d}
      }
    END;
  END Inv;

PROCEDURE Det(READONLY m: T): LONGREAL =
  BEGIN
    WITH
      m00 = m[0,0], m01 = m[0,1],
      m10 = m[1,0], m11 = m[1,1],
      
      d = + m00 * m11 - m01 * m10
    DO
      RETURN d
    END;
  END Det;  

PROCEDURE Print(
    wr: Wr.T; 
    READONLY m: T; 
    olp := "( "; osep := "\n  "; orp := " )";
    ilp := "("; isep := " "; irp := ")";
    fmt: PROCEDURE (x: ElemT): TEXT := LR2x2Rep.DefaultElemFmt;
    align: BOOLEAN := TRUE;
  ) RAISES {Wr.Failure, Thread.Alerted} =
  BEGIN
    LRMxN.Print(wr, m, ilp, olp, osep, orp, isep, irp, fmt, align)
  END Print;

PROCEDURE ToText(READONLY m: T): TEXT =
  BEGIN
    RETURN LRMxN.ToText(m)
  END ToText;

(* LR2x2Extra *)

PROCEDURE Rotation(angle: LONGREAL; scale: LONGREAL): T =
  BEGIN
    WITH 
      ct = scale * Math.cos(angle),
      st = scale * Math.sin(angle)
    DO
      RETURN T{
        LR2Rep.T{+ct, +st},
        LR2Rep.T{-st, +ct}
      }
    END
  END Rotation;

BEGIN 
END LR2x2.