-- test complex series against complex_elementary_functions
-- just a coarse test
-- ON_BRANCH_CUT_xxx used to avoid that issue here

-- This generic is attempted to be instantiated two times for the
-- types FLOAT and LONG_FLOAT.
-- Link and execute the instantiations that successfully compile:
-- TEST_COMPLEX_SERIES
-- TEST_LONG_COMPLEX_SERIES


with text_io; use text_io;
with generic_real_arrays;
with generic_complex_types;
with generic_complex_arrays;
with complex_io;
with generic_complex_elementary_functions;
with generic_complex_constraints;
with generic_elementary_functions;

generic
  type real is digits <>;
procedure test_generic_complex_series;

procedure test_generic_complex_series is

  package real_arrays is new generic_real_arrays(real);
  use real_arrays;

  package complex_types is new generic_complex_types(real);
  use complex_types;

  package complex_arrays is new generic_complex_arrays
                            (real, real_vector, real_matrix, complex);
  use complex_arrays;

  package cx_io is new complex_io(real, complex);
  use cx_io;

  package complex_constraints is new generic_complex_constraints
                                     ( real, complex );
  use complex_constraints;

  package complex_elementary_functions is new
          generic_complex_elementary_functions( real, complex, imaginary);
  use complex_elementary_functions;

  package elementary_functions is new generic_elementary_functions(real);
  use elementary_functions;

  package real_io is new text_io.float_io(real);
  use real_io;

  PI : constant := 3.14159_26535_89793_23846_26433_83279_50288_41971 ;
  PI_2 : constant := PI/2.0;
  SQRT_2 : constant := 1.41421_35623_73095_04880_16887_24209_69807_85696 ;

  complex_i : complex := compose_from_cartesian(0.0,1.0);
  z : complex;
  zz : complex;
  zr : complex;
  z_function : complex;
  z_series : complex;
  x : real;
  y : real;
  epsilon : real;
  tolerance : real;
  width : natural := real'digits+1;

  data : array(1..15) of real:=
        ( 0.0 , 1.0/32768.0, -1.0/32768.0, 1.0/1024.0, -1.0/1024.0,
          1.0/128.0, -1.0/128.0, 1.0/32.0, -1.0/32.0,
          0.125, -0.125, 0.25, -0.25, 0.5, -0.5);

  fact_1 : constant := 1.0;
  fact_2 : constant := 2.0 * fact_1;
  fact_3 : constant := 3.0 * fact_2;
  fact_4 : constant := 4.0 * fact_3;
  fact_5 : constant := 5.0 * fact_4;
  fact_6 : constant := 6.0 * fact_5;
  fact_7 : constant := 7.0 * fact_6;
  fact_8 : constant := 8.0 * fact_7;
  fact_9 : constant := 9.0 * fact_8;
  fact_10 : constant := 10.0 * fact_9;
  fact_11 : constant := 11.0 * fact_10;

  minimum_print : boolean := true; -- true if print only first time
                                   -- function differs from series

  subtype function_name_string_type is string(1..20);
  type function_name_record is
    record
      function_name : function_name_string_type;
      zp : complex; -- point argument
      zf : complex; -- GCEF function
      zs : complex; -- series
      max_error : real; -- vector error for reference
      excl : integer; -- values excluded, on branch cut
  end record;

  type    function_name_list_type is array(positive range <>) of
                    function_name_record;
  function_name : function_name_string_type;
  function_name_list : function_name_list_type(1..50);
  function_name_list_length : natural := 0;

  function printed(zp,zf,zs:complex;
                    name:string) return boolean is
    max_error : real;
  begin
    begin
      max_error := abs(zf-zs)/abs(zs);
    exception
      when constraint_error =>
        max_error := abs(zf-zs);
    end;
    function_name := name & (1..function_name'length-name'length=>' ');
    for i in 1..function_name_list_length loop
      if function_name = function_name_list(i).function_name then
        if function_name_list(i).max_error < max_error then
          function_name_list(i).zp := zp;
          function_name_list(i).zf := zf;
          function_name_list(i).zs := zs;
          function_name_list(i).max_error := max_error;
        end if;
        return true and minimum_print;
      end if;
    end loop;
    function_name_list_length := function_name_list_length+1;
    function_name_list(function_name_list_length).function_name :=function_name;
    function_name_list(function_name_list_length).zp := zp;
    function_name_list(function_name_list_length).zf := zf;
    function_name_list(function_name_list_length).zs := zs;
    function_name_list(function_name_list_length).max_error := max_error;
    return false;
  end printed;

  procedure print(zp,zf,zs:complex;
                  name:string) is
  begin
    if not printed(zp,zf,zs,name) then
      put(zf,aft=>width); put_line(" = gcef " & name);
      put(zs,aft=>width); put_line(" = series");
      put("    "); put(zp,aft=>width); put_line(" = z");
      new_line;
    end if;
  end print;

  procedure final_print is
  begin
    new_line;
    put_line("Final print");
    new_line;
    for i in 1..function_name_list_length loop
      put(function_name_list(i).zf,aft=>width);
      put_line(" = gcef " & function_name_list(i).function_name);
      put(function_name_list(i).zs,aft=>width); put_line(" = series");
      put("    "); put(function_name_list(i).zp,aft=>width); put_line(" = z");
      put("        "); put(function_name_list(i).max_error);
      put_line(" = max relative vector error");
      new_line;
    end loop;
  end final_print;

  procedure check(zp,zf,zs:complex;
                  name:string) is
  begin
    if RE(zs)=0.0 then
      if RE(zf)/=0.0 then
        print(zp,zf,zs,name);
      end if;
    else
      if abs((RE(zf)-RE(zs))/RE(zs))>tolerance then
        print(zp,zf,zs,name);
      end if;
    end if;
    if IM(zs)=0.0 then
      if IM(zf)/=0.0 then
        print(zp,zf,zs,name);
      end if;
    else
      if abs((IM(zf)-IM(zs))/IM(zs))>tolerance then
        print(zp,zf,zs,name);
      end if;
    end if;
  exception
    when others =>
      put(zp); put_line(" z caused exception in check of " & name);
  end check;

  function square(z:complex) return complex is
  begin
    x := (RE(z)-IM(z)) * (RE(z)+IM(z));
    y := RE(z)*IM(z);
    return compose_from_cartesian(x,y+y);
  end square;

