Re: how to round integers (Figured it out!)

["Robert I. Eachus" <rieachus@attbi.com>]

Jeffrey Carter wrote:

> Then your original algorithm would be fine.

I would say my elaboration of David C. Hoos version.  But it was a point 
worth making, and since I have written this sort of stuff several times, 
might as well put up the "right" answer for the general case.

But what I really should do is write the one-liner package that converts 
to float to do the divide, then converts back to the integer type (done):

   generic
     type Int is range <>;
   function Rounding2(Dividend, Divisor: Int) return Int;

   function Rounding2(Dividend, Divisor: Int) return Int is
   begin
     return Int(Long_Float(Dividend)/Long_Float(Divisor));
   end Rounding2;

Note that this package will round away from zero.  But the interesting 
thing to do of course, is compare the two packages on timing (program at 
end):

E:\Ada\Bug>time_rounding
time_rounding
  Iteration  1 Round = 49.000 Round2 = 9.640
  Iteration  2 Round = 48.979 Round2 = 9.658
  Iteration  3 Round = 49.010 Round2 = 9.661
  Iteration  4 Round = 49.001 Round2 = 9.643
  Iteration  5 Round = 49.381 Round2 = 9.645

  Averages     Round = 49.074 Round2 = 9.650

First note that your milage may definitely vary, this was on a 700 MHz 
Athlon Classic.  Second, the program as written will sit and churn, in 
my case for over 5 minutes.  Third, no those results are not surprising. 
  By using the floating point divide, ON THIS PROCESSOR, the divides 
have one execution pipe to themselves, which is unused in the other 
version.  Other architectures have different designs.  For example some 
use the same pipe for integer and floating point divides.  And finally 
version 2 works out to 2.59 iterations per microsecond.  Not bad for all 
the work done there.  (I didn't play with other options, including 
unrolling loops and inlining the code.)

-----------------------------------------------------------------------
with Ada.Text_IO; use Ada.Text_IO;
with Ada.Calendar; use Ada.Calendar;
with Ada.Numerics.Discrete_Random;
procedure Time_Rounding is

   subtype Int is Integer range -Integer'Last..Integer'Last;
   package Int_Rand is new Ada.Numerics.Discrete_Random(Int);
   -- I want random numbers that are closed over division.
   -- Integer'First/-1 overflows.

   Gen: Int_Rand.Generator;
   function Rand(G: in Int_Rand.Generator := Gen)
       return Integer renames Int_Rand.Random;

   type Int_Array is array(1..5_000) of Integer;
   A,B, Results : Int_Array;
   Temp: Integer;
   Maximum: Integer;
   Number_Of_Iterations: constant Integer := 5;
   Times: array(Integer range 1..2, 1..Number_Of_Iterations) of Duration;

   Start, Stop: Time;

   package Duration_IO is new Fixed_IO(Duration);

   procedure Put_Time(Seconds: Duration;
             Fore: Field := 1;
             Aft: Field := 3;
             Exp: Field := 0) renames Duration_IO.Put;

   generic
     type Int is range <>;
   function Rounding(Dividend, Divisor: Int) return Int;

   generic
     type Int is range <>;
   function Rounding2(Dividend, Divisor: Int) return Int;

   function Rounding(Dividend, Divisor: Int) return Int is
   begin
     if (Dividend >= 0 and Divisor > 0) or else
        (Dividend <= 0 and Divisor < 0)
     then return (Dividend + Dividend + Divisor) /
           (Divisor + Divisor);
     elsif Dividend <= 0 and Divisor > 0
     then return (Dividend + Dividend - Divisor + 1) /
          (Divisor + Divisor);
     else return (Dividend + Dividend - Divisor - 1) /
          (Divisor + Divisor);
     end if;
   exception
     when Constraint_Error =>
       if Divisor = 0 then raise Constraint_Error; end if;
       -- get the case we can't handle out of the way...
       declare
         Temp: Int;
       begin
         Temp := abs Dividend rem Divisor;
         if Divisor < 0 and Dividend > 0 then
           if Temp > -(Int'First/2)
               or else (Temp-1) * 2 >= -(Divisor + 1)
           then return Dividend/Divisor - 1;
           else return Dividend/Divisor;
           end if;
         elsif Divisor < 0 and Dividend < 0 then
           if Temp > -(Int'First/2)
               or else (Temp-1) * 2 >= -(Divisor + 2)
           then return Dividend/Divisor + 1;
           else return Dividend/Divisor;
           end if;
         elsif Divisor > 0 and Dividend > 0 then
           if Temp > Int'Last/2 or else Temp * 2 >= Divisor
           then return Dividend/Divisor + 1;
           else return Dividend/Divisor;
           end if;
        else -- Divisor > 0 and Dividend < 0 then
           if Temp > Int'Last/2 or else (Temp-1) * 2 >= Divisor - 1
           then return Dividend/Divisor - 1;
           else return Dividend/Divisor;
           end if;
         end if;
       end;
   end Rounding;

   function Round is new Rounding(Integer);

   function Rounding2(Dividend, Divisor: Int) return Int is
   begin
     return Int(Long_Float(Dividend)/Long_Float(Divisor));
   end Rounding2;

   function Round2 is new Rounding2(Integer);

begin

  -- Initialize Random Number Generator;

  Int_Rand.Reset(Gen, 1234); --  I want repeatable results.

  -- Initialize Arrays

   for I in Int_Array'Range loop
     A(I) := Rand;
     B(I) := Rand;
     while B(I) = 0 loop
       B(I) := Rand;
     end loop; -- Don't want divide by zero errors.

   end loop;

  -- Test rounding algorithms, and as a side effect load caches.

   for I in Int_Array'Range loop
     Maximum := Integer'First;
     for J in Int_Array'Range loop
       Temp := Round(A(I), B(J));
       if Round2(A(I), B(J)) /= Temp then
         Put_Line("For " & Integer'Image(A(I)) & Integer'Image(B(I)) & " 
the results are "
            & Integer'Image(Temp) & " and " & Integer'Image(Round2(A(I), 
B(J))) & '.');

       -- There could be some differences.  Assuming the values are 
independent and
       -- uniformly distributed the rounding will be different if:
       -- 1. A(I) and B(J) have opposite signs.
       -- 2. A(I) is even.
       -- 3. B(J) mod A(I) = A(I)/2;
       -- The first two constraints eliminate 3/4 of the cases.  For the 
others, there
       -- will be Integer'Last/A(I) + Integer'Last mod A(I) > A(I)/2 
(math, not Ada!)
       -- Easy enough to write a program to get the probability per 
pair: 1.16415E-10.
       -- Of course, if I wanted to check the computations, I could 
choose A & B in a
       -- way that increases the odds.  For example, make B small and 
always even.

       end if;
       if Temp > Maximum then Maximum := Temp; end if;
     end loop;
     Results(I) := Maximum;
   end loop;

   -- Now do it for timings.

   for Iteration in 1..Number_Of_Iterations loop
     Start := Clock;
       for I in Int_Array'Range loop
       Maximum := Integer'First;
       for J in Int_Array'Range loop
         Temp := Round(A(I), B(J));
         if Temp > Maximum then Maximum := Temp; end if;
       end loop;
       Results(I) := Maximum;
     end loop;
     Stop := Clock;
     Times(1, Iteration) := Stop-Start;

     Start := Clock;
       for I in Int_Array'Range loop
       Maximum := Integer'First;
       for J in Int_Array'Range loop
         Temp := Round2(A(I), B(J));
         if Temp > Maximum then Maximum := Temp; end if;
       end loop;
       Results(I) := Maximum;
     end loop;
     Stop := Clock;
     Times(2, Iteration) := Stop-Start;

   end loop;

   declare
     Sum1, Sum2: Duration := 0.0;
   begin
     for I in 1..Number_Of_Iterations loop
       Sum1 := Sum1 + Times(1,I);
       Sum2 := Sum2 + Times(2,I);
       Put(" Iteration " & Integer'Image(I) & " Round = "); 
Put_Time(Times(1,I));
       Put(" Round2 = "); Put_Time(Times(2,I)); New_Line;
     end loop;
     New_Line;

     Put(" Averages    " & " Round = "); 
Put_Time(Sum1/Number_Of_Iterations);
     Put(" Round2 = "); Put_Time(Sum2/Number_Of_Iterations); New_Line;
   end;

end Time_Rounding;