-- LONG_COMPLEX_ROOT_F.ADA    change the function F calculation to your function
--                            change REGION_MAX and REGION_MIN to bound area
--                            change ACCURACY to bound answer
--                            change MAX_ROOTS value to maximum nunber of roots

with LONG_COMPLEX_TYPES; use LONG_COMPLEX_TYPES;
package LONG_COMPLEX_ROOT_F is

  REGION_MIN : constant COMPLEX := (-10.0, -0.001) ;
  REGION_MAX : constant COMPLEX := ( 10.0,  0.001) ;
  ACCURACY   : constant COMPLEX := ( 0.000001, 0.0000000001) ;
  MAX_ROOTS  : constant INTEGER      := 20 ;

  function F( Z : COMPLEX ) return COMPLEX ;

end LONG_COMPLEX_ROOT_F;


with TEXT_IO; use TEXT_IO;
with LONG_COMPLEX_ELEMENTARY_FUNCTIONS; use LONG_COMPLEX_ELEMENTARY_FUNCTIONS;
package body LONG_COMPLEX_ROOT_F is
  function F( Z : COMPLEX ) return COMPLEX is
  begin
    return (1.0/Z);
  end F;
begin
  PUT_LINE("(1.0/Z)");
end LONG_COMPLEX_ROOT_F;

-- LONG_COMPLEX_ROOTS.ADA

with LONG_COMPLEX_ROOT_F; use LONG_COMPLEX_ROOT_F; 
                                            -- change this package for new F

with LONG_COMPLEX_TYPES; use LONG_COMPLEX_TYPES;
with LONG_COMPLEX_ARRAYS; use LONG_COMPLEX_ARRAYS;
with LONG_REAL_ARRAYS; use LONG_REAL_ARRAYS;
with MATHEMATICAL_CONSTANTS;
with TEXT_IO; use TEXT_IO;
with COMPLEX_IO;
with REAL_ARRAYS_IO;
with COMPLEX_ARRAYS_IO;

