Return to ts_cmpnum.adb CVS log

Up to [local] / OpenXM_contrib / PHC / Ada / Math_Lib / Numbers

Import the second public release of PHCpack.

OKed by Jan Verschelde.

with text_io,integer_io;                 use text_io,integer_io;
with Standard_Complex_Numbers;           use Standard_Complex_Numbers;
with Standard_Complex_Numbers_io;        use Standard_Complex_Numbers_io;
with Standard_Complex_Numbers_Polar;     use Standard_Complex_Numbers_Polar;
with Multprec_Floating_Numbers;          use Multprec_Floating_Numbers;
with Multprec_Floating_Numbers_io;       use Multprec_Floating_Numbers_io;
with Multprec_Complex_Numbers;           use Multprec_Complex_Numbers;
with Multprec_Complex_Numbers_io;        use Multprec_Complex_Numbers_io;
with Standard_Random_Numbers;            use Standard_Random_Numbers;
with Multprec_Random_Numbers;            use Multprec_Random_Numbers;

procedure ts_cmpnum is

-- DESCRIPTION :
--   Interactive/Random testing on standard/multi-precision complex arithmetic.

  procedure Test_Standard_io is

    c : Standard_Complex_Numbers.Complex_Number;
    use Standard_Complex_Numbers;

  begin
    new_line;
    put_line("Testing input/output for standard complex numbers.");
    new_line;
    put("Give a complex number c : "); get(c);
    put("-> c : "); put(c); new_line;
    put("-> 1/c : "); put(1.0/c); new_line;
    put("-> 1/c : "); put(Create(1.0)/c); new_line;
  end Test_Standard_io;

  procedure Test_Multprec_io is

    c : Multprec_Complex_Numbers.Complex_Number;

  begin
    new_line;
    put_line("Testing input/output for multi-precision complex numbers.");
    new_line;
    put("Give a complex number c : "); get(c);
    put("-> c : "); put(c); new_line;
  end Test_Multprec_io;

  procedure Test_Roots is

    d,k : natural;
    a,c,prod : Standard_Complex_Numbers.Complex_Number;
    ans : character;

  begin
    new_line;
    put_line("Solving x^d - c = 0, with c a standard complex number.");
    new_line;
    put("Give the degree d : "); get(d);
    put("Give the constant c : "); get(c);
    loop
      put("Which root do you want ? "); get(k);
      a := Root(c,d,k);
      put("The root is "); put(a); new_line;
      prod := a;
      for j in 2..d loop
        prod := prod*a;
      end loop;
      put("root^d  =   "); put(prod); new_line;
      if Equal(prod,c)
       then put_line("root^d = c, test is successful.");
       else put_line("root^d /= c, bug detected? ");
            put("Difference : "); put(prod-c); new_line;
      end if;
      put("Do you want other roots ? (y/n) "); get(ans);
      exit when ans /= 'y';
    end loop;
  end Test_Roots;

  function Random ( sz : natural; low,upp : integer ) return Floating_Number is

  -- DESCRIPTION :
  --   Generates a random number of the given size, with exponent between
  --   the bounds low and upp.

    res : Floating_Number := Random(sz);
    exp : integer := Random(low,upp);

  begin
    if exp > 0
     then for i in 1..exp loop
            Mul(res,10.0);
          end loop;
     elsif exp < 0
         then for i in 1..(-exp) loop
                Div(res,10.0);
              end loop;
    end if;
    return res;
  end Random;

  function Random ( sz : natural; low,upp : integer )
                  return Multprec_Complex_Numbers.Complex_Number is

  -- DESCRIPTION :
  --   Generates a random number of the given size, with exponents for real
  --   and imaginary parts between low and upp.

  begin
    return Create(Random(sz,low,upp),Random(sz,low,upp));
  end Random;

  procedure Standard_Random_Addition_and_Subtraction is

  -- DESCRIPTION :
  --   Three tests are performed:
  --   1) n1+n2-n2 = n1, with "+" and "-".
  --   2) Add(n1,n2) is the same as n1 := n1+n2?
  --   3) Sub(n1+n2,n1) leads to n2?

    n1,n2,sum1,sum2 : Standard_Complex_Numbers.Complex_Number;

    procedure Report_Bug is
    begin
      new_line;
      put("  n1 : "); put(n1); new_line;
      put("  n2 : "); put(n2); new_line;
    end Report_Bug;

  begin
    n1 := Random;
    n2 := Random;
    sum1 := n1+n2;
    sum2 := sum1-n2;
    if Equal(sum2,n1)
     then put("n1+n2-n2 okay");
     else put("n1+n2-n2 Bug?"); Report_Bug;
          put("diff : "); put(sum2-n1); new_line;
    end if;
    Add(sum2,n2);
    if Equal(sum2,sum1)
     then put("  Add okay");
     else put("  Add Bug?"); Report_Bug;
          put("diff : "); put(sum2-sum1); new_line;
    end if;
    Sub(sum2,n1);
    if Equal(sum2,n2)
     then put("  Sub okay"); new_line;
     else put("  Sub Bug?"); Report_Bug;
          put("diff : "); put(sum2-n2); new_line;
    end if;
  exception
    when CONSTRAINT_ERROR => put_line("input caused exception:");
                             Report_Bug; raise;
  end Standard_Random_Addition_and_Subtraction;

  procedure Standard_Additions_and_Subtractions_on_Randoms is

  -- DESCRIPTION :
  --   Generates a number of random floats and performs repeated
  --   additions and subtractions with checks on consistencies.

    nb : natural;

  begin
    put("Give the number of tests : "); get(nb);
    for i in 1..nb loop
      Standard_Random_Addition_and_Subtraction;
    end loop;
  end Standard_Additions_and_Subtractions_on_Randoms;

  procedure Multprec_Random_Addition_and_Subtraction
              ( sz1,sz2 : in natural; low,upp : in integer ) is

  -- DESCRIPTION :
  --   Three tests are performed:
  --   1) n1+n2-n2 = n1, with "+" and "-".
  --   2) Add(n1,n2) is the same as n1 := n1+n2?
  --   3) Sub(n1+n2,n1) leads to n2?

    n1,n2,sum1,sum2 : Multprec_Complex_Numbers.Complex_Number;

    procedure Report_Bug is
    begin
      new_line;
      put("  n1 : "); put(n1); new_line;
      put("  n2 : "); put(n2); new_line;
    end Report_Bug;

  begin
    n1 := Random(sz1,low,upp);
    n2 := Random(sz2,low,upp);
    sum1 := n1+n2;
    sum2 := sum1-n2;
    if Equal(sum2,n1)
     then put("n1+n2-n2 okay");
     else put("n1+n2-n2 Bug?"); Report_Bug;
          put("diff : "); put(sum2-n1); new_line;
    end if;
    Add(sum2,n2);
    if Equal(sum2,sum1)
     then put("  Add okay");
     else put("  Add Bug?"); Report_Bug;
          put("diff : "); put(sum2-sum1); new_line;
    end if;
    Sub(sum2,n1);
    if Equal(sum2,n2)
     then put("  Sub okay"); new_line;
     else put("  Sub Bug?"); Report_Bug;
          put("diff : "); put(sum2-n2); new_line;
    end if;
    Clear(n1); Clear(n2);
    Clear(sum1); Clear(sum2);
  exception
    when CONSTRAINT_ERROR => put_line("input caused exception:");
                             Report_Bug; raise;
  end Multprec_Random_Addition_and_Subtraction;

  procedure Multprec_Additions_and_Subtractions_on_Randoms is

  -- DESCRIPTION :
  --   Generates a number of random floats and performs repeated
  --   additions and subtractions with checks on consistencies.

    nb,sz1,sz2 : natural;
    low,upp : integer;

  begin
    put("Give the number of tests : "); get(nb);
    put("Give the size of the 1st number : "); get(sz1);
    put("Give the size of the 2nd number : "); get(sz2);
    put("Give lower bound on exponent : "); get(low);
    put("Give upper bound on exponent : "); get(upp);
    for i in 1..nb loop
      Multprec_Random_Addition_and_Subtraction(sz1,sz2,low,upp);
    end loop;
  end Multprec_Additions_and_Subtractions_on_Randoms;

  procedure Interactive_Multiplication_and_Division is

    n1,n2,prod,quot : Multprec_Complex_Numbers.Complex_Number;
    ans : character;

  begin
    loop
      put("Give 1st number : "); get(n1);
      put("-> n1 : "); put(n1); new_line;
      put("Give 2nd number : "); get(n2);
      put("-> n2 : "); put(n2); new_line;
      prod := n1*n2;
      put("n1*n2 : "); put(prod); new_line;
      quot := prod/n2;
      put("(n1*n2)/n2  : "); put(quot); new_line;
      Clear(n1); Clear(n2); Clear(prod); Clear(quot);
      put("Do you want more tests ? (y/n) "); get(ans);
      exit when (ans /= 'y');
    end loop;
  end Interactive_Multiplication_and_Division;

  procedure Random_Multiplication_and_Division
               ( sz1,sz2 : in natural; low,upp : in integer ) is

  -- DESCRIPTION :
  --   Three tests are performed :
  --   1) n1*n2/n2 = n1, with "*" and "/".
  --   2) Mul(n1,n2) is the same as n1 := n1*n2 ?
  --   3) Div(n1*n2,n1) leads to n2 ?

    n1,n2,prod,quot : Multprec_Complex_Numbers.Complex_Number;

    procedure Report_Bug is
    begin
      new_line;
      put("  n1 : "); put(n1); new_line;
      put("  n2 : "); put(n2); new_line;
    end Report_Bug;

  begin
    n1 := Random(sz1,low,upp);
    n2 := Random(sz2,low,upp);
    prod := n1*n2;
    quot := prod/n2;
    if Equal(quot,n1)
     then put("n1*n2/n2 okay");
     else put("n1*n2/n2 Bug?"); Report_Bug;
          put("Diff : "); put(quot-n1); new_line;
    end if;
    Mul(quot,n2);
    if Equal(prod,quot)
     then put("  Mul okay");
     else put("  Mul Bug?"); Report_Bug;
          put("Diff : "); put(quot-prod); new_line;
    end if;
    Div(prod,n1);
    if Equal(prod,n2)
     then put("  Div okay"); new_line;
     else put("  Div Bug?"); Report_Bug;
          put("Diff : "); put(prod-n2); new_line;
    end if;
    Clear(n1); Clear(n2);
    Clear(prod); Clear(quot);
  exception
    when CONSTRAINT_ERROR => put_line("input caused exception :");
                             Report_Bug; raise;
  end Random_Multiplication_and_Division;

  procedure Multiplications_and_Divisions_on_Randoms is

  -- DESCRIPTION :
  --   Generates a number of random floats and performs repeated
  --   multiplications and divisions with checks on consistencies.

    nb,sz1,sz2 : natural;
    low,upp : integer;

  begin
    put("Give the number of tests : "); get(nb);
    put("Give the size of the 1st number : "); get(sz1);
    put("Give the size of the 2nd number : "); get(sz2);
    put("Give lower bound on exponent : "); get(low);
    put("Give upper bound on exponent : "); get(upp);
    for i in 1..nb loop
      Random_Multiplication_and_Division(sz1,sz2,low,upp);
    end loop;
  end Multiplications_and_Divisions_on_Randoms;

  procedure Main is

    ans : character;

  begin
    new_line;
    put_line("Interactive testing of standard and multi-precision "
                 & "complex numbers.");
    loop
      new_line;
      put_line("Choose one of the following : ");
      put_line("  0. Exit this program.                                     ");
      put_line("  1. Input/Output of standard complex numbers.              ");
      put_line("  2. Addition/subtraction on random standard numbers.       ");
      put_line("  3. Compute roots of unity of standard complex numbers.    ");
      put_line("  4. Input/Output of multi-precision complex numbers.       ");
      put_line("  5. Addition/subtraction on random multi-precision numbers.");
      put_line("  6. Multiplication/division/remainder on random "
                                              & "multi-precision numbers.   ");
      put_line("  7. Multiplication/division on user-given numbers.         ");
      put("Type in your choice (0,1,2,3,4,5,6, or 7) : "); get(ans);
      exit when (ans = '0');
      new_line;
      case ans is
        when '1' => Test_Standard_io;
        when '2' => Standard_Additions_and_Subtractions_on_Randoms;
        when '3' => Test_Roots;
        when '4' => Test_Multprec_io;
        when '5' => Multprec_Additions_and_Subtractions_on_Randoms;
        when '6' => Multiplications_and_Divisions_on_Randoms;
        when '7' => Interactive_Multiplication_and_Division;
        when others => null;
      end case;
    end loop;
  end Main;

begin
  Main;
end ts_cmpnum;