home *** CD-ROM | disk | FTP | other *** search
/ Geek Gadgets 1 / ADE-1.bin / ade-dist / gnat-2.06-src.tgz / tar.out / fsf / gnat / ada / eval_fat.adb < prev    next >
Text File  |  1996-09-28  |  12KB  |  500 lines
  1. ------------------------------------------------------------------------------
  2. --                                                                          --
  3. --                         GNAT COMPILER COMPONENTS                         --
  4. --                                                                          --
  5. --                             E V A L _ F A T                              --
  6. --                                                                          --
  7. --                                 B o d y                                  --
  8. --                                                                          --
  9. --                            $Revision: 1.2 $                              --
  10. --                                                                          --
  11. -- GNAT is free software;  you can  redistribute it  and/or modify it under --
  12. -- terms of the  GNU General Public License as published  by the Free Soft- --
  13. -- ware  Foundation;  either version 2,  or (at your option) any later ver- --
  14. -- sion.  GNAT is distributed in the hope that it will be useful, but WITH- --
  15. -- OUT ANY WARRANTY;  without even the  implied warranty of MERCHANTABILITY --
  16. -- or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License --
  17. -- for  more details.  You should have  received  a copy of the GNU General --
  18. -- Public License  distributed with GNAT;  see file COPYING.  If not, write --
  19. -- to the Free Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA. --
  20. --                                                                          --
  21. ------------------------------------------------------------------------------
  22.  
  23. with Stand;  use Stand;
  24. with Ttypef; use Ttypef;
  25.  
  26. package body Eval_Fat is
  27.  
  28.    Radix : Uint renames Uint_2;
  29.    --  This code is currently only correct for the radix 2 case
  30.  
  31.    function Float_Radix return Ureal renames Ureal_2;
  32.    --  Radix expressed in real form
  33.  
  34.    function Float_Radix_Inv return Ureal renames Ureal_Half;
  35.    --  Inverse of radix expressed in real form
  36.  
  37.    -----------------------
  38.    -- Local Subprograms --
  39.    -----------------------
  40.  
  41.    procedure Decompose
  42.      (RT       : R;
  43.       X        : in T;
  44.       Fraction : out T;
  45.       Exponent : out UI);
  46.    --  Decomposes a floating-point number into fraction and exponent parts
  47.  
  48.    function Machine_Mantissa (RT : R) return UI;
  49.    --  Get value of machine mantissa
  50.  
  51.    function Radix_To_M_Minus_1 (RT : R) return T;
  52.    --  Returns value of radix to power of mantissa minus one
  53.  
  54.    --------------
  55.    -- Adjacent --
  56.    --------------
  57.  
  58.    function Adjacent (RT : R; X, Towards : T) return T is
  59.    begin
  60.       if Towards = X then
  61.          return X;
  62.  
  63.       elsif Towards > X then
  64.          return Succ (RT, X);
  65.  
  66.       else
  67.          return Pred (RT, X);
  68.       end if;
  69.    end Adjacent;
  70.  
  71.    -------------
  72.    -- Ceiling --
  73.    -------------
  74.  
  75.    function Ceiling (RT : R; X : T) return T is
  76.       XT : constant T := Truncation (RT, X);
  77.  
  78.    begin
  79.       if UR_Is_Negative (X) then
  80.          return XT;
  81.  
  82.       elsif X = XT then
  83.          return X;
  84.  
  85.       else
  86.          return XT + Ureal_1;
  87.       end if;
  88.    end Ceiling;
  89.  
  90.    -------------
  91.    -- Compose --
  92.    -------------
  93.  
  94.    function Compose (RT : R; Fraction : T; Exponent : UI) return T is
  95.       Arg_Frac : T;
  96.       Arg_Exp  : UI;
  97.  
  98.    begin
  99.       Decompose (RT, Fraction, Arg_Frac, Arg_Exp);
  100.       return Scaling (RT, Arg_Frac, Exponent);
  101.    end Compose;
  102.  
  103.    ---------------
  104.    -- Copy_Sign --
  105.    ---------------
  106.  
  107.    function Copy_Sign (RT : R; Value, Sign : T) return T is
  108.       Result : T;
  109.  
  110.    begin
  111.       Result := abs Value;
  112.  
  113.       if UR_Is_Negative (Sign) then
  114.          return -Result;
  115.       else
  116.          return Result;
  117.       end if;
  118.    end Copy_Sign;
  119.  
  120.    ---------------
  121.    -- Decompose --
  122.    ---------------
  123.  
  124.    procedure Decompose
  125.      (RT       : R;
  126.       X        : in T;
  127.       Fraction : out T;
  128.       Exponent : out UI)
  129.    is
  130.       Factor   : T;
  131.       My_Fract : T;
  132.       Sign_X   : T;
  133.       Abs_X    : T;
  134.       Unit     : UI;
  135.       My_Exp   : UI;
  136.       Scale    : UI;
  137.  
  138.    begin
  139.       if X = Ureal_0 then
  140.          Fraction := Ureal_0;
  141.          Exponent := Uint_0;
  142.          return;
  143.       end if;
  144.  
  145.       Abs_X := abs X;
  146.  
  147.       if UR_Is_Negative (X) then
  148.          Sign_X := -Ureal_1;
  149.       else
  150.          Sign_X := Ureal_1;
  151.       end if;
  152.  
  153.       if Abs_X >= Ureal_1 then
  154.          Factor := Float_Radix_Inv;
  155.          Unit := Uint_1;
  156.  
  157.       else
  158.          Factor := Float_Radix;
  159.          Unit := -Uint_1;
  160.       end if;
  161.  
  162.       My_Fract := Abs_X;
  163.       My_Exp := Uint_0;
  164.  
  165.       while My_Fract >= Ureal_1
  166.         or else My_Fract < Ureal_Half
  167.       loop
  168.          My_Fract := My_Fract * Factor;
  169.          My_Exp := My_Exp + Unit;
  170.       end loop;
  171.  
  172.       --  The remaining step is to truncate the fraction to the appropriate
  173.       --  number of digits, since we have been doing this in full precision.
  174.  
  175.       Scale := Radix ** Machine_Mantissa (RT);
  176.       My_Fract := UR_From_Uint (UR_Trunc (My_Fract * Scale)) / Scale;
  177.  
  178.       Fraction := Sign_X * My_Fract;
  179.       Exponent := My_Exp;
  180.       return;
  181.  
  182.    end Decompose;
  183.  
  184.    --------------
  185.    -- Exponent --
  186.    --------------
  187.  
  188.    function Exponent (RT : R; X : T) return UI is
  189.       X_Frac : T;
  190.       X_Exp  : UI;
  191.  
  192.    begin
  193.       Decompose (RT, X, X_Frac, X_Exp);
  194.       return X_Exp;
  195.    end Exponent;
  196.  
  197.    -----------
  198.    -- Floor --
  199.    -----------
  200.  
  201.    function Floor (RT : R; X : T) return T is
  202.       XT : constant T := Truncation (RT, X);
  203.  
  204.    begin
  205.       if UR_Is_Positive (X) then
  206.          return XT;
  207.  
  208.       elsif XT = X then
  209.          return X;
  210.  
  211.       else
  212.          return XT - Ureal_1;
  213.       end if;
  214.    end Floor;
  215.  
  216.    --------------
  217.    -- Fraction --
  218.    --------------
  219.  
  220.    function Fraction (RT : R; X : T) return T is
  221.       X_Frac : T;
  222.       X_Exp  : UI;
  223.  
  224.    begin
  225.       Decompose (RT, X, X_Frac, X_Exp);
  226.       return X_Frac;
  227.    end Fraction;
  228.  
  229.    ------------------
  230.    -- Leading_Part --
  231.    ------------------
  232.  
  233.    function Leading_Part (RT : R; X : T; Radix_Digits : UI) return T is
  234.       L    : UI;
  235.       Y, Z : T;
  236.  
  237.    begin
  238.       if Radix_Digits >= Machine_Mantissa (RT) then
  239.          return X;
  240.  
  241.       else
  242.          L := Exponent (RT, X) - Radix_Digits;
  243.          Y := Truncation (RT, Scaling (RT, X, -L));
  244.          Z := Scaling (RT, Y, L);
  245.          return Z;
  246.       end if;
  247.  
  248.    end Leading_Part;
  249.  
  250.    -------------
  251.    -- Machine --
  252.    -------------
  253.  
  254.    function Machine (RT : R; X : T) return T is
  255.       X_Frac : T;
  256.       X_Exp  : UI;
  257.  
  258.    begin
  259.       Decompose (RT, X, X_Frac, X_Exp);
  260.       return Compose (RT, X_Frac, X_Exp);
  261.    end Machine;
  262.  
  263.    ----------------------
  264.    -- Machine_Mantissa --
  265.    ----------------------
  266.  
  267.    function Machine_Mantissa (RT : R) return UI is
  268.    begin
  269.       if RT = Standard_Short_Float then
  270.          return UI_From_Int (Short_Float_Attr_Machine_Mantissa);
  271.  
  272.       elsif RT = Standard_Float then
  273.          return UI_From_Int (Float_Attr_Machine_Mantissa);
  274.  
  275.       elsif RT = Standard_Long_Float then
  276.          return UI_From_Int (Long_Float_Attr_Machine_Mantissa);
  277.  
  278.       else
  279.          pragma Assert (RT = Standard_Long_Long_Float);
  280.          return UI_From_Int (Long_Long_Float_Attr_Machine_Mantissa);
  281.       end if;
  282.    end Machine_Mantissa;
  283.  
  284.    -----------
  285.    -- Model --
  286.    -----------
  287.  
  288.    function Model (RT : R; X : T) return T is
  289.       X_Frac : T;
  290.       X_Exp  : UI;
  291.  
  292.    begin
  293.       Decompose (RT, X, X_Frac, X_Exp);
  294.       return Compose (RT, X_Frac, X_Exp);
  295.    end Model;
  296.  
  297.    ----------
  298.    -- Pred --
  299.    ----------
  300.  
  301.    function Pred (RT : R; X : T) return T is
  302.       X_Frac : T;
  303.       X_Exp  : UI;
  304.       Mach   : T;
  305.  
  306.    begin
  307.       Decompose (RT, X, X_Frac, X_Exp);
  308.       Mach := Compose (RT, X_Frac, X_Exp);
  309.  
  310.       --  If number was not a machine number, then just return the machine
  311.       --  number below it (compose always returns the machine number below)
  312.  
  313.       if X /= Mach then
  314.          return Mach;
  315.  
  316.       else
  317.          --  Subtract from the given number a number equivalent to the value
  318.          --  of its least significant bit. Given that the most significant bit
  319.          --  represents a value of 1.0 * radix ** (exp - 1), the value we want
  320.          --  is obtained by shifting this by (mantissa-1) bits to the right,
  321.          --  i.e. decreasing the exponent by that amount.
  322.  
  323.          return X - Scaling (RT, Ureal_1, X_Exp - Machine_Mantissa (RT));
  324.       end if;
  325.    end Pred;
  326.  
  327.    ------------------------
  328.    -- Radix_To_M_Minus_1 --
  329.    ------------------------
  330.  
  331.    function Radix_To_M_Minus_1 (RT : R) return T is
  332.    begin
  333.       return Float_Radix ** (Machine_Mantissa (RT) - 1);
  334.    end Radix_To_M_Minus_1;
  335.  
  336.    ---------------
  337.    -- Remainder --
  338.    ---------------
  339.  
  340.    function Remainder (RT : R; X, Y : T) return T is
  341.       A        : T;
  342.       B        : T;
  343.       Arg      : T;
  344.       P        : T;
  345.       Arg_Frac : T;
  346.       P_Frac   : T;
  347.       Sign_X   : T;
  348.       IEEE_Rem : T;
  349.       Arg_Exp  : UI;
  350.       P_Exp    : UI;
  351.       K        : UI;
  352.       P_Even   : Boolean;
  353.  
  354.    begin
  355.       if UR_Is_Positive (X) then
  356.          Sign_X :=  Ureal_1;
  357.       else
  358.          Sign_X := -Ureal_1;
  359.       end if;
  360.  
  361.       Arg := abs X;
  362.       P   := abs Y;
  363.  
  364.       if Arg < P then
  365.          P_Even := True;
  366.          IEEE_Rem := Arg;
  367.          P_Exp := Exponent (RT, P);
  368.  
  369.       else
  370.          Decompose (RT, Arg, Arg_Frac, Arg_Exp);
  371.          Decompose (RT, P,   P_Frac,   P_Exp);
  372.  
  373.          P := Compose (RT, P_Frac, Arg_Exp);
  374.          K := Arg_Exp - P_Exp;
  375.          P_Even := True;
  376.          IEEE_Rem := Arg;
  377.  
  378.          for Cnt in reverse 0 .. UI_To_Int (K) loop
  379.             if IEEE_Rem >= P then
  380.                P_Even := False;
  381.                IEEE_Rem := IEEE_Rem - P;
  382.             else
  383.                P_Even := True;
  384.             end if;
  385.  
  386.             P := P * Ureal_Half;
  387.          end loop;
  388.       end if;
  389.  
  390.       --  That completes the calculation of modulus remainder. The final step
  391.       --  is get the IEEE remainder. Here we compare Rem with (abs Y) / 2.
  392.  
  393.       if P_Exp >= 0 then
  394.          A := IEEE_Rem;
  395.          B := abs Y * Ureal_Half;
  396.  
  397.       else
  398.          A := IEEE_Rem * Ureal_2;
  399.          B := abs Y;
  400.       end if;
  401.  
  402.       if A > B or else (A = B and then not P_Even) then
  403.          IEEE_Rem := IEEE_Rem - abs Y;
  404.       end if;
  405.  
  406.       return Sign_X * IEEE_Rem;
  407.  
  408.    end Remainder;
  409.  
  410.    --------------
  411.    -- Rounding --
  412.    --------------
  413.  
  414.    function Rounding (RT : R; X : T) return T is
  415.       Result : T;
  416.       Tail   : T;
  417.  
  418.    begin
  419.       Result := Truncation (RT, abs X);
  420.       Tail   := abs X - Result;
  421.  
  422.       if Tail > Ureal_Half  then
  423.          Result := Result + Ureal_1;
  424.       end if;
  425.  
  426.       if UR_Is_Negative (X) then
  427.          return -Result;
  428.       else
  429.          return Result;
  430.       end if;
  431.  
  432.    end Rounding;
  433.  
  434.    -------------
  435.    -- Scaling --
  436.    -------------
  437.  
  438.    function Scaling (RT : R; X : T; Adjustment : UI) return T is
  439.    begin
  440.       return X * Uint_2 ** Adjustment;
  441.    end Scaling;
  442.  
  443.    ----------
  444.    -- Succ --
  445.    ----------
  446.  
  447.    function Succ (RT : R; X : T) return T is
  448.       X_Frac : T;
  449.       X_Exp  : UI;
  450.  
  451.    begin
  452.       Decompose (RT, X, X_Frac, X_Exp);
  453.  
  454.       --  Note: this gives either X or the machine number just below it, so
  455.       --  we want the macnine number above this in either case. Compuation
  456.       --  of this is similar to that in Pred, except that we add the value
  457.       --  of the least significant bit instead of subtracting it.
  458.  
  459.       return X + Scaling (RT, Ureal_1, X_Exp - Machine_Mantissa (RT));
  460.  
  461.    end Succ;
  462.  
  463.    ----------------
  464.    -- Truncation --
  465.    ----------------
  466.  
  467.    function Truncation (RT : R; X : T) return T is
  468.    begin
  469.       return UR_From_Uint (UR_Trunc (X));
  470.    end Truncation;
  471.  
  472.    -----------------------
  473.    -- Unbiased_Rounding --
  474.    -----------------------
  475.  
  476.    function Unbiased_Rounding (RT : R; X : T) return T is
  477.       Result : T;
  478.       Tail   : T;
  479.  
  480.    begin
  481.       Result := Truncation (RT, abs X);
  482.       Tail   := abs X - Result;
  483.  
  484.       if Tail > Ureal_Half  then
  485.          Result := Result + Ureal_1;
  486.  
  487.       elsif Tail = Ureal_Half then
  488.          Result := Ureal_2 * Rounding (RT, Result / Ureal_2);
  489.       end if;
  490.  
  491.       if UR_Is_Negative (X) then
  492.          return -Result;
  493.       else
  494.          return Result;
  495.       end if;
  496.  
  497.    end Unbiased_Rounding;
  498.  
  499. end Eval_Fat;
  500.