procedure LONG_COMPLEX_ROOTS is

  package CX_IO is new COMPLEX_IO(LONG_FLOAT, COMPLEX);
  use CX_IO;
  package CX_ARR_IO is new COMPLEX_ARRAYS_IO
                          (LONG_FLOAT, COMPLEX, COMPLEX_VECTOR, COMPLEX_MATRIX);
  use CX_ARR_IO;

  PI : constant := MATHEMATICAL_CONSTANTS.PI;

  Z_INTEGRAL : COMPLEX;
  LOCAL_MAX  : COMPLEX;
  LOCAL_MIN  : COMPLEX;
  TEST       : Boolean := False;

  function INTEGRATE(LOCAL_MIN, LOCAL_MAX : COMPLEX) return COMPLEX is
    A   : COMPLEX;
    B   : COMPLEX;
    X   : COMPLEX;
    SUM : COMPLEX := (0.0,0.0);
    INT : COMPLEX := (0.0,0.0);
    N : constant := 9;
    Z : REAL_VECTOR(0..N) := (
        -0.97390_65285_17172,
        -0.86506_33666_88985,
        -0.67940_95682_99024,
        -0.43339_53941_29247,
        -0.14887_43389_81631,
        +0.14887_43389_81631,
        +0.43339_53941_29247,
        +0.67940_95682_99024,
        +0.86506_33666_88985,
        +0.97390_65285_17172);
    W : REAL_VECTOR(0..N) := (
         0.06667_13443_08688,
         0.14945_13491_50581,
         0.21908_63625_15982,
         0.26926_67193_09996,
         0.29552_42247_14753,
         0.29552_42247_14753,
         0.26926_67193_09996,
         0.21908_63625_15982,
         0.14945_13491_50581,
         0.06667_13443_08688);
  begin
    A := (RE(LOCAL_MIN),IM(LOCAL_MIN));
    B := (RE(LOCAL_MAX),IM(LOCAL_MIN));
    for I in 0..N loop
      X := ((B-A)*Z(I)+(A+B))/2.0;
      SUM := SUM + W(I)/F(X);  -- poles of inverse are roots of function
    end loop;
    INT := SUM*(B-A)/2.0;
    if TEST then
      PUT("Integral from "); PUT(A); PUT(" to "); PUT(B); PUT(" is ");
      PUT(INT); NEW_LINE;
    end if; 
    SUM := (0.0,0.0);
    A := B;
    B := (RE(LOCAL_MAX),IM(LOCAL_MAX));
    for I in 0..N loop
      X := ((B-A)*Z(I)+(A+B))/2.0;
      SUM := SUM + W(I)/F(X);  -- poles of inverse are roots of function
    end loop;
    if TEST then
      PUT("Integral from "); PUT(A); PUT(" to "); PUT(B); PUT(" is ");
      PUT(SUM*(B-A)/2.0); NEW_LINE;
    end if; 
    INT := INT + SUM*(B-A)/2.0;
    SUM := (0.0,0.0);
    A := B;
    B := (RE(LOCAL_MIN),IM(LOCAL_MAX));
    for I in 0..N loop
      X := ((B-A)*Z(I)+(A+B))/2.0;
      SUM := SUM + W(I)/F(X);  -- poles of inverse are roots of function
    end loop;
    if TEST then
      PUT("Integral from "); PUT(A); PUT(" to "); PUT(B); PUT(" is ");
      PUT(SUM*(B-A)/2.0); NEW_LINE;
    end if; 
    INT := INT + SUM*(B-A)/2.0;
    SUM := (0.0,0.0);
    A := B;
    B := (RE(LOCAL_MIN),IM(LOCAL_MIN));
    for I in 0..N loop
      X := ((B-A)*Z(I)+(A+B))/2.0;
      SUM := SUM + W(I)/F(X);  -- poles of inverse are roots of function
    end loop;
    if TEST then
      PUT("Integral from "); PUT(A); PUT(" to "); PUT(B); PUT(" is ");
      PUT(SUM*(B-A)/2.0); NEW_LINE;
    end if; 
    INT := INT + SUM*(B-A)/2.0;
    return INT/(2.0*PI);
  exception
    when others =>
      PUT(X); PUT_LINE(" = Hit root");
      return (0.0,0.0); -- special case
  end INTEGRATE;

  function INTEGRATE2 return COMPLEX is
    INT : COMPLEX := (0.0,0.0);
    A, B : COMPLEX;
  begin
    A := (0.3,0.0);
    B := (0.2,0.2);
    INT := INT + ((1.0/F(A))+(1.0/F(B)))*((B-A)/2.0);
    A := B;
    B := (0.0,0.3);
    INT := INT + ((1.0/F(A))+(1.0/F(B)))*((B-A)/2.0);
    A := B;
    B := (-0.2,0.2);
    INT := INT + ((1.0/F(A))+(1.0/F(B)))*((B-A)/2.0);
    A := B;
    B := (-0.3,0.0);
    INT := INT + ((1.0/F(A))+(1.0/F(B)))*((B-A)/2.0);
    A := B;
    B := (-0.2,-0.2);
    INT := INT + ((1.0/F(A))+(1.0/F(B)))*((B-A)/2.0);
    A := B;
    B := (0.0,-0.3);
    INT := INT + ((1.0/F(A))+(1.0/F(B)))*((B-A)/2.0);
    A := B;
    B := (0.2,-0.2);
    INT := INT + ((1.0/F(A))+(1.0/F(B)))*((B-A)/2.0);
    A := B;
    B := (0.3,0.0);
    INT := INT + ((1.0/F(A))+(1.0/F(B)))*((B-A)/2.0);
    return INT/(2.0*PI);
  end INTEGRATE2;

begin
  PUT_LINE("Long Complex Roots");
  Z_INTEGRAL := INTEGRATE2;
  PUT(Z_INTEGRAL); PUT_LINE(" = trap rule about 0,0");
  NEW_LINE;

  TEST := True;
--  PUT(COMPLEX'(REGION_MIN)); PUT_LINE(" =min");
--  PUT(COMPLEX'(REGION_MAX)); PUT_LINE(" =max");
--  Z_INTEGRAL := INTEGRATE( REGION_MIN, REGION_MAX);
--  PUT(Z_INTEGRAL); PUT_LINE(" = integral");
--  NEW_LINE;

  PUT(COMPLEX'(-0.01,-0.01)); PUT_LINE(" =min");
  PUT(COMPLEX'(10.0,10.0)); PUT_LINE(" =max");
  Z_INTEGRAL := INTEGRATE( (-0.01,-0.01), (10.0,10.0));
  PUT(Z_INTEGRAL); PUT_LINE(" = integral");
  NEW_LINE;

--  PUT(COMPLEX'(2.0,2.0)); PUT_LINE(" =min");
--  PUT(COMPLEX'(-2.0,-2.0)); PUT_LINE(" =max");
--  Z_INTEGRAL := INTEGRATE( (2.0,2.0), (-2.0,-2.0));
--  PUT(Z_INTEGRAL); PUT_LINE(" = integral");
--  NEW_LINE;

--  PUT(COMPLEX'(-2.0,-2.0)); PUT_LINE(" =min");
--  PUT(COMPLEX'(2.0,2.0)); PUT_LINE(" =max");
--  Z_INTEGRAL := INTEGRATE( (-2.0,-2.0), (2.0,2.0));
--  PUT(Z_INTEGRAL); PUT_LINE(" = integral");
--  NEW_LINE;

--  PUT(COMPLEX'(-4.0,-4.0)); PUT_LINE(" =min");
--  PUT(COMPLEX'(4.0,4.0)); PUT_LINE(" =max");
--  Z_INTEGRAL := INTEGRATE( (-4.0,-4.0), (4.0,4.0));
--  PUT(Z_INTEGRAL); PUT_LINE(" = integral");
--  NEW_LINE;

  PUT(COMPLEX'(-1.0,-1.0)); PUT_LINE(" =min");
  PUT(COMPLEX'(1.0,1.0)); PUT_LINE(" =max");
  Z_INTEGRAL := INTEGRATE( (-1.0,-1.0), (1.0,1.0));
  PUT(Z_INTEGRAL); PUT_LINE(" = integral");
  NEW_LINE;

--  PUT(COMPLEX'(1.1,0.1)); PUT_LINE(" =min");
--  PUT(COMPLEX'(0.9,-0.1)); PUT_LINE(" =max");
--  Z_INTEGRAL := INTEGRATE( (1.1,0.1), (0.9,-0.1));
--  PUT(Z_INTEGRAL); PUT_LINE(" = integral");
--  NEW_LINE;

--  PUT(COMPLEX'(-1.5,-1.0)); PUT_LINE(" =min");
--  PUT(COMPLEX'(-0.5,1.0)); PUT_LINE(" =max");
--  Z_INTEGRAL := INTEGRATE( (-1.5,-1.0), (-0.5,1.0));
--  PUT(Z_INTEGRAL); PUT_LINE(" = integral");
--  NEW_LINE;

--  PUT(COMPLEX'(-1.0,-1.0)); PUT_LINE(" =min");
--  PUT(COMPLEX'(-0.5,1.0)); PUT_LINE(" =max");
--  Z_INTEGRAL := INTEGRATE( (-1.0,-1.0), (-0.5,1.0));
--  PUT(Z_INTEGRAL); PUT_LINE(" = integral");
--  NEW_LINE;

--  PUT(LONG_COMPLEX'(1.5,0.0)); PUT_LINE(" =min");
--  PUT(LONG_COMPLEX'(0.5,1.0)); PUT_LINE(" =max");
--  Z_INTEGRAL := INTEGRATE( (1.5,0.0), (0.5,1.0));
--  PUT(Z_INTEGRAL); PUT_LINE(" = integral");
--  NEW_LINE;

  PUT(COMPLEX'(-0.5,-0.5)); PUT_LINE(" =min");
  PUT(COMPLEX'(0.5,0.5)); PUT_LINE(" =max");
  Z_INTEGRAL := INTEGRATE( (-0.5,-0.5), (0.5,0.5));
  PUT(Z_INTEGRAL); PUT_LINE(" = integral");
  NEW_LINE;

--  PUT(COMPLEX'(-0.5,0.5)); PUT_LINE(" =min");
--  PUT(COMPLEX'(0.5,1.5)); PUT_LINE(" =max");
--  Z_INTEGRAL := INTEGRATE( (-0.5,0.5), (0.5,1.5));
--  PUT(Z_INTEGRAL); PUT_LINE(" = integral");
--  NEW_LINE;

--  PUT(COMPLEX'(-6.0,-2.0)); PUT_LINE(" =min");
--  PUT(COMPLEX'(3.0,1.5)); PUT_LINE(" =max");
--  Z_INTEGRAL := INTEGRATE( (-6.0,-2.0), (3.0,1.5));
--  PUT(Z_INTEGRAL); PUT_LINE(" = integral");
end LONG_COMPLEX_ROOTS;