--
with TEXT_IO ;

package body RATIONAL_ARRAYS is

-- RATIONAL MATRIX ARITHMETIC PACKAGE BODY
--
--    Purpose :  To provide the basic set of RATIONAL matrix arithmetic
--               capability in Ada
--
--    Method :   This package defines the types
--                       rational_matrix  and  rational_vector
--               this package also defines the overloaded procedure
--                PUT  for the types  rational_matrix  and  rational_vector
--
--               this package then defines the overloaded operators
--
--                 +   -   *    unary +   unary -
--
--                 for combining types  rational_matrix  rational_vector  and  RATIONAL
--
--    WRITTEN BY : JON SQUIRE , 27 FEB 1983
--
--    COPYRIGHT  1983  WESTINGHOUSE ELECTRIC CORP.
--    COPYRIGHT  1984  WESTINGHOUSE ELECTRIC CORP. -- additions
--
--    for the type  rational_matrix  the first subscript is the row number,
--                               the second subscript is the column number
--

  function SQRT ( X : RATIONAL ) return RATIONAL is
    Y : RATIONAL := FRACTION ( 1 , 2 ) ;
  begin
    for I in 1 .. 10 loop
      Y := FRACTION ( 1 , 2 ) * ( Y + X / Y ) ;
    end loop ;
    return Y ;
  end SQRT ;

  function "+" ( A , B : RATIONAL_MATRIX ) return RATIONAL_MATRIX is
    C : RATIONAL_MATRIX ( A'FIRST( 1 ) .. A'LAST ( 1 ) , A'FIRST ( 2 ) .. A'
      LAST ( 2 )) ;
  begin

-- CHECK THAT MATRICES ARE COMPATIBLE FOR ADDING
    if A'LENGTH ( 1 ) /= B'LENGTH ( 1 ) or
       A'LENGTH ( 2 ) /= B'LENGTH ( 2 ) then
      raise MATRIX_ERROR ;
    end if ;
    if A'FIRST ( 1 ) = B'FIRST ( 1 ) and
       A'FIRST ( 2 ) = B'FIRST ( 2 ) then

-- SIMPLE INDEXING
      for I in A'RANGE ( 1 ) loop
        for J in A'RANGE ( 2 ) loop
          C ( I , J ) := A ( I , J ) + B ( I , J ) ;
        end loop ;
      end loop ;
    else

-- MATRICES HAVE DIFFERENT INDICES
      for I in A'RANGE ( 1 ) loop
        for J in A'RANGE ( 2 ) loop
          C ( I , J ) := A ( I , J ) + B ( I - A'FIRST( 1 ) + B'FIRST ( 1 ) ,
             J - A'FIRST ( 2 ) + B'FIRST ( 2 )) ;
        end loop ;
      end loop ;
    end if ;
    return C ;
  end ;

  function "-" ( A , B : RATIONAL_MATRIX ) return RATIONAL_MATRIX is
    C : RATIONAL_MATRIX ( A'FIRST( 1 ) .. A'LAST ( 1 ) , A'FIRST ( 2 ) .. A'
      LAST ( 2 )) ;
  begin

-- CHECK THAT MATRICES ARE COMPATIBLE FOR SUBTRACTING
    if A'LENGTH ( 1 ) /= B'LENGTH ( 1 ) or
       A'LENGTH ( 2 ) /= B'LENGTH ( 2 ) then
      raise MATRIX_ERROR ;
    end if ;
    if A'FIRST ( 1 ) = B'FIRST ( 1 ) and
       A'FIRST ( 2 ) = B'FIRST ( 2 ) then

-- SIMPLE INDEXING
      for I in A'RANGE ( 1 ) loop
        for J in A'RANGE ( 2 ) loop
          C ( I , J ) := A ( I , J ) - B ( I , J ) ;
        end loop ;
      end loop ;
    else

