Return to ts_fltnum.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 Multprec_Natural_Numbers;           use Multprec_Natural_Numbers;
with Multprec_Natural_Numbers_io;        use Multprec_Natural_Numbers_io;
with Multprec_Integer_Numbers;           use Multprec_Integer_Numbers;
with Multprec_Integer_Numbers_io;        use Multprec_Integer_Numbers_io;
with Standard_Floating_Numbers;          use Standard_Floating_Numbers;
with Standard_Floating_Numbers_io;       use Standard_Floating_Numbers_io;
with Standard_Mathematical_Functions;    use Standard_Mathematical_Functions;
with Multprec_Floating_Numbers;          use Multprec_Floating_Numbers;
with Multprec_Floating_Numbers_io;       use Multprec_Floating_Numbers_io;
with Standard_Random_Numbers;            use Standard_Random_Numbers;
with Multprec_Random_Numbers;            use Multprec_Random_Numbers;

procedure ts_fltnum is

  tol : constant double_float := 10.0**(-8);

  procedure Read ( f : in out Floating_Number; name : in string ) is

    n : natural;

  begin
    put("Give " & name & " : "); get(f);
    put("Current size is "); put(Size_Fraction(f),1);
    put(".  Give expansion factor : "); get(n);
    if n > 0
     then Expand(f,n);
    end if;
  end Read;

  procedure Formatted_Output ( f : in Floating_Number ) is

  -- DESCRIPTION :
  --   Reads the format parameters and writes the floating-point number
  --   accordingly.

    fore,aft,exp : natural;
    
  begin
    put("Give the number of places before the decimal point : "); get(fore);
    put("Give the number of places after the decimal point : ");  get(aft);
    put("Give the number of places of the exponent : ");          get(exp);
    put("-> formatted : "); put(f,fore,aft,exp); new_line;
  end Formatted_Output;

  procedure Test_io is

  -- DESCRIPTION :
  --   Reads and writes a floating-point number.

    f,abf : Floating_Number;
    ans : character;

  begin
    put_line("Testing input/output for multi-precision floating numbers.");
    loop
      put("Give a floating number : "); get(f);
      put("-> your floating : "); put(f); new_line;
      abf := AbsVal(f);
      put("-> its absolute value : "); put(abf); new_line;
      put("-> #decimal places in fraction : ");
      put(Decimal_Places_Fraction(f),1); new_line;
      put("-> #decimal places in exponent : ");
      put(Decimal_Places_Exponent(f),1); new_line;
      put("Do you want formatted output ? (y/n) "); get(ans);
      if ans = 'y'
       then Formatted_Output(f);
      end if;
      Clear(f); Clear(abf);
      put("Do you want more tests ? (y/n) "); get(ans);
      exit when (ans /= 'y');
    end loop;
  end Test_io;

  function Truncate ( f : in double_float ) return integer is

    i : integer := integer(f);

  begin
    if i >= 0
     then if double_float(i) > f + tol
           then i := i-1;
          end if;
     else if double_float(i) < f - tol
           then i := i+1;
          end if;
    end if;
    return i;
  end Truncate;

  procedure Test_Creation is

    f : Floating_Number;
    d,fd : double_float;
    i : integer;
    ans : character;

  begin
    put_line("Testing the creation of multi-precision floating numbers.");
    loop
      put("Give an integer : "); get(i);
      put("-> your integer : "); put(i,1); new_line;
      f := Create(i);
      put("-> as floating number : "); put(f); new_line;
      put("Give a standard float : "); get(d);
      put("-> your float : "); put(d); new_line;
      f := Create(d);
      put("-> as floating number : "); put(f); new_line;
      fd := Round(f);
      put("-> rounded as standard float : "); put(fd); new_line;
      if d = fd
       then put_line("Creation/Rounding test is successful.");
       else put_line("Difference up to working precision ?");
            put("d - Round(Create(d)) : "); put(f-fd); new_line;
      end if;
      put("Give a floating number : "); get(f);
      put("-> your floating number : "); put(f); new_line;
      d := Round(f);
      put("-> rounded as float     :"); put(d); new_line;
      put("Do you want more tests ? (y/n) "); get(ans);
      exit when (ans /= 'y');
    end loop;
  end Test_Creation;

  procedure Test_Compare ( f1 : in Floating_Number; f2 : in double_float ) is
  begin
    if Equal(f1,f2)
     then put_line("The numbers are equal.");
     else put_line("The numbers are different.");
    end if;
    if f1 < f2
     then put_line("First number is less than second number.");
     else put_line("First number is not less than second number.");
    end if;
    if f1 > f2
     then put_line("First number is greater than second number.");
     else put_line("First number is not greater than second number.");
    end if;
  end Test_Compare;

  procedure Test_Compare ( f1,f2 : in Floating_Number ) is
  begin
    if Equal(f1,f2)
     then put_line("The numbers are equal.");
     else put_line("The numbers are different.");
    end if;
    if f1 < f2
     then put_line("First number is less than second number.");
     else put_line("First number is not less than second number.");
    end if;
    if f1 > f2
     then put_line("First number is greater than second number.");
     else put_line("First number is not greater than second number.");
    end if;
  end Test_Compare;

  procedure Zero_Test ( f : in Floating_Number ) is
  begin
    if Equal(f,0.0)
     then put_line(" equals zero.");
     else put_line(" is different from zero.");
    end if;
  end Zero_Test;

  procedure Test_Comparison is

    f1,f2 : Floating_Number;
   -- f2 : double_float;
    ans : character;

  begin
    put_line("Testing comparison/copying for multi-precision floats.");
    loop
      put("Give 1st number f1 : "); get(f1);
      put(" f1 : "); put(f1); 
      Zero_Test(f1);
      put("Give 2nd number f2 : "); get(f2);
      put(" f2 : "); put(f2);
      Zero_Test(f2);
      Test_Compare(f1,f2);
     -- Copy(f1,f2);
     -- put_line("After copy :");
     -- Test_Compare(f1,f2);
      put("Do you want more tests ? (y/n) "); get(ans);
      exit when ans /= 'y';
    end loop;
  end Test_Comparison;

  procedure Test_Size is

    f,mf : Floating_Number;
    ans : character;
    factor : integer;
    rnd : boolean;

  begin
    put_line("Testing trunc/round/expand for multi-precision floats");
    loop
      put("Give a floating number : "); get(f);
      put("-> your floating : "); put(f); new_line;
      put("The size of the fraction : "); put(Size_Fraction(f),1); new_line;
      loop
        put("Give size modificator : "); get(factor);
        if factor <= 0
         then put("Do you want to truncate or to round ? (t/r) "); get(ans);
              rnd := (ans = 'r');
        end if;
        if factor > 0
         then -- mf := Expand(f,factor);
              Expand(f,factor);
              put("expanded : "); put(f); -- put(mf);
              new_line;
         elsif factor < 0
             then if rnd
                   then -- mf := Round(f,-factor);
                        Round(f,-factor); put("rounded : ");
                   else -- mf := Trunc(f,-factor);
                        Trunc(f,-factor); put("truncated : ");
                  end if;
                  put(f); -- put(mf);
                  new_line;
             else if rnd
                   then -- mf := Round(f,factor);
                        Round(f,factor); put("rounded : ");
                   else -- mf := Trunc(f,factor);
                        Trunc(f,factor); put("truncated : ");
                  end if;
                  put(f); -- put(mf);
                  new_line;
                  -- mf := Expand(f,factor);
                  Expand(f,factor);
                  put("expanded : "); put(f); -- put(mf);
                  new_line;
        end if;
        put("Do you want other size modificators ? (y/n) "); get(ans);
        exit when (ans /= 'y');
      end loop;
      put("Do you want more tests ? (y/n) "); get(ans);
      exit when (ans /= 'y');
    end loop;
  end Test_Size;

  procedure Test_Addition is

    ans : character;
    f1,f2,sum1,sum2 : Floating_Number;

  begin
    put_line("Testing the addition operations.");
    loop
      Read(f1,"f1");
     -- put("Give 1st number f1 : "); get(f1);
      put("-> f1 : "); put(f1); new_line;
      Read(f2,"f2");
     -- put("Give 2nd number f2 : "); get(f2);
      put("-> f2 : "); put(f2); new_line;
      sum1 := f1+f2;
      put("f1+f2 : "); put(sum1); new_line;
      sum2 := f2+f1;
      put("f2+f1 : "); put(sum2); new_line;
      if Equal(sum1,sum2)
       then put_line("Test on commutativity is successful.");
       else put_line("Failure, bug detected.");
      end if;
      put("Do you want more tests ? (y/n) "); get(ans);
      exit when ans /= 'y';
    end loop;
  end Test_Addition;

  procedure Test_Subtraction is

    ans : character;
    f1,f2,diff : Floating_Number;

  begin
    put_line("Testing the subtraction operations.");
    loop
      Read(f1,"f1");
     -- put("Give 1st number f1 : "); get(f1);
      put("-> f1 : "); put(f1); new_line;
      Read(f2,"f2");
     -- put("Give 2nd number f2 : "); get(f2);
      put("-> f2 : "); put(f2); new_line;
      diff := f1-f2;
      put("f1 - f2 : "); put(diff); new_line;
      Add(diff,f2);
      put("(f1-f2)+f2 : "); put(diff); new_line;
      if Equal(diff,f1)
       then put_line("Test of subtraction is successful.");
       else put_line("Failure, bug detected.");
      end if;
      put("Do you want more tests ? (y/n) "); get(ans);
      exit when ans /= 'y';
    end loop;
  end Test_Subtraction;

  procedure Test_Multiplication is

    ans : character;
    f1,f2,prod1,prod2 : Floating_Number;

  begin
    put_line("Testing the multiplication operations.");
    loop
      Read(f1,"f1");
     -- put("Give 1st number : "); get(f1);
      put("-> f1 : "); put(f1); new_line;
      Read(f2,"f2");
     -- put("Give 2nd number : "); get(f2);
      put("-> f2 : "); put(f2); new_line;
      prod1 := f1*f2;
      put("Product f1*f2 : "); put(prod1); new_line;
      prod2 := f2*f1;
      put("Product f2*f1 : "); put(prod2); new_line;
      if Equal(prod1,prod2)
       then put_line("Test on commutativity is successful.");
       else put_line("Failure, product not commutative: bug!");
      end if;
      put("Do you want more tests ? (y/n) "); get(ans);
      exit when ans /= 'y';
    end loop;
  end Test_Multiplication;

  procedure Test_Exponentiation is

    ans : character;
    e1,e2 : Integer_Number;
    f,exp1,exp2,prod,expo : Floating_Number;

  begin
    put_line("Testing the exponentiation operations.");
    loop
      Read(f,"f");
     -- put("Give a number : "); get(f);
      put("-> your number f : "); put(f); new_line;
      put("Give 1st exponent : "); get(e1);
      put("-> your 1st exponent e1 : "); put(e1); new_line;
      exp1 := f**e1;
      put("f**e1 : "); put(exp1); new_line;
      put("Give 2nd exponent : "); get(e2);
      put("-> your 2nd exponent e2 : "); put(e2); new_line;
      exp2 := f**e2;
      put("f**e2 : "); put(exp2); new_line;
      prod := exp1*exp2;
      put("(f**e1)*(f**e2) : "); put(prod); new_line;
      expo := f**(e1+e2);
      put("f**(e1+e2)      : "); put(expo); new_line;
      if Equal(prod,expo)
       then put_line("Test of exponentiation is successful.");
       else put_line("Failure, bug detected.");
      end if;
      put("Do you want more tests ? (y/n) "); get(ans);
      exit when ans /= 'y';
    end loop;
  end Test_Exponentiation;

  procedure Test_Division is

    ans : character;
    f1,f2,quot,prod,diff : Floating_Number;

  begin
    put_line("Testing the division operations.");
    loop
      Read(f1,"f1");
     -- put("Give 1st number f1 : "); get(f1);
      put("-> f1 : "); put(f1); new_line;
      Read(f2,"f2");
     -- put("Give 2nd number f2 : "); get(f2);
      put("-> f2 : "); put(f2); new_line;
      prod := f1*f2;
      put("f1*f2 : "); put(prod); new_line;
      quot := prod/f2;
      put("(f1*f2)/f2 : "); put(quot); new_line;
      if Equal(quot,f1)
       then put_line("Test of division is successful.");
       else put("Failure, bug detected?");
            put_line("  Difference up to working precision?");
            diff := quot - f1;
            put("(f1*f2)/f2 - f1 : "); put(diff); new_line;
      end if;
      Copy(f1,quot);
      Div(quot,f2);    put("f1/f2 : "); put(quot); new_line;
      prod := quot*f2; put("(f1/f2)*f2 : "); put(prod); new_line;
                       put(" f1        : "); put(f1); new_line;
      if Equal(prod,f1)
       then put_line("Test of division/remainder computation is successful.");
       else put("Failure, bug detected?");
            put_line("  Difference up to working precision?");
            if prod > f1
             then diff := prod - f1;
             else diff := f1 - prod;
            end if;
            put("(f1/f2)*f2 - f1 : "); put(diff); new_line;
      end if;
      put("Do you want more tests ? (y/n) "); get(ans);
      exit when ans /= 'y';
    end loop;
  end Test_Division;

  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;

  procedure 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,tmp : Floating_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 : "); tmp := sum2-n1; put(tmp); new_line;
          Clear(tmp);
    end if;
    Add(sum2,n2);
    if Equal(sum2,sum1)
     then put("  Add okay");
     else put("  Add Bug?"); Report_Bug;
          put("diff : "); tmp := sum2-sum1; put(tmp); new_line;
          Clear(tmp);
    end if;
    Sub(sum2,n1);
    if Equal(sum2,n2)
     then put("  Sub okay"); new_line;
     else put("  Sub Bug?"); Report_Bug;
          put("diff : "); tmp := sum2-n2; put(tmp); new_line;
          Clear(tmp);
    end if;
    Clear(n1); Clear(n2);
    Clear(sum1); Clear(sum2);
  exception
    when others => put_line("input caused exception:"); Report_Bug; raise;
  end Random_Addition_and_Subtraction;

  procedure 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
      Random_Addition_and_Subtraction(sz1,sz2,low,upp);
    end loop;
  end Additions_and_Subtractions_on_Randoms;

  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,tmp : Floating_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 : "); tmp := quot-n1; put(tmp); new_line;
          Clear(tmp);
    end if;
    Mul(quot,n2);
    if Equal(prod,quot)
     then put("  Mul okay");
     else put("  Mul Bug?"); Report_Bug;
          put("Diff : "); tmp := quot-prod; put(tmp); new_line;
          Clear(tmp);
    end if;
    Div(prod,n1);
    if Equal(prod,n2)
     then put("  Div okay"); new_line;
     else put("  Div Bug?"); Report_Bug;
          put("Diff : "); tmp := prod-n2; put(tmp); new_line;
          Clear(tmp);
    end if;
    Clear(n1); Clear(n2);
    Clear(prod); Clear(quot);
  exception
    when others => 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 multi-precision floating numbers.");
    loop
      new_line;
      put_line("Choose one of the following : ");
      put_line("  0. exit program      1. Input/Output     2. Creation      ");
      put_line("  3. Comparison/Copy   4. Addition         5. Subtraction   ");
      put_line("  6. Multiplication    7. Exponentiation   8. Division      ");
      put_line("  9. Truncate/Round/Expand                                  ");
      put_line("  A. Addition/subtraction on randomly generated numbers.    ");
      put_line("  B. Multiplication/division/remainder on random numbers.   ");
      put("Type in your choice (0,1,2,3,4,5,6,7,8,9,A, or B) : "); get(ans);
      exit when (ans = '0');
      new_line;
      case ans is
        when '1' => Test_io;
        when '2' => Test_Creation;
        when '3' => Test_Comparison;
        when '4' => Test_Addition;
        when '5' => Test_Subtraction;
        when '6' => Test_Multiplication;
        when '7' => Test_Exponentiation;
        when '8' => Test_Division;
        when '9' => Test_Size;
        when 'A' => Additions_and_Subtractions_on_Randoms;
        when 'B' => Multiplications_and_Divisions_on_Randoms;
        when others => null;
      end case;
    end loop;
  end Main;

begin
  Main;
end ts_fltnum;