--
with TEXT_IO ; use TEXT_IO ;

package body GENERIC_DIGITS_ARRAYS is


-- DIGITS MATRIX ARITHMETIC PACKAGE BODY
--
--    Purpose :  To provide the basic set of D_NUMBER matrix arithmetic
--               capability in ada
--
--    Method :   This package defines the types
--                       DIGITS_MATRIX  and  DIGITS_VECTOR
--               this package also defines the overloaded procedure
--                put  for the types  DIGITS_MATRIX  and  DIGITS_VECTOR
--
--               this package then defines the overloaded operators
--
--                 +   -   *    unary +   unary -
--
--                 for combining types  DIGITS_MATRIX  DIGITS_VECTOR  and  D_NUMBER
--
--
--    WRITTEN BY : JON SQUIRE , 27 FEB 1983
--    MODIFIED 10 JUNE 1987 FOR DECIMAL FLOATING POINT
--
--    COPYRIGHT  1983  WESTINGHOUSE ELECTRIC CORP.
--    COPYRIGHT  1984  WESTINGHOUSE ELECTRIC CORP. -- additions
--    COPYRIGHT  1987  WESTINGHOUSE ELECTRIC CORP. -- modification
--
--    FOR THE TYPE  DIGITS_MATRIX  THE FIRST SUBSCRIPT IS THE ROW NUMBER,
--                               THE SECOND SUBSCRIPT IS THE COLUMN NUMBER
--

  CONSTANT_ZERO : D_NUMBER := DIGITS_ZERO ;
  CONSTANT_ONE : D_NUMBER := DIGITS_ONE ;

  function "+" ( A , B : DIGITS_MATRIX ) return DIGITS_MATRIX is

    C : DIGITS_MATRIX ( A'FIRST( 1 ) .. A'LAST ( 1 ) , A'FIRST ( 2 ) .. A'LAST (
        2 )) ;
  begin

-- CHECK THAT MATRICIES 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

-- MATRICIES HAVE DIFFERENT INDICIES
      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 : DIGITS_MATRIX ) return DIGITS_MATRIX is

    C : DIGITS_MATRIX ( A'FIRST( 1 ) .. A'LAST ( 1 ) , A'FIRST ( 2 ) .. A'LAST (
        2 )) ;
  begin

-- CHECK THAT MATRICIES 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

-- MATRICIES HAVE DIFFERENT INDICIES
      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 : DIGITS_MATRIX ) return DIGITS_MATRIX is

    C : DIGITS_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 ) := CONSTANT_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 : DIGITS_VECTOR ) return DIGITS_VECTOR is

    C : DIGITS_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 : DIGITS_VECTOR ) return DIGITS_VECTOR is

    C : DIGITS_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 : DIGITS_VECTOR ;
                 A : DIGITS_MATRIX ) return DIGITS_VECTOR is

    CV : DIGITS_VECTOR ( A'FIRST( 2 ) .. A'LAST ( 2 )) ;
    SUM : D_NUMBER ;
  begin
    if A'LENGTH ( 1 ) /= V'LENGTH then
      raise MATRIX_ERROR ;
    end if ;
    for J in A'RANGE ( 2 ) loop
      SUM := CONSTANT_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 : DIGITS_MATRIX ;
                 V : DIGITS_VECTOR ) return DIGITS_VECTOR is

    CV : DIGITS_VECTOR ( A'FIRST( 1 ) .. A'LAST ( 1 )) ;
    SUM : D_NUMBER ;
  begin
    if A'LENGTH ( 2 ) /= V'LENGTH then
      raise MATRIX_ERROR ;
    end if ;
    for I in A'RANGE ( 1 ) loop
      SUM := CONSTANT_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 : DIGITS_VECTOR ;
                         Q : DIGITS_VECTOR ) return D_NUMBER is

    F : D_NUMBER := CONSTANT_ZERO ;
    TEMP : D_NUMBER ;
  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 : DIGITS_VECTOR ;
                           Q : DIGITS_VECTOR ) return DIGITS_VECTOR is

    CV : DIGITS_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 : DIGITS_MATRIX ) return DIGITS_VECTOR is

    CV : DIGITS_VECTOR ( A'FIRST( 2 ) .. A'LAST ( 2 )) ;
    N : INTEGER := A'LENGTH ( 2 ) ;

-- RETURN N DIMENSIONAL VECTOR ORTHAGONAL 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 DIGITS_MATRIX is

    A : DIGITS_MATRIX ( 1 .. N , 1 .. N ) ;
  begin
    for I in A'RANGE ( 1 ) loop
      for J in A'RANGE ( 2 ) loop
        A ( I , J ) := CONSTANT_ZERO ;
      end loop ;
      A ( I , I ) := CONSTANT_ONE ;
    end loop ;
    return A ;
  end IDENTITY_MATRIX ;

  function "*" ( F : D_NUMBER ;
                 A : DIGITS_MATRIX ) return DIGITS_MATRIX is

    C : DIGITS_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 : DIGITS_MATRIX ;
                 F : D_NUMBER ) return DIGITS_MATRIX is

    C : DIGITS_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 : D_NUMBER ;
                 A : DIGITS_VECTOR ) return DIGITS_VECTOR is

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

  function "*" ( A : DIGITS_VECTOR ;
                 F : D_NUMBER ) return DIGITS_VECTOR is

    C : DIGITS_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 : DIGITS_MATRIX ) return DIGITS_VECTOR is

    V : DIGITS_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 : DIGITS_VECTOR ;
                               COLUMN : INTEGER ;
                               A : in out DIGITS_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 : DIGITS_MATRIX ) return DIGITS_VECTOR is

    V : DIGITS_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 : DIGITS_VECTOR ;
                            ROW : INTEGER ;
                            A : in out DIGITS_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 : DIGITS_MATRIX ) return DIGITS_MATRIX is

    SUB_MATRIX : DIGITS_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 : DIGITS_MATRIX ;
                                ROW , COLUMN : INTEGER ;
                                A : in out DIGITS_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 : DIGITS_MATRIX ) is

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

  procedure PUT ( V : DIGITS_VECTOR ) is

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

  procedure GET ( A : out DIGITS_MATRIX ) is

  begin
    for I in A'RANGE ( 1 ) loop
      for J in A'RANGE ( 2 ) loop
        GET ( A( I , J )) ;
      end loop ;
    end loop ;
  end GET ;

  procedure GET ( V : out DIGITS_VECTOR ) is

  begin
    for I in V'RANGE loop
      GET ( V( I )) ;
    end loop ;
  end GET ;
end GENERIC_DIGITS_ARRAYS ;