begin

  new_line;
  put_line("Test Series against complex functions");

  epsilon := real'base'epsilon;
  while epsilon/real(real'machine_radix) + 1.0 /= 1.0 loop
    epsilon := epsilon/real(real'machine_radix);
  end loop;
  tolerance := 32.0*epsilon;
  put(tolerance);
  put_line(" = component relative error bound being used.");
  new_line;

  for i in data'range loop
    for j in data'range loop

      z := compose_from_cartesian(data(i),data(j));
      zz := square(z);
      zr := sqrt(z);

      if not on_branch_cut_sqrt(z) then
        z_function := SQRT(1.0+z);
        z_series := 1.0 + (1.0/2.0 - (1.0/8.0 - (1.0/16.0 - (5.0/128.0
                    - (7.0/256.0 - (21.0/1024.0 - (33.0/2048.0
                    - (429.0/32768.0 - (715.0/65536.0
                    - (2431.0/262144.0 - (4199.0/524288.0)
                    * z ) * z ) * z ) * z ) * z ) * z ) * z )
                    * z ) * z ) * z ) * z;
        check(1.0+z,z_function,z_series,"SQRT(1+z)");
      end if;

      if not on_branch_cut_log(z) then
        z_function := LOG(1.0+z);
        z_series := (1.0 - (1.0/2.0 - (1.0/3.0 - (1.0/4.0 - (1.0/5.0 - (1.0/6.0
                    - (1.0/7.0 - (1.0/8.0 - (1.0/9.0 - (1.0/10.0 - (1.0/11.0
                    - (1.0/12.0 - (1.0/13.0 - (1.0/14.0)
                    * z ) * z ) * z ) * z ) * z ) * z ) * z )
                    * z ) * z ) * z ) * z ) * z ) * z ) * z;
        check(1.0+z,z_function,z_series,"LOG(1+z)");
      end if;

      z_function := EXP(z);
      z_series := 1.0 + (1.0 + (1.0/fact_2 + (1.0/fact_3 + (1.0/fact_4
                  + (1.0/fact_5 + (1.0/fact_6 + (1.0/fact_7
                  + (1.0/fact_8 + (1.0/fact_9 + (1.0/fact_10
                  + (1.0/fact_11) * z ) * z ) * z ) * z ) * z ) * z )
                  * z ) * z ) * z ) * z ) * z;
      check(z,z_function,z_series,"EXP");

      z_function := SIN(z);
      z_series := (1.0 - (1.0/fact_3 - (1.0/fact_5 - (1.0/fact_7 - (1.0/fact_9
                  - (1.0/fact_11) * zz ) * zz ) * zz ) * zz ) * zz ) * z;
      check(z,z_function,z_series,"SIN");

      z_function := COS(z);
      z_series := 1.0 - (1.0/fact_2 - (1.0/fact_4 - (1.0/fact_6 - (1.0/fact_8
                  - (1.0/fact_10) * zz ) * zz ) * zz ) * zz ) * zz;
      check(z,z_function,z_series,"COS");

      z_function := TAN(z);
      z_series := (1.0 + (1.0/3.0 + (2.0/15.0 + (17.0/315.0
                  + (62.0/2835.0 + (1392.0/155925.0)
                  * zz ) * zz ) * zz ) * zz ) * zz ) * z;
      check(z,z_function,z_series,"TAN");

      if not on_branch_cut_cot(z) then
        z_function := COT(z);
        z_series := 1.0/z - (1.0/3.0 + (1.0/45.0 + (2.0/945.0
                    + (1.0/4725.0 + (2.0/93555.0 + (1382.0/638512875.0)
                    * zz ) * zz ) * zz ) * zz ) * zz ) * z;
        check(z,z_function,z_series,"COT");
      end if;

      if not on_branch_cut_arcsin(z) then
        z_function := ARCSIN(z);
        z_series := (1.0 + (1.0/6.0 + (3.0/40.0 + (5.0/112.0
                    + (35.0/1152.0 + (63.0/2816.0)
                    * zz ) * zz ) * zz ) * zz ) * zz ) * z;
        check(z,z_function,z_series,"ARCSIN");

        z_function := ARCSIN(z-1.0);
        z_series := -PI_2 + SQRT_2*(sqrt(z) + sqrt(square(z)*z)/12.0
                    + 3.0*sqrt(square(square(z))*z)/160.0
                    + 5.0*sqrt(square(square(z))*square(z)*z)/896.0
                    +35.0*sqrt(square(square(square(z)))*z)/18432.0
          +63.0*sqrt(square(square(square(z)))*square(z)*z)/90112.0
        +231.0*sqrt(square(square(square(z)))*square(square(z))*z)/851968.0);
        check(z,z_function,z_series,"ARCSIN(z-1)");

        z_function := ARCSIN(z+1.0);
        z_series := PI_2
                    + complex_i*SQRT_2*(-zr
                    + zr*z/12.0
                    - 3.0*zr*zz/160.0
                    + 5.0*zr*zz*z/896.0
                    - 35.0*zr*zz*zz/18432.0
                    + 63.0*zr*zz*zz*z/90112.0
                    - 231.0*zr*zz*zz*zz/851968.0);
        if (IM(z_series) < 0.0 and IM(z) > 0.0) or
           (IM(z_series) > 0.0 and IM(z) < 0.0) then
          set_IM(z_series,-IM(z_series));
        end if;
        check(z,z_function,z_series,"ARCSIN(z+1)");

        z_function := ARCSIN(z-(0.0,1.0));
        z_series := - complex_i*log(1.0+sqrt_2)
                    + sqrt_2*z/2.0
                    - complex_i*sqrt_2*square(z)/8.0
                    - sqrt_2*square(z)*z/48.0;
        check(z,z_function,z_series,"ARCSIN(z-i)");

        z_function := ARCSIN(z+(0.0,1.0));
        z_series := complex_i*log(1.0+sqrt_2)
                    + sqrt_2*z/2.0
                    + complex_i*sqrt_2*square(z)/8.0
                    - sqrt_2*square(z)*z/48.0;
        check(z,z_function,z_series,"ARCSIN(z+i)");
      end if;

      if not on_branch_cut_arccos(z) then
        z_function := ARCCOS(z);
        z_series := PI/2.0 - (1.0 + (1.0/6.0 + (3.0/40.0 + (5.0/112.0
                    + (35.0/1152.0 + (63.0/2816.0)
                    * zz ) * zz ) * zz ) * zz ) * zz ) * z;
        check(z,z_function,z_series,"ARCCOS");

        z_function := ARCCOS(z-1.0);
        z_series := PI - SQRT_2*(sqrt(z) + sqrt(square(z)*z)/12.0
                    + 3.0*sqrt(square(square(z))*z)/160.0
                    + 5.0*sqrt(square(square(z))*square(z)*z)/896.0
                    +35.0*sqrt(square(square(square(z)))*z)/18432.0
          +63.0*sqrt(square(square(square(z)))*square(z)*z)/90112.0
        +231.0*sqrt(square(square(square(z)))*square(square(z))*z)/851968.0);
        check(z,z_function,z_series,"ARCCOS(z-1)");

        z_function := ARCCOS(z+1.0);
        z_series := complex_i*SQRT_2*(zr
                    - zr*z/12.0
                    + 3.0*zr*zz/160.0
                    - 5.0*zr*zz*z/896.0
                    + 35.0*zr*zz*zz/18432.0
                    - 63.0*zr*zz*zz*z/90112.0
                    + 231.0*zr*zz*zz*zz/851968.0);
        check(z,z_function,z_series,"ARCCOS(z+1)");

        z_function := ARCCOS(z-(0.0,1.0));
        z_series := pi_2 + complex_i*log(1.0+SQRT_2)
                    - sqrt_2*z/2.0
                    + complex_i*sqrt_2*square(z)/8.0
                    + sqrt_2*square(z)*z/48.0;
        check(z,z_function,z_series,"ARCCOS(z-i)");

        z_function := ARCCOS(z+(0.0,1.0));
        z_series := pi_2 - complex_i*log(1.0+sqrt_2)
                    - sqrt_2*z/2.0
                    - complex_i*sqrt_2*square(z)/8.0
                    + sqrt_2*square(z)*z/48.0;
        check(z,z_function,z_series,"ARCCOS(z+i)");
      end if;

      if not on_branch_cut_arctan(z) then
        z_function := ARCTAN(z);
        z_series := (1.0 - (1.0/3.0 - (1.0/5.0 - (1.0/7.0 - (1.0/9.0
                    - (1.0/11.0 - (1.0/13.0 - (1.0/15.0)
                    * zz ) * zz ) * zz ) * zz ) * zz ) * zz ) * zz ) * z;
        check(z,z_function,z_series,"ARCTAN");
      end if;

      if not on_branch_cut_arccot(z) then
        z_function := ARCCOT(z);
        z_series := PI/2.0 - (1.0 - (1.0/3.0 - (1.0/5.0 - (1.0/7.0 - (1.0/9.0
                    - (1.0/11.0 - (1.0/13.0 - (1.0/15.0 - (1.0/17.0)
                    * zz ) * zz ) * zz ) * zz ) * zz ) * zz ) * zz ) * zz ) * z;
        check(z,z_function,z_series,"ARCCOT");
      end if;

      z_function := SINH(z);
      z_series := (1.0 + (1.0/fact_3 + (1.0/fact_5 + (1.0/fact_7 + (1.0/fact_9
                  + (1.0/fact_11)
                  * zz ) * zz ) * zz ) * zz ) * zz ) * z;
      check(z,z_function,z_series,"SINH");

      z_function := COSH(z);
      z_series := 1.0 + (1.0/fact_2 + (1.0/fact_4 + (1.0/fact_6 + (1.0/fact_8
                  + (1.0/fact_10)
                  * zz ) * zz ) * zz ) * zz ) * zz;
      check(z,z_function,z_series,"COSH");

      z_function := TANH(z);
      z_series := (1.0 - (1.0/3.0 - (2.0/15.0 - (17.0/315.0
                  - (62.0/2835.0 - (1382.0/155925.0)
                  * zz ) * zz ) * zz ) * zz ) * zz ) * z;
      check(z,z_function,z_series,"TANH");

      if not on_branch_cut_coth(z) then
        z_function := COTH(z);
        z_series := 1.0/z + (1.0/3.0 - (1.0/45.0 - (2.0/945.0
                    - (1.0/4725.0 - (2.0/93555.0 - (1382.0/638512875.0)
                    * zz ) * zz ) * zz ) * zz ) * zz ) * z;
        check(z,z_function,z_series,"COTH");
      end if;

      if not on_branch_cut_arcsinh(z) then
        z_function := ARCSINH(z);
        z_series := (1.0 - (1.0/6.0 - (3.0/40.0 - (5.0/112.0
                    - (35.0/1152.0 - (63.0/2816.0)
                    * zz ) * zz ) * zz ) * zz ) * zz ) * z;
       check(z,z_function,z_series,"ARCSINH");
      end if;

      if not on_branch_cut_arccosh(z) then
        z_function := ARCCOSH(z);
        z_series := complex_i*( (-pi/2.0) + (1.0+ (1.0/6.0 + (3.0/40.0
                    + (5.0/112.0 + (35.0/1152.0 + (63.0/2816.0)
                    * zz ) * zz ) * zz ) * zz ) * zz ) * z);
        if RE(z_series) <= 0.0 then
          z_series := - z_series;
        end if;
        check(z,z_function,z_series,"ARCCOSH");
      end if;

      if not on_branch_cut_arctanh(z) then
        z_function := ARCTANH(z);
        z_series := (1.0 + (1.0/3.0 + (1.0/5.0 + (1.0/7.0 + (1.0/9.0
                    + (1.0/11.0 + (1.0/13.0)
                    * zz ) * zz ) * zz ) * zz ) * zz ) * zz ) * z;
        check(z,z_function,z_series,"ARCTANH");
      end if;

      if not on_branch_cut_arccoth(z) then
        z_function := ARCCOTH(z);
        z_series := complex_i*PI/2.0 + (1.0 + (1.0/3.0 + (1.0/5.0 + (1.0/7.0
                    + (1.0/9.0 + (1.0/11.0 + (1.0/13.0)
                    * zz ) * zz ) * zz ) * zz ) * zz ) * zz ) * z;
        check(z,z_function,z_series,"ARCCOTH");
      end if;
    end loop;
  end loop;
  put_line("End First Pass of Test Series");
  final_print;
  put_line("End Test Series");
end test_generic_complex_series;

with test_generic_complex_series;
procedure test_complex_series is new
          test_generic_complex_series(float);

with test_generic_complex_series;
procedure test_long_complex_series is new
          test_generic_complex_series(long_float);