-- MATRICES HAVE DIFFERENT INDICES
      for I in A'RANGE ( 1 ) loop
        for J in A'RANGE ( 2 ) loop
          C ( I , J ) := A ( I , J ) - B ( I - A'FIRST( 1 ) + B'FIRST ( 1 ) ,
             J - A'FIRST ( 2 ) + B'FIRST ( 2 )) ;
        end loop ;
      end loop ;
    end if ;
    return C ;
  end ;

  function "*" ( A , B : RATIONAL_MATRIX ) return RATIONAL_MATRIX is
    C : RATIONAL_MATRIX ( A'FIRST( 1 ) .. A'LAST ( 1 ) , B'FIRST ( 2 ) .. B'
      LAST ( 2 )) ;
  begin
    if A'LENGTH ( 2 ) /= B'LENGTH ( 1 ) then
      raise MATRIX_ERROR ;
    end if ;
    for I in A'RANGE ( 1 ) loop
      for J in B'RANGE ( 2 ) loop
        C ( I , J ) := RATIONAL_ZERO ;
        for K in A'RANGE ( 2 ) loop
          C ( I , J ) := C ( I , J ) + A ( I , K ) * B ( K - A'FIRST( 2 ) + B
            'FIRST ( 1 ) , J) ;
        end loop ;
      end loop ;
    end loop ;
    return C ;
  end ;

  function "+" ( A , B : RATIONAL_VECTOR ) return RATIONAL_VECTOR is
    C : RATIONAL_VECTOR ( A'FIRST .. A'LAST ) ;
  begin
    if A'LENGTH /= B'LENGTH then
      raise MATRIX_ERROR ;
    end if ;
    for I in A'RANGE loop
      C ( I ) := A ( I ) + B ( I - A'FIRST + B'FIRST ) ;
    end loop ;
    return C ;
  end ;

  function "-" ( A , B : RATIONAL_VECTOR ) return RATIONAL_VECTOR is
    C : RATIONAL_VECTOR ( A'FIRST .. A'LAST ) ;
  begin
    if A'LENGTH /= B'LENGTH then
      raise MATRIX_ERROR ;
    end if ;
    for I in A'RANGE loop
      C ( I ) := A ( I ) - B ( I - A'FIRST + B'FIRST ) ;
    end loop ;
    return C ;
  end ;

  function "*" ( V : RATIONAL_VECTOR ;
                 A : RATIONAL_MATRIX ) return RATIONAL_VECTOR is
    CV : RATIONAL_VECTOR ( A'FIRST( 2 ) .. A'LAST ( 2 )) ;
    SUM : RATIONAL ;
  begin
    if A'LENGTH ( 1 ) /= V'LENGTH then
      raise MATRIX_ERROR ;
    end if ;
    for J in A'RANGE ( 2 ) loop
      SUM := RATIONAL_ZERO ;
      for I in A'RANGE ( 1 ) loop
        SUM := SUM + V ( I - A'FIRST( 1 ) + V'FIRST) * A ( I , J ) ;
      end loop ;
      CV ( J ) := SUM ;
    end loop ;
    return CV ;
  end ;

  function "*" ( A : RATIONAL_MATRIX ;
                 V : RATIONAL_VECTOR ) return RATIONAL_VECTOR is
    CV : RATIONAL_VECTOR ( A'FIRST( 1 ) .. A'LAST ( 1 )) ;
    SUM : RATIONAL ;
  begin
    if A'LENGTH ( 2 ) /= V'LENGTH then
      raise MATRIX_ERROR ;
    end if ;
    for I in A'RANGE ( 1 ) loop
      SUM := RATIONAL_ZERO ;
      for J in A'RANGE ( 2 ) loop
        SUM := SUM + V ( J - A'FIRST( 2 ) + V'FIRST) * A ( I , J ) ;
      end loop ;
      CV ( I ) := SUM ;
    end loop ;
    return CV ;
  end ;

  function DOT_PRODUCT ( P : RATIONAL_VECTOR ;
                         Q : RATIONAL_VECTOR ) return RATIONAL is
    F : RATIONAL := RATIONAL_ZERO ;
    TEMP : RATIONAL ;
  begin
    if P'LENGTH /= Q'LENGTH then
      raise MATRIX_ERROR ;
    end if ;
    for I in P'RANGE loop
      TEMP := P ( I ) * Q ( I - P'FIRST + Q'FIRST ) ;
      F := F + TEMP * TEMP ;
    end loop ;
    return SQRT ( F ) ;
  end DOT_PRODUCT ;

  function CROSS_PRODUCT ( P : RATIONAL_VECTOR ;
                           Q : RATIONAL_VECTOR ) return RATIONAL_VECTOR is
    CV : RATIONAL_VECTOR ( P'FIRST .. P'LAST ) ;
  begin
    if P'LENGTH /= 3 or
       Q'LENGTH /= 3 then
      raise MATRIX_ERROR ;
    end if ;
    CV ( P'FIRST ) := P ( P'FIRST + 1 ) * Q ( Q'FIRST + 2 ) - P ( P'FIRST + 2
       ) * Q ( Q'FIRST + 1 ) ;
    CV ( P'FIRST + 1 ) := P ( P'FIRST + 2 ) * Q ( Q'FIRST + 0 ) - P ( P'FIRST
       + 0 ) * Q ( Q'FIRST + 2 ) ;
    CV ( P'FIRST + 2 ) := P ( P'FIRST + 0 ) * Q ( Q'FIRST + 1 ) - P ( P'FIRST
       + 1 ) * Q ( Q'FIRST + 0 ) ;
    return CV ;
  end CROSS_PRODUCT ;

  function CROSS_PRODUCT ( A : RATIONAL_MATRIX ) return RATIONAL_VECTOR is
    CV : RATIONAL_VECTOR ( A'FIRST( 2 ) .. A'LAST ( 2 )) ;
    N : INTEGER := A'LENGTH ( 2 ) ;

-- RETURN N DIMENSIONAL VECTOR ORTHOGONAL TO N-1 VECTORS IN A
  begin
    if A'LENGTH ( 1 ) + 1 /= A'LENGTH ( 2 ) then
      raise MATRIX_ERROR ;
    end if ;
    return CV ;
  end CROSS_PRODUCT ;

  function IDENTITY_MATRIX ( N : INTEGER ) return RATIONAL_MATRIX is
    A : RATIONAL_MATRIX ( 1 .. N , 1 .. N ) ;
  begin
    for I in A'RANGE ( 1 ) loop
      for J in A'RANGE ( 2 ) loop
        A ( I , J ) := RATIONAL_ZERO ;
      end loop ;
      A ( I , I ) := RATIONAL_ONE ;
    end loop ;
    return A ;
  end IDENTITY_MATRIX ;

  function "*" ( F : RATIONAL ;
                 A : RATIONAL_MATRIX ) return RATIONAL_MATRIX is
    C : RATIONAL_MATRIX ( A'FIRST( 1 ) .. A'LAST ( 1 ) , A'FIRST ( 2 ) .. A'
      LAST ( 2 )) ;
  begin
    for I in A'RANGE ( 1 ) loop
      for J in A'RANGE ( 2 ) loop
        C ( I , J ) := A ( I , J ) * F ;
      end loop ;
    end loop ;
    return C ;
  end ;

  function "*" ( A : RATIONAL_MATRIX ;
                 F : RATIONAL ) return RATIONAL_MATRIX is
    C : RATIONAL_MATRIX ( A'FIRST( 1 ) .. A'LAST ( 1 ) , A'FIRST ( 2 ) .. A'
      LAST ( 2 )) ;
  begin
    for I in A'RANGE ( 1 ) loop
      for J in A'RANGE ( 2 ) loop
        C ( I , J ) := A ( I , J ) * F ;
      end loop ;
    end loop ;
    return C ;
  end ;

  function "*" ( F : RATIONAL ;
                 A : RATIONAL_VECTOR ) return RATIONAL_VECTOR is
    C : RATIONAL_VECTOR ( A'FIRST .. A'LAST ) ;
  begin
    for I in A'RANGE loop
      C ( I ) := A ( I ) * F ;
    end loop ;
    return C ;
  end "*" ;

  function "*" ( A : RATIONAL_VECTOR ;
                 F : RATIONAL ) return RATIONAL_VECTOR is
    C : RATIONAL_VECTOR ( A'FIRST .. A'LAST ) ;
  begin
    for I in A'RANGE loop
      C ( I ) := A ( I ) * F ;
    end loop ;
    return C ;
  end "*" ;

  function COLUMN_TO_VECTOR ( COLUMN : INTEGER ;
                              A : RATIONAL_MATRIX ) return RATIONAL_VECTOR is
    V : RATIONAL_VECTOR ( A'FIRST( 1 ) .. A'LAST ( 1 )) ;
  begin
    for J in A'RANGE ( 1 ) loop
      V ( J ) := A ( J , COLUMN ) ;
    end loop ;
    return V ;
  end COLUMN_TO_VECTOR ;

  procedure VECTOR_TO_COLUMN ( V : RATIONAL_VECTOR ;
                               COLUMN : INTEGER ;
                               A : in out RATIONAL_MATRIX ) is
  begin
    if V'LENGTH /= A'LENGTH ( 1 ) or
       COLUMN < A'FIRST ( 2 ) or
       COLUMN > A'LAST ( 2 ) then
      raise MATRIX_ERROR ;
    end if ;
    for I in A'RANGE ( 1 ) loop
      A ( I , COLUMN ) := V ( I - A'FIRST( 1 ) + V'FIRST) ;
    end loop ;
  end VECTOR_TO_COLUMN ;

  function ROW_TO_VECTOR ( ROW : INTEGER ;
                           A : RATIONAL_MATRIX ) return RATIONAL_VECTOR is
    V : RATIONAL_VECTOR ( A'FIRST( 2 ) .. A'LAST ( 2 )) ;
  begin
    if ROW < A'FIRST ( 1 ) or
       ROW > A'LAST ( 1 ) then
      raise MATRIX_ERROR ;
    end if ;
    for J in A'RANGE ( 2 ) loop
      V ( J ) := A ( ROW , J ) ;
    end loop ;
    return V ;
  end ROW_TO_VECTOR ;

  procedure VECTOR_TO_ROW ( V : RATIONAL_VECTOR ;
                            ROW : INTEGER ;
                            A : in out RATIONAL_MATRIX ) is
  begin
    if V'LENGTH /= A'LENGTH ( 2 ) or
       ROW < A'FIRST ( 1 ) or
       ROW > A'LAST ( 1 ) then
      raise MATRIX_ERROR ;
    end if ;
    for J in A'RANGE ( 2 ) loop
      A ( ROW , J ) := V ( J - A'FIRST( 2 ) + V'FIRST) ;
    end loop ;
  end VECTOR_TO_ROW ;

  function EXTRACT_SUB_MATRIX ( ROW_1 , ROW_2 , COL_1 , COL_2 : INTEGER ;
                                A : RATIONAL_MATRIX ) return RATIONAL_MATRIX
     is
    SUB_MATRIX : RATIONAL_MATRIX ( ROW_1 .. ROW_2 , COL_1 .. COL_2 ) ;
  begin
    if ROW_2 < ROW_1 or
       COL_2 < COL_1 then
      raise MATRIX_ERROR ;
    end if ;
    if ROW_1 < A'FIRST ( 1 ) or
       ROW_2 > A'LAST ( 1 ) or
       COL_1 < A'FIRST ( 2 ) or
       COL_2 > A'LAST ( 2 ) then
      raise MATRIX_ERROR ;
    end if ;
    for I in ROW_1 .. ROW_2 loop
      for J in COL_1 .. COL_2 loop
        SUB_MATRIX ( I , J ) := A ( I - ROW_1 + A'FIRST( 1 ) , J - COL_1 + A'
          FIRST ( 2 )) ;
      end loop ;
    end loop ;
    return SUB_MATRIX ;
  end EXTRACT_SUB_MATRIX ;

  procedure INSERT_SUB_MATRIX ( SUB_MATRIX : RATIONAL_MATRIX ;
                                ROW , COLUMN : INTEGER ;
                                A : in out RATIONAL_MATRIX ) is
  begin
    if ROW < A'FIRST ( 1 ) or
       COLUMN < A'FIRST ( 2 ) or
       ROW + SUB_MATRIX'LENGTH ( 1 ) - 1 > A'LAST ( 1 ) or
       COLUMN + SUB_MATRIX'LENGTH ( 2 ) - 1 > A'LAST ( 2 ) then
      raise MATRIX_ERROR ;
    end if ;
    for I in SUB_MATRIX'RANGE ( 1 ) loop
      for J in SUB_MATRIX'RANGE ( 2 ) loop
        A ( I - SUB_MATRIX'FIRST( 1 ) + ROW , J - SUB_MATRIX'FIRST ( 2 ) +
           COLUMN) := SUB_MATRIX ( I , J ) ;
      end loop ;
    end loop ;
  end INSERT_SUB_MATRIX ;

  procedure PUT ( A : RATIONAL_MATRIX ) is

  begin
    for I in A'RANGE ( 1 ) loop
      for J in A'RANGE ( 2 ) loop
        TEXT_IO.PUT ( " MAT(" ) ;
        TEXT_IO.PUT ( INTEGER'IMAGE(I) ) ;
        TEXT_IO.PUT ( " ," ) ;
        TEXT_IO.PUT ( INTEGER'IMAGE(J) ) ;
        TEXT_IO.PUT ( " ) = " ) ;
        PUT ( A( I , J )) ;
        TEXT_IO.NEW_LINE ;
      end loop ;
    end loop ;
  end PUT ;

  procedure PUT ( V : RATIONAL_VECTOR ) is

  begin
    for I in V'RANGE loop
      TEXT_IO.PUT ( " VEC(" ) ;
      TEXT_IO.PUT ( INTEGER'IMAGE(I) ) ;
      TEXT_IO.PUT ( " ) = " ) ;
      PUT ( V( I )) ;
      TEXT_IO.NEW_LINE ;
    end loop ;
  end PUT ;


end RATIONAL_ARRAYS ;