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 / s-arit64.adb < prev    next >
Text File  |  1996-09-28  |  17KB  |  667 lines
  1. ------------------------------------------------------------------------------
  2. --                                                                          --
  3. --                         GNAT COMPILER COMPONENTS                         --
  4. --                                                                          --
  5. --                      S Y S T E M . A R I T H _ 6 4                       --
  6. --                                                                          --
  7. --                                 B o d y                                  --
  8. --                                                                          --
  9. --                            $Revision: 1.2 $                              --
  10. --                                                                          --
  11. -- The GNAT library is free software; you can redistribute it and/or modify --
  12. -- it under terms of the GNU Library General Public License as published by --
  13. -- the Free Software  Foundation; either version 2, or (at your option) any --
  14. -- later version.  The GNAT library is distributed in the hope that it will --
  15. -- be useful, but WITHOUT ANY WARRANTY;  without even  the implied warranty --
  16. -- of MERCHANTABILITY  or  FITNESS FOR  A PARTICULAR PURPOSE.  See the  GNU --
  17. -- Library  General  Public  License for  more  details.  You  should  have --
  18. -- received  a copy of the GNU  Library  General Public License  along with --
  19. -- the GNAT library;  see the file  COPYING.LIB.  If not, write to the Free --
  20. -- Free Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA.        --
  21. --                                                                          --
  22. ------------------------------------------------------------------------------
  23.  
  24. with Text_IO; use Text_IO;
  25. with Interfaces; use Interfaces;
  26. with Unchecked_Conversion;
  27.  
  28. package body System.Arith_64 is
  29.  
  30.    pragma Suppress (Overflow_Check);
  31.    pragma Suppress (Range_Check);
  32.  
  33.    subtype Uns64 is Unsigned_64;
  34.    function To_Uns is new Unchecked_Conversion (Int64, Uns64);
  35.    function To_Int is new Unchecked_Conversion (Uns64, Int64);
  36.  
  37.    subtype Uns32 is Unsigned_32;
  38.  
  39.    -----------------------
  40.    -- Local Subprograms --
  41.    -----------------------
  42.  
  43.    function "+" (A, B : Uns32) return Uns64;
  44.    function "+" (A : Uns64; B : Uns32) return Uns64;
  45.    pragma Inline ("+");
  46.    --  Length doubling additions
  47.  
  48.    function "-" (A : Uns64; B : Uns32) return Uns64;
  49.    pragma Inline ("-");
  50.    --  Length doubling subtraction
  51.  
  52.    function "*" (A, B : Uns32) return Uns64;
  53.    function "*" (A : Uns64; B : Uns32) return Uns64;
  54.    pragma Inline ("*");
  55.    --  Length doubling multiplications
  56.  
  57.    function "/" (A : Uns64; B : Uns32) return Uns64;
  58.    pragma Inline ("/");
  59.    --  Length doubling division
  60.  
  61.    function "rem" (A : Uns64; B : Uns32) return Uns64;
  62.    pragma Inline ("rem");
  63.    --  Length doubling remainder
  64.  
  65.    function "&" (Hi, Lo : Uns32) return Uns64;
  66.    pragma Inline ("&");
  67.    --  Concatenate hi, lo values to form 64-bit result
  68.  
  69.    function Lo (A : Uns64) return Uns32;
  70.    pragma Inline (Lo);
  71.    --  Low order half of 64-bit value
  72.  
  73.    function Hi (A : Uns64) return Uns32;
  74.    pragma Inline (Hi);
  75.    --  High order half of 64 bit value
  76.  
  77.    ---------
  78.    -- "+" --
  79.    ---------
  80.  
  81.    function "+" (A, B : Uns32) return Uns64 is
  82.    begin
  83.       return Uns64 (A) + Uns64 (B);
  84.    end "+";
  85.  
  86.    function "+" (A : Uns64; B : Uns32) return Uns64 is
  87.    begin
  88.       return A + Uns64 (B);
  89.    end "+";
  90.  
  91.    ---------
  92.    -- "-" --
  93.    ---------
  94.  
  95.    function "-" (A : Uns64; B : Uns32) return Uns64 is
  96.    begin
  97.       return A - Uns64 (B);
  98.    end "-";
  99.  
  100.    ---------
  101.    -- "*" --
  102.    ---------
  103.  
  104.    function "*" (A, B : Uns32) return Uns64 is
  105.    begin
  106.       return Uns64 (A) * Uns64 (B);
  107.    end "*";
  108.  
  109.    function "*" (A : Uns64; B : Uns32) return Uns64 is
  110.    begin
  111.       return A * Uns64 (B);
  112.    end "*";
  113.  
  114.    ---------
  115.    -- "/" --
  116.    ---------
  117.  
  118.    function "/" (A : Uns64; B : Uns32) return Uns64 is
  119.    begin
  120.       return A / Uns64 (B);
  121.    end "/";
  122.  
  123.    -----------
  124.    -- "rem" --
  125.    -----------
  126.  
  127.    function "rem" (A : Uns64; B : Uns32) return Uns64 is
  128.    begin
  129.       return A rem Uns64 (B);
  130.    end "rem";
  131.  
  132.    ---------
  133.    -- "&" --
  134.    ---------
  135.  
  136.    function "&" (Hi, Lo : Uns32) return Uns64 is
  137.    begin
  138.       return Shift_Left (Uns64 (Hi), 32) or Uns64 (Lo);
  139.    end "&";
  140.  
  141.    --------
  142.    -- Hi --
  143.    --------
  144.  
  145.    function Hi (A : Uns64) return Uns32 is
  146.    begin
  147.       return Uns32 (Shift_Right (A, 32));
  148.    end Hi;
  149.  
  150.    --------
  151.    -- Lo --
  152.    --------
  153.  
  154.    function Lo (A : Uns64) return Uns32 is
  155.    begin
  156.       return Uns32 (A and 16#FFFF_FFFF#);
  157.    end Lo;
  158.  
  159.    --------------------------
  160.    -- Add_With_Ovflo_Check --
  161.    --------------------------
  162.  
  163.    function Add_With_Ovflo_Check (X, Y : Int64) return Int64 is
  164.       R : constant Int64 := To_Int (To_Uns (X) + To_Uns (Y));
  165.  
  166.    begin
  167.       if X >= 0 then
  168.          if Y < 0 or else R >= 0 then
  169.             return R;
  170.          end if;
  171.  
  172.       else -- X < 0
  173.          if Y > 0 or else R < 0 then
  174.             return R;
  175.          end if;
  176.       end if;
  177.  
  178.       raise Constraint_Error;
  179.    end Add_With_Ovflo_Check;
  180.  
  181.    -----------------------------
  182.    -- Divide_With_Ovflo_Check --
  183.    -----------------------------
  184.  
  185.    function Divide_With_Ovflo_Check (X, Y : Int64) return Int64 is
  186.    begin
  187.       if Y = 0
  188.         or else (Y = (-1) and then X = (-(2 ** 63)))
  189.       then
  190.          raise Constraint_Error;
  191.       else
  192.          return X / Y;
  193.       end if;
  194.    end Divide_With_Ovflo_Check;
  195.  
  196.    -------------------
  197.    -- Double_Divide --
  198.    -------------------
  199.  
  200.    procedure Double_Divide
  201.      (X, Y, Z : Int64;
  202.       Q, R    : out Int64;
  203.       Round   : Boolean)
  204.    is
  205.       Xu  : constant Uns64 := To_Uns (abs X);
  206.       Xhi : constant Uns32 := Hi (Xu);
  207.       Xlo : constant Uns32 := Lo (Xu);
  208.  
  209.       Yu  : constant Uns64 := To_Uns (abs Y);
  210.       Yhi : constant Uns32 := Hi (Yu);
  211.       Ylo : constant Uns32 := Lo (Yu);
  212.  
  213.       Zu  : constant Uns64 := To_Uns (abs Z);
  214.       Zhi : constant Uns32 := Hi (Zu);
  215.       Zlo : constant Uns32 := Lo (Zu);
  216.  
  217.       T1, T2     : Uns64;
  218.       Du, Qu, Ru : Uns64;
  219.       Den_Pos    : Boolean;
  220.  
  221.    begin
  222.       if Yu = 0 or else Zu = 0 then
  223.          raise Constraint_Error;
  224.       end if;
  225.  
  226.       --  Compute Y * Z. Note that if the result overflows 64 bits unsigned,
  227.       --  then the rounded result is clearly zero (since the dividend is at
  228.       --  most 2**63 - 1, the extra bit of precision is nice here!)
  229.  
  230.       if Yhi /= 0 then
  231.          if Zhi /= 0 then
  232.             Q := 0;
  233.             R := X;
  234.             return;
  235.          else
  236.             T2 := Yhi * Zlo;
  237.          end if;
  238.  
  239.       else
  240.          if Zhi /= 0 then
  241.             T2 := Ylo * Zhi;
  242.          else
  243.             T2 := 0;
  244.          end if;
  245.       end if;
  246.  
  247.       T1 := Xlo * Ylo;
  248.       T2 := T2 + Hi (T1);
  249.  
  250.       if Hi (T2) /= 0 then
  251.          Q := 0;
  252.          R := X;
  253.          return;
  254.       end if;
  255.  
  256.       Du := Lo (T2) & Lo (T1);
  257.       Qu := Xu / Du;
  258.       Ru := Xu rem Du;
  259.  
  260.       --  Deal with rounding case
  261.  
  262.       if Round and then Ru > (Du - 1) / 2 then
  263.          Qu := Qu + 1;
  264.       end if;
  265.  
  266.       --  Set final signs (RM 4.5.5(27-30))
  267.  
  268.       Den_Pos := (Y < 0) = (Z < 0);
  269.  
  270.       --  Case of dividend (X) sign positive
  271.  
  272.       if X >= 0 then
  273.          R := To_Int (Ru);
  274.  
  275.          if Den_Pos then
  276.             Q := To_Int (Qu);
  277.          else
  278.             Q := -To_Int (Qu);
  279.          end if;
  280.  
  281.       --  Case of dividend (X) sign negative
  282.  
  283.       else
  284.          R := -To_Int (Ru);
  285.  
  286.          if Den_Pos then
  287.             Q := -To_Int (Qu);
  288.          else
  289.             Q := To_Int (Qu);
  290.          end if;
  291.       end if;
  292.    end Double_Divide;
  293.  
  294.    -------------------------------
  295.    -- Multiply_With_Ovflo_Check --
  296.    -------------------------------
  297.  
  298.    function Multiply_With_Ovflo_Check (X, Y : Int64) return Int64 is
  299.       Xu  : constant Uns64 := To_Uns (abs X);
  300.       Xhi : constant Uns32 := Hi (Xu);
  301.       Xlo : constant Uns32 := Lo (Xu);
  302.  
  303.       Yu  : constant Uns64 := To_Uns (abs Y);
  304.       Yhi : constant Uns32 := Hi (Yu);
  305.       Ylo : constant Uns32 := Lo (Yu);
  306.  
  307.       T1, T2 : Uns64;
  308.  
  309.    begin
  310.       if Xhi /= 0 then
  311.          if Yhi /= 0 then
  312.             raise Constraint_Error;
  313.          else
  314.             T2 := Xhi * Ylo;
  315.          end if;
  316.  
  317.       else
  318.          if Yhi /= 0 then
  319.             T2 := Xlo * Yhi;
  320.          else
  321.             return X * Y;
  322.          end if;
  323.       end if;
  324.  
  325.       T1 := Xlo * Ylo;
  326.       T2 := T2 + Hi (T1);
  327.  
  328.       if Hi (T2) /= 0 then
  329.          raise Constraint_Error;
  330.       end if;
  331.  
  332.       T2 := Lo (T2) & Lo (T1);
  333.  
  334.       if X >= 0 then
  335.          if Y >= 0 then
  336.             return To_Int (T2);
  337.          else
  338.             return -To_Int (T2);
  339.          end if;
  340.       else -- X < 0
  341.          if Y < 0 then
  342.             return To_Int (T2);
  343.          else
  344.             return -To_Int (T2);
  345.          end if;
  346.       end if;
  347.  
  348.    end Multiply_With_Ovflo_Check;
  349.  
  350.    -------------------
  351.    -- Scaled_Divide --
  352.    -------------------
  353.  
  354.    procedure Scaled_Divide
  355.      (X, Y, Z : Int64;
  356.       Q, R    : out Int64;
  357.       Round   : Boolean)
  358.    is
  359.       Xu  : constant Uns64 := To_Uns (abs X);
  360.       Xhi : constant Uns32 := Hi (Xu);
  361.       Xlo : constant Uns32 := Lo (Xu);
  362.  
  363.       Yu  : constant Uns64 := To_Uns (abs Y);
  364.       Yhi : constant Uns32 := Hi (Yu);
  365.       Ylo : constant Uns32 := Lo (Yu);
  366.  
  367.       Zu  : Uns64 := To_Uns (abs Z);
  368.       Zhi : Uns32 := Hi (Zu);
  369.       Zlo : Uns32 := Lo (Zu);
  370.  
  371.       D1, D2, D3, D4 : Uns32;
  372.       --  The dividend, four digits (D1 is high order)
  373.  
  374.       Q1, Q2 : Uns32;
  375.       --  The quotient, two digits (Q1 is high order)
  376.  
  377.       S1, S2, S3 : Uns32;
  378.       --  Value to subtract, three digits (S1 is high order)
  379.  
  380.       Qu : Uns64;
  381.       Ru : Uns64;
  382.       --  Unsigned quotient and remainder
  383.  
  384.       Scale : Natural;
  385.       --  Scaling factor used for multiple-precision divide. Dividend and
  386.       --  Divisor are multiplied by 2 ** Scale, and the final remainder
  387.       --  is divided by the scaling factor. The reason for this scaling
  388.       --  is to allow more accurate estimation of quotient digits.
  389.  
  390.       T1, T2, T3 : Uns64;
  391.       --  Temporary values
  392.  
  393.    begin
  394.       --  First do the multiplication, giving the four digit dividend
  395.  
  396.       T1 := Xlo * Ylo;
  397.       D4 := Lo (T1);
  398.       D3 := Hi (T1);
  399.  
  400.       if Yhi /= 0 then
  401.          T1 := Xlo * Yhi;
  402.          T2 := D3 + Lo (T1);
  403.          D3 := Lo (T2);
  404.          D2 := Hi (T1) + Hi (T2);
  405.  
  406.          if Xhi /= 0 then
  407.             T1 := Xhi * Ylo;
  408.             T2 := D3 + Lo (T1);
  409.             D3 := Lo (T2);
  410.             T3 := D2 + Hi (T1);
  411.             T3 := T3 + Hi (T2);
  412.             D2 := Lo (T3);
  413.             D1 := Hi (T3);
  414.  
  415.             T1 := (D1 & D2) + (Xhi * Yhi);
  416.             D1 := Hi (T1);
  417.             D2 := Lo (T1);
  418.  
  419.          else
  420.             D1 := 0;
  421.          end if;
  422.  
  423.       else
  424.          if Xhi /= 0 then
  425.             T1 := Xhi * Ylo;
  426.             T2 := D3 + Lo (T1);
  427.             D3 := Lo (T2);
  428.             D2 := Hi (T1) + Hi (T2);
  429.  
  430.          else
  431.             D2 := 0;
  432.          end if;
  433.  
  434.          D1 := 0;
  435.       end if;
  436.  
  437.       --  Now it is time for the dreaded multiple precision division. First
  438.       --  an easy case, check for the simple case of a one digit divisor.
  439.  
  440.       if Zhi = 0 then
  441.          if D1 /= 0 or else D2 >= Zlo then
  442.             raise Constraint_Error;
  443.  
  444.          --  Here we are dividing at most three digits by one digit
  445.  
  446.          else
  447.             T1 := D2 & D3;
  448.             T2 := Lo (T1 rem Zlo) & D4;
  449.  
  450.             Qu := Lo (T1 / Zlo) & Lo (T2 / Zlo);
  451.             Ru := T2 rem Zlo;
  452.          end if;
  453.  
  454.       --  If divisor is double digit and too large, raise error
  455.  
  456.       elsif (D1 & D2) >= Zu then
  457.          raise Constraint_Error;
  458.  
  459.       --  This is the complex case where we definitely have a double digit
  460.       --  divisor and a dividend of at least three digits. We use the classical
  461.       --  multiple division algorithm (see  section (4.3.1) of Knuth's "The Art
  462.       --  of Computer Programming", Vol. 2 for a description (algorithm D).
  463.  
  464.       else
  465.          --  First normalize the divisor so that it has the leading bit on.
  466.          --  We do this by finding the appropriate left shift amount.
  467.  
  468.          Scale := 0;
  469.  
  470.          if (Zhi and 16#FFFF0000#) = 0 then
  471.             Scale := 16;
  472.             Zu := Shift_Left (Zu, 16);
  473.          end if;
  474.  
  475.          if (Hi (Zu) and 16#FF00_0000#) = 0 then
  476.             Scale := Scale + 8;
  477.             Zu := Shift_Left (Zu, 8);
  478.          end if;
  479.  
  480.          if (Hi (Zu) and 16#F000_0000#) = 0 then
  481.             Scale := Scale + 4;
  482.             Zu := Shift_Left (Zu, 4);
  483.          end if;
  484.  
  485.          if (Hi (Zu) and 16#C000_0000#) = 0 then
  486.             Scale := Scale + 2;
  487.             Zu := Shift_Left (Zu, 2);
  488.          end if;
  489.  
  490.          if (Hi (Zu) and 16#8000_0000#) = 0 then
  491.             Scale := Scale + 1;
  492.             Zu := Shift_Left (Zu, 1);
  493.          end if;
  494.  
  495.          Zhi := Hi (Zu);
  496.          Zlo := Lo (Zu);
  497.  
  498.          --  Note that when we scale up the dividend, it still fits in four
  499.          --  digits, since we already tested for overflow, and scaling does
  500.          --  not change the invariant that (D1 & D2) >= Zu.
  501.  
  502.          T1 := Shift_Left (D1 & D2, Scale);
  503.          D1 := Hi (T1);
  504.          T2 := Shift_Left (0 & D3, Scale);
  505.          D2 := Lo (T1) or Hi (T2);
  506.          T3 := Shift_Left (0 & D4, Scale);
  507.          D3 := Lo (T2) or Hi (T3);
  508.          D4 := Lo (T3);
  509.  
  510.          --  Compute first quotient digit. We have to divide three digits by
  511.          --  two digits, and we estimate the quotient by dividing the leading
  512.          --  two digits by the leading digit. Given the scaling we did above
  513.          --  which ensured the first bit of the divisor is set, this gives an
  514.          --  estimate of the quotient that is at most two too high.
  515.  
  516.          if D1 = Zhi then
  517.             Q1 := 2 ** 32 - 1;
  518.          else
  519.             Q1 := Lo ((D1 & D2) / Zhi);
  520.          end if;
  521.  
  522.          --  Compute amount to subtract
  523.  
  524.          T1 := Q1 * Zlo;
  525.          T2 := Q1 * Zhi;
  526.          S3 := Lo (T1);
  527.          T1 := Hi (T1) + Lo (T2);
  528.          S2 := Lo (T1);
  529.          S1 := Hi (T1) + Hi (T2);
  530.  
  531.          --  Adjust quotient digit if it was too high
  532.  
  533.          loop
  534.             exit when S1 < D1;
  535.  
  536.             if S1 = D1 then
  537.                exit when S2 < D2;
  538.  
  539.                if S2 = D2 then
  540.                   exit when S3 <= D3;
  541.                end if;
  542.             end if;
  543.  
  544.             Q1 := Q1 - 1;
  545.  
  546.             T1 := (S2 & S3) - Zlo;
  547.             S3 := Lo (T1);
  548.             T1 := (S1 & S2) - Zhi;
  549.             S2 := Lo (T1);
  550.             S1 := Hi (T1);
  551.          end loop;
  552.  
  553.          --  Subtract from dividend (note: do not bother to set D1 to
  554.          --  zero, since it is no longer needed in the calculation).
  555.  
  556.          T1 := (D2 & D3) - S3;
  557.          D3 := Lo (T1);
  558.          T1 := (D1 & Hi (T1)) - S2;
  559.          D2 := Lo (T1);
  560.  
  561.          --  Compute second quotient digit in same manner
  562.  
  563.          if D2 = Zhi then
  564.             Q2 := 2 ** 32 - 1;
  565.          else
  566.             Q2 := Lo ((D2 & D3) / Zhi);
  567.          end if;
  568.  
  569.          T1 := Q2 * Zlo;
  570.          T2 := Q2 * Zhi;
  571.          S3 := Lo (T1);
  572.          T1 := Hi (T1) + Lo (T2);
  573.          S2 := Lo (T1);
  574.          S1 := Hi (T1) + Hi (T2);
  575.  
  576.          loop
  577.             exit when S1 < D2;
  578.  
  579.             if S1 = D2 then
  580.                exit when S2 < D3;
  581.  
  582.                if S2 = D3 then
  583.                   exit when S3 <= D4;
  584.                end if;
  585.             end if;
  586.  
  587.             Q2 := Q2 - 1;
  588.  
  589.             T1 := (S2 & S3) - Zlo;
  590.             S3 := Lo (T1);
  591.             T1 := (S1 & S2) - Zhi;
  592.             S2 := Lo (T1);
  593.             S1 := Hi (T1);
  594.          end loop;
  595.  
  596.          T1 := (D3 & D4) - S3;
  597.          D4 := Lo (T1);
  598.          T1 := (D2 & Hi (T1)) - S2;
  599.          D3 := Lo (T1);
  600.  
  601.          --  The two quotient digits are now set, and the remainder of the
  602.          --  scaled division is in (D3 & D4). To get the remainder for the
  603.          --  original unscaled division, we rescale this dividend.
  604.  
  605.          Qu := Q1 & Q2;
  606.          Ru := Shift_Right (D3 & D4, Scale);
  607.       end if;
  608.  
  609.       --  Deal with rounding case
  610.  
  611.       if Round and then Ru > (Zu - 1) / 2 then
  612.          Qu := Qu + 1;
  613.       end if;
  614.  
  615.       --  Set final signs (RM 4.5.5(27-30))
  616.  
  617.       --  Case of dividend (X * Y) sign positive
  618.  
  619.       if (X >= 0 and then Y >= 0)
  620.         or else (X < 0 and then Y < 0)
  621.       then
  622.          R := To_Int (Ru);
  623.  
  624.          if Z > 0 then
  625.             Q := To_Int (Qu);
  626.          else
  627.             Q := -To_Int (Qu);
  628.          end if;
  629.  
  630.       --  Case of dividend (X * Y) sign negative
  631.  
  632.       else
  633.          R := -To_Int (Ru);
  634.  
  635.          if Z > 0 then
  636.             Q := -To_Int (Qu);
  637.          else
  638.             Q := To_Int (Qu);
  639.          end if;
  640.       end if;
  641.  
  642.    end Scaled_Divide;
  643.  
  644.    -------------------------------
  645.    -- Subtract_With_Ovflo_Check --
  646.    -------------------------------
  647.  
  648.    function Subtract_With_Ovflo_Check (X, Y : Int64) return Int64 is
  649.       R : constant Int64 := To_Int (To_Uns (X) - To_Uns (Y));
  650.  
  651.    begin
  652.       if X >= 0 then
  653.          if Y > 0 or else R >= 0 then
  654.             return R;
  655.          end if;
  656.  
  657.       else -- X < 0
  658.          if Y <= 0 or else R < 0 then
  659.             return R;
  660.          end if;
  661.       end if;
  662.  
  663.       raise Constraint_Error;
  664.    end Subtract_With_Ovflo_Check;
  665.  
  666. end System.Arith_64;
  667.