Return to multprec_mathematical_functions.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;                            use text_io;
with Multprec_Floating_Numbers_io;       use Multprec_Floating_Numbers_io;

with Standard_Floating_Numbers;          use Standard_Floating_Numbers;
with Standard_Mathematical_Functions;    use Standard_Mathematical_Functions;

package body Multprec_Mathematical_Functions is

-- EXPONENTIAL AND LOGARITHMIC FUNCTIONS :

  function "**" ( x,y : Floating_Number ) return Floating_Number is

    res : Floating_Number;
    sx : double_float := Round(x);
    sy : double_float := Round(y);

  begin
    sx := sx**sy;
    res := Create(sx);
    return res;
  end "**";

  function LOG2 ( x : Floating_Number ) return Floating_Number is

    res : Floating_Number;
    sx : double_float := Round(x);

  begin
    sx := LOG2(sx);
    res := Create(sx);
    return res;
  end LOG2;

  function LOG10 ( x : Floating_Number ) return Floating_Number is

    res : Floating_Number;
    sx : double_float := Round(x);

  begin
    sx := LOG10(sx);
    res := Create(sx);
    return res;
  end LOG10;

  function SQRT ( x : Floating_Number ) return Floating_Number is

    res : Floating_Number;
    sx : double_float := Round(x);
    sizex : natural := Size_Fraction(x);
    sizeres : natural;

    procedure Iterate ( a : in Floating_Number ) is

      oneres,twores : Floating_Number;

    begin
      for i in 1..3 loop
        Copy(res,oneres);      put("init   : "); put(oneres); new_line;
        twores := 2.0*oneres;  put("2*init : "); put(twores); new_line;
        Mul(res,res);
        Sub(res,a);       put("residual : "); put(res); new_line;
        Div(res,twores);  put("inc : "); put(res); new_line;
        twores := oneres - res;
        Copy(twores,res);
        Clear(twores); Clear(oneres);
      end loop;
    end Iterate;

  begin
    sx := SQRT(sx);
    res := Create(sx); 
    sizeres := Size_Fraction(res);
    if (sx /= 0.0) and (sizeres < sizex)
     then Expand(res,sizex-sizeres);
          if sx >= 1.0
           then Iterate(x);
           else declare
                  invx : Floating_Number := 1.0/x;
                  invres : Floating_Number := 1.0/res;
                begin
                  Copy(invres,res);  Clear(invres);
                  Iterate(invx);
                  invres := 1.0/res;
                  Copy(invres,res);  Clear(invres);
                end;
          end if;
    end if;
    return res;
  end SQRT;

-- TRIGONOMETRIC FUNCTIONS :

  function SIN ( x : Floating_Number ) return Floating_Number is

    res : Floating_Number;
    sx : double_float := Round(x);

  begin
    sx := SIN(sx);
    res := Create(sx);
    return res;
  end SIN;

  function COS ( x : Floating_Number ) return Floating_Number is

    res : Floating_Number;
    sx : double_float := Round(x);

  begin
    sx := COS(sx);
    res := Create(sx);
    return res;
  end COS;

  function TAN ( x : Floating_Number ) return Floating_Number is

    res : Floating_Number;
    sx : double_float := Round(x);

  begin
    sx := TAN(sx);
    res := Create(sx);
    return res;
  end TAN;

  function ARCSIN ( x : Floating_Number ) return Floating_Number is

    res : Floating_Number;
    sx : double_float := Round(x);

  begin
    sx := ARCSIN(sx);
    res := Create(sx);
    return res;
  end ARCSIN;

  function ARCCOS ( x : Floating_Number ) return Floating_Number is

    res : Floating_Number;
    sx : double_float := Round(x);

  begin
    sx := ARCCOS(sx);
    res := Create(sx);
    return res;
  end ARCCOS;

  function ARCTAN ( x : Floating_Number ) return Floating_Number is

    res : Floating_Number;
    sx : double_float := Round(x);

  begin
    sx := ARCTAN(sx);
    res := Create(sx);
    return res;
  end ARCTAN;

  function Radius ( x,y : Floating_Number ) return Floating_Number is

    res : Floating_Number;
    sx : double_float := Round(x);
    sy : double_float := Round(y);

  begin
    sx := Radius(sx,sy);
    res := Create(sx);
    return res;
  end Radius;

  function Angle  ( x,y : Floating_Number ) return Floating_Number is

    res : Floating_Number;
    sx : double_float := Round(x);
    sy : double_float := Round(y);

  begin
    sx := Angle(sx,sy);
    res := Create(sx);
    return res;
  end Angle;

end Multprec_Mathematical_Functions;