-- Package body GENERIC_COMPLEX_TYPES
--
--Copyright (C) 1993 Free Software Foundation
--written by Jon Squire with assistance from many others
--
--This file is part of the numerics library.  This library
--is free software; you can redistribute it and/or modify it under the
--terms of the GNU Library General Public License as published by the Free
--Software Foundation; either version 2 of the License, or (at your
--option) any later version.  This library is distributed in the hope
--that it will be useful, but WITHOUT ANY WARRANTY; without even the
--implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
--PURPOSE.  See the GNU Library General Public License for more details.
--You should have received a copy of the GNU Library General Public
--License along with this library; if not, write to the Free Software
--Foundation, 675 Mass Ave, Cambridge, MA 02139, USA.
--

with GENERIC_ELEMENTARY_FUNCTIONS ; -- prior standard

package body GENERIC_COMPLEX_TYPES is

  package REAL_ELEMENTARY_FUNCTIONS is new
                                GENERIC_ELEMENTARY_FUNCTIONS ( REAL ) ;
  use REAL_ELEMENTARY_FUNCTIONS ;


-- Subprograms for type COMPLEX

  -- Selection, conversion and composition operations

  function RE ( X : COMPLEX ) return REAL is
  begin
    return X.RE ;
  end RE ;

  function IM ( X : COMPLEX ) return REAL is
  begin
    return X.IM ;
  end IM ;

  function IM ( X : IMAGINARY ) return REAL is
  begin
    return REAL(X) ;
  end IM ;

  procedure SET_RE ( X : in out COMPLEX ;
                     RE : in REAL ) is
  begin
    X.RE := RE ;
  end SET_RE ;

  procedure SET_IM ( X : in out COMPLEX ;
                     IM : in REAL ) is
  begin
    X.IM := IM ;
  end SET_IM ;

  procedure SET_IM ( X :    out IMAGINARY ;
                     IM : in REAL ) is
  begin
    X := IMAGINARY(IM) ;
  end SET_IM ;

  function COMPOSE_FROM_CARTESIAN ( RE : REAL ) return COMPLEX is
  begin
    return ( RE , 0.0 ) ;
  end COMPOSE_FROM_CARTESIAN ;

  function COMPOSE_FROM_CARTESIAN ( RE, IM : REAL ) return COMPLEX is
  begin
    return ( RE , IM ) ;
  end COMPOSE_FROM_CARTESIAN ;

  function COMPOSE_FROM_CARTESIAN ( IM : IMAGINARY ) return COMPLEX is
  begin
    return ( 0.0 , REAL(IM) ) ;
  end COMPOSE_FROM_CARTESIAN ;

  function COMPOSE_FROM_POLAR ( MODULUS, ARGUMENT : REAL )
                                                            return COMPLEX is
  begin
    return ( MODULUS * COS (ARGUMENT) , MODULUS * SIN (ARGUMENT) ) ;
  end COMPOSE_FROM_POLAR ;

  function COMPOSE_FROM_POLAR ( MODULUS, ARGUMENT, CYCLE : REAL )
                                                           return COMPLEX is
  begin
    if CYCLE <= 0.0 then
      raise ARGUMENT_ERROR ;
    end if ;
    return ( MODULUS * COS (ARGUMENT, CYCLE) ,
             MODULUS * SIN (ARGUMENT, CYCLE) ) ;
  end COMPOSE_FROM_POLAR ;

  function MODULUS ( X : COMPLEX ) return REAL is
  begin
    if X.RE = 0.0 and then X.IM = 0.0 then
      return 0.0 ;
    end if;
    if abs X.RE > abs X.IM then
      return abs X.RE * SQRT ( 1.0 + (X.IM/X.RE)*(X.IM/X.RE) );
    else
      return abs X.IM * SQRT ( 1.0 + (X.RE/X.IM)*(X.RE/X.IM) );
    end if;
  end MODULUS ;

  function ARGUMENT ( X : COMPLEX ) return REAL is
  begin
    if X.RE = 0.0 and then X.IM = 0.0 then
      return 0.0 ; -- by definition, a zero length vector returns ARGUMENT 0.0
    end if ;
    return ARCTAN ( X.IM , X.RE ) ;
  end ARGUMENT ;

  function ARGUMENT ( X : COMPLEX ;
                      CYCLE : REAL ) return REAL is
  begin
    if CYCLE <= 0.0 then
      raise ARGUMENT_ERROR ;
    end if ;
    if X.RE = 0.0 and then X.IM = 0.0 then
      return 0.0 ;
    end if ;
    return ARCTAN ( X.IM , X.RE , CYCLE ) ;
  end ARGUMENT ;

  -- pure IMAGINARY arithmetic operations

  function "<" (LEFT, RIGHT : IMAGINARY) return BOOLEAN is
  begin
    return REAL(LEFT) < REAL(RIGHT);
  end "<";

  function "<=" (LEFT, RIGHT : IMAGINARY) return BOOLEAN is
  begin
    return REAL(LEFT) <= REAL(RIGHT);
  end "<=";

  function ">" (LEFT, RIGHT : IMAGINARY) return BOOLEAN is
  begin
    return REAL(LEFT) > REAL(RIGHT);
  end ">";

  function ">=" (LEFT, RIGHT : IMAGINARY) return BOOLEAN is
  begin
    return REAL(LEFT) >= REAL(RIGHT);
  end ">=";

  function "+" (RIGHT : IMAGINARY) return IMAGINARY is
  begin
    return RIGHT;
  end "+";


  function "-" (RIGHT : IMAGINARY) return IMAGINARY is
  begin
    return IMAGINARY(-REAL(RIGHT));
  end "-";


  function "+" (LEFT, RIGHT : IMAGINARY) return IMAGINARY is
  begin
    return IMAGINARY(REAL(LEFT) + REAL(RIGHT));
  end "+";


  function "-" (LEFT, RIGHT : IMAGINARY) return IMAGINARY is
  begin
    return IMAGINARY(REAL(LEFT) - REAL(RIGHT));
  end "-";


  -- COMPLEX arithmetic operations

  function "+" ( RIGHT : COMPLEX ) return COMPLEX is
  begin
    return RIGHT ;
  end "+";


  function "-" ( RIGHT : COMPLEX ) return COMPLEX is
  begin
    return ( - RIGHT.RE , - RIGHT.IM ) ;
  end "-";


  function CONJUGATE ( X : COMPLEX ) return COMPLEX is
  begin
    return ( X.RE , - X.IM ) ;
  end CONJUGATE ;


  function "+" ( LEFT , RIGHT : COMPLEX ) return COMPLEX is
  begin
    return ( LEFT.RE + RIGHT.RE , LEFT.IM + RIGHT.IM ) ;
  end "+";


  function "-" ( LEFT , RIGHT : COMPLEX ) return COMPLEX is
  begin
    return ( LEFT.RE - RIGHT.RE , LEFT.IM - RIGHT.IM ) ;
  end "-";


  function "*" ( LEFT , RIGHT : COMPLEX ) return COMPLEX is
  begin
    if LEFT.RE = 0.0 and then LEFT.IM = 0.0 then
      return (0.0,0.0) ;
    elsif RIGHT.RE = 0.0 and then RIGHT.IM = 0.0 then
      return (0.0,0.0) ;
    end if ;
    if  LEFT.RE * RIGHT.RE - LEFT.IM * RIGHT.IM  /= 0.0 or else
        LEFT.RE * RIGHT.IM + LEFT.IM * RIGHT.RE  /= 0.0 then
      return ( LEFT.RE * RIGHT.RE - LEFT.IM * RIGHT.IM , LEFT.RE * RIGHT.IM + LEFT.IM * RIGHT.RE ) ;
    end if;
    return (((2.0*LEFT.RE) * RIGHT.RE - (2.0*LEFT.IM) * RIGHT.IM) /
              REAL'(2.0) ,
            ((2.0*LEFT.RE) * RIGHT.IM + (2.0*LEFT.IM) * RIGHT.RE) /
              REAL'(2.0) ) ;
  exception
    when CONSTRAINT_ERROR | NUMERIC_ERROR =>
      return (((LEFT.RE/2.0) * RIGHT.RE - (LEFT.IM/2.0) * RIGHT.IM) *
                REAL'(2.0) ,
              ((LEFT.RE/2.0) * RIGHT.IM + (LEFT.IM/2.0) * RIGHT.RE) *
                REAL'(2.0) ) ;
  end "*";


  function "/" ( LEFT , RIGHT : COMPLEX ) return COMPLEX is
    R : REAL ;
  begin
    if abs RIGHT.RE > abs RIGHT.IM then
      R := 1.0 + (RIGHT.IM/RIGHT.RE)*(RIGHT.IM/RIGHT.RE);
      return ( (LEFT.RE/RIGHT.RE + (LEFT.IM/RIGHT.RE)*(RIGHT.IM/RIGHT.RE))/R ,
               (LEFT.IM/RIGHT.RE - (LEFT.RE/RIGHT.RE)*(RIGHT.IM/RIGHT.RE))/R ) ;
    else
      R := (RIGHT.RE/RIGHT.IM)*(RIGHT.RE/RIGHT.IM) + 1.0;
      return ( ((LEFT.RE/RIGHT.IM)*(RIGHT.RE/RIGHT.IM) + LEFT.IM/RIGHT.IM)/R ,
               ((LEFT.IM/RIGHT.IM)*(RIGHT.RE/RIGHT.IM) - LEFT.RE/RIGHT.IM)/R ) ;
    end if;
  exception
    when CONSTRAINT_ERROR | NUMERIC_ERROR =>
      if abs RIGHT.RE > abs RIGHT.IM then
        R := 1.0 + (RIGHT.IM/RIGHT.RE)*(RIGHT.IM/RIGHT.RE);
        return ( (LEFT.RE/R)/RIGHT.RE + ((LEFT.IM/R)/RIGHT.RE)*(RIGHT.IM/RIGHT.RE) ,
                 (LEFT.IM/R)/RIGHT.RE - ((LEFT.RE/R)/RIGHT.RE)*(RIGHT.IM/RIGHT.RE) ) ;
      else
        R := (RIGHT.RE/RIGHT.IM)*(RIGHT.RE/RIGHT.IM) + 1.0;
        return ( ((LEFT.RE/R)/RIGHT.IM)*(RIGHT.RE/RIGHT.IM) + (LEFT.IM/R)/RIGHT.IM ,
                 ((LEFT.IM/R)/RIGHT.IM)*(RIGHT.RE/RIGHT.IM) - (LEFT.RE/R)/RIGHT.IM ) ;
      end if;
  end "/";


  -- REAL/IMAGINARY arithmetic operations --

  function "*" (LEFT , RIGHT : IMAGINARY) return REAL is
  begin
    return -REAL(LEFT)*REAL(RIGHT);
  end "*";


  function "*" (LEFT : IMAGINARY;
                RIGHT : REAL) return IMAGINARY is
  begin
    return IMAGINARY(REAL(LEFT)*RIGHT);
  end "*";


  function "*" (LEFT : REAL;
                RIGHT : IMAGINARY) return IMAGINARY is
  begin
    return IMAGINARY(LEFT*REAL(RIGHT));
  end "*";


  function "/" (LEFT , RIGHT : IMAGINARY) return REAL is
  begin
    return REAL(LEFT)/REAL(RIGHT);
  end "/";


  function "/" (LEFT : IMAGINARY;
                RIGHT : REAL) return IMAGINARY is
  begin
    return IMAGINARY(REAL(LEFT)/RIGHT);
  end "/";


  function "/" (LEFT : REAL;
                RIGHT : IMAGINARY) return IMAGINARY is
  begin
    return IMAGINARY(-LEFT/REAL(RIGHT));
  end "/";


  -- REAL/COMPLEX arithmetic operations --

  function "+" (LEFT : REAL;
                RIGHT : COMPLEX) return COMPLEX is
  begin
    return ( LEFT + RIGHT.RE , RIGHT.IM ) ;
  end "+" ;


  function "+" (LEFT : COMPLEX;
                RIGHT : REAL) return COMPLEX is
  begin
    return ( LEFT.RE + RIGHT , LEFT.IM ) ;
  end "+" ;


  function "-" (LEFT : REAL;
                RIGHT : COMPLEX) return COMPLEX is
  begin
    return ( LEFT - RIGHT.RE , - RIGHT.IM ) ;
  end "-" ;


  function "-" (LEFT : COMPLEX;
                RIGHT : REAL) return COMPLEX is
  begin
    return ( LEFT.RE - RIGHT , LEFT.IM ) ;
  end "-" ;


  function "*" (LEFT : REAL;
                RIGHT : COMPLEX) return COMPLEX is
  begin
    return ( LEFT * RIGHT.RE , LEFT * RIGHT.IM ) ;
  end "*" ;


  function "*" (LEFT : COMPLEX;
                RIGHT : REAL) return COMPLEX is
  begin
    return ( LEFT.RE * RIGHT , LEFT.IM * RIGHT ) ;
  end "*" ;


  function "/" (LEFT : REAL;
                RIGHT : COMPLEX) return COMPLEX is
    R : REAL ;
  begin
    if abs RIGHT.RE > abs RIGHT.IM then
      R := 1.0 + (RIGHT.IM/RIGHT.RE)*
                 (RIGHT.IM/RIGHT.RE);
      return ( (LEFT/RIGHT.RE)/R ,
             - ((LEFT/RIGHT.RE)*(RIGHT.IM/RIGHT.RE))/R ) ;
    else
      R := (RIGHT.RE/RIGHT.IM)*
           (RIGHT.RE/RIGHT.IM) + 1.0;
      return (
               ((LEFT/RIGHT.IM)*(RIGHT.RE/RIGHT.IM))/R ,
              - (LEFT/RIGHT.IM)/R ) ;
    end if;
  exception
    when CONSTRAINT_ERROR | NUMERIC_ERROR =>
      if abs RIGHT.RE > abs RIGHT.IM then
        R := 1.0 + (RIGHT.IM/RIGHT.RE)*
                   (RIGHT.IM/RIGHT.RE);
        return ( (LEFT/RIGHT.RE)/R ,
              - ((LEFT/R)/RIGHT.RE)*(RIGHT.IM/RIGHT.RE) ) ;
      else
        R := (RIGHT.RE/RIGHT.IM)*
             (RIGHT.RE/RIGHT.IM) + 1.0;
        return (
                 ((LEFT/R)/RIGHT.IM)*(RIGHT.RE/RIGHT.IM) ,
                - (LEFT/R)/RIGHT.IM ) ;
      end if;
  end "/";


  function "/" (LEFT : COMPLEX;
                RIGHT : REAL) return COMPLEX is
  begin
    return ( LEFT.RE / RIGHT , LEFT.IM / RIGHT ) ;
  end "/" ;


  -- IMAGINARY/COMPLEX arithmetic operations --

  function "+" (LEFT : IMAGINARY;
                RIGHT : COMPLEX) return COMPLEX is
  begin
    return (RIGHT.RE , REAL(LEFT) + RIGHT.IM ) ;
  end "+";


  function "+" (LEFT : COMPLEX;
                RIGHT : IMAGINARY) return COMPLEX is
  begin
    return (LEFT.RE , LEFT.IM + REAL(RIGHT));
  end "+";


  function "-" (LEFT : IMAGINARY;
                RIGHT : COMPLEX) return COMPLEX is
  begin
    return (-RIGHT.RE , REAL(LEFT) - RIGHT.IM ) ;
  end "-";


  function "-" (LEFT : COMPLEX;
                RIGHT : IMAGINARY) return COMPLEX is
  begin
    return (LEFT.RE , LEFT.IM - REAL(RIGHT));
  end "-";


  function "*" (LEFT : IMAGINARY;
                RIGHT : COMPLEX) return COMPLEX is
  begin
    return (-REAL(LEFT)*RIGHT.IM , REAL(LEFT)*RIGHT.RE);
  end "*";


  function "*" (LEFT : COMPLEX;
                RIGHT : IMAGINARY) return COMPLEX is
  begin
    return (-REAL(RIGHT)*LEFT.IM , REAL(RIGHT)*LEFT.RE);
  end "*";


  function "/" (LEFT : IMAGINARY;
                RIGHT : COMPLEX) return COMPLEX is
    C : COMPLEX := ( 0.0 , REAL(LEFT) ) ;
  begin
    return C / RIGHT;
  end "/";


  function "/" (LEFT : COMPLEX;
                RIGHT : IMAGINARY) return COMPLEX is
  begin
    return ( LEFT.IM/REAL(RIGHT) , -LEFT.RE/REAL(RIGHT) );
  end "/";


  -- Mixed REAL/IMAGINARY/COMPLEX

  function "+" (LEFT : REAL;
                RIGHT : IMAGINARY) return COMPLEX is
  begin
    return ( LEFT , REAL(RIGHT) ) ;
  end "+";


  function "+" (LEFT : IMAGINARY;
                RIGHT : REAL) return COMPLEX is
  begin
    return ( RIGHT , REAL(LEFT) ) ;
  end "+";


  function "-" (LEFT : REAL;
                RIGHT : IMAGINARY) return COMPLEX is
  begin
    return ( LEFT , -REAL(RIGHT) ) ;
  end "-";


  function "-" (LEFT : IMAGINARY;
                RIGHT : REAL) return COMPLEX is
  begin
    return ( -RIGHT , REAL(LEFT) ) ;
  end "-";


-- exponentiation for COMPLEX and IMAGINARY

  function "**" ( LEFT : COMPLEX ;
                  RIGHT : INTEGER ) return COMPLEX is
    VALUE : COMPLEX ;
  begin
    if RIGHT = 0 then -- by definition
      return (1.0,0.0) ;
    elsif REAL'(abs LEFT)  = 0.0 then -- by definition
      return (0.0,0.0) ;
    elsif LEFT = (1.0,0.0) then -- by definition
      return (1.0,0.0) ;
    end if ;
    VALUE := LEFT;
    if RIGHT > 0 then
      for I in 1 .. RIGHT-1 loop
        VALUE := VALUE * LEFT ;
      end loop ;
      return VALUE ;
    else
      begin
        for I in 1 .. abs(RIGHT+1) loop
          VALUE := VALUE * LEFT ;
        end loop ;
      exception -- overflow here (the denominator) means underflow
        when NUMERIC_ERROR | CONSTRAINT_ERROR =>
          return (0.0,0.0) ;
      end ;
      VALUE := (1.0,0.0) / VALUE ; -- overflow here takes care of itself
    end if ;
    return VALUE ;
  end "**" ;


  function "**" (LEFT : IMAGINARY;
                 RIGHT : INTEGER) return COMPLEX is
  begin
    if RIGHT = 0 then -- by definition
      return (1.0,0.0) ;
    elsif REAL(LEFT)  = 0.0 then -- by definition
      return (0.0,0.0) ;
    end if ;
    if RIGHT mod 4 = 1 then
      return (0.0, REAL(LEFT)**RIGHT) ;
    elsif RIGHT mod 4 = 2 then
      return (-REAL(LEFT)**RIGHT, 0.0) ;
    elsif RIGHT mod 4 = 3 then
      return (0.0, -REAL(LEFT)**RIGHT) ;
    else
      return (REAL(LEFT)**RIGHT, 0.0) ;
    end if;
  end "**" ;

  function "abs" (RIGHT : IMAGINARY) return IMAGINARY is
  begin
    return IMAGINARY(abs REAL(RIGHT)) ;
  end "abs";

  function "abs" (RIGHT : IMAGINARY) return REAL is
  begin
    return abs REAL(RIGHT) ;
  end "abs";

end GENERIC_COMPLEX_TYPES ;