Red-faced professor gets bitten in search for portability

Red-faced professor gets bitten in search for portability

[Michael Feldman]

Asked by students whether Ada has an equivalent of the Pascal "Trunc" 
operation, which just returns the integer part of its argument, 
I put together a function based on the old trick we used
in the Fortran days.

  FUNCTION Trunc (X: Float) RETURN Integer IS
  BEGIN
    RETURN Integer(X - 0.5);
  END Trunc;

Since conversion of Float to Integer is, in Ada, a rounding operation,
this looks like a good solution, right? WRONG! The trouble is that,
according to the LRM, the result of the conversion is implementation-
dependent if the fractional part of the float quantity lies just 
between the two integers (that is, = 0.5).  (LRM sect. 4.6)

To see the hidden nastiness here, suppose Y has an integral value.
If Y = 10.0 (say), then Trunc(Y) returns 10 on compilers where the 0.5 
case rounds _up_, and returns 9 on compilers where the 0.5 case rounds 
_down_. Oops! A simple bit of code falls right into a portability trap.

The only solution seems to lie in forcing an integer division, since integer
division indeed truncates, portably. So our new function is

  FUNCTION Trunc (X: Float) RETURN Integer IS
  BEGIN
    RETURN Integer(2.0 * X) / 2;
  END Trunc;

which, I suppose, can be optimized efficiently, but seems like a
kludgy way to do a simple truncation. I think this is an example of
an operation that's much more easily done as an intrinsic than as
a programmer-defined function.

Any thoughts in net-land?

By the way - I am all the more red-faced because I let the dangerous
Trunc monster loose into a textbook. 'Course, I haven't run into 
another Ada text that even discusses the question...

Always in search of the portable Ada program, I remain

Mike

-------------------------------------------------------------------------------
Michael B. Feldman
Visiting Professor 1991-92               Professor
Dept. of Comp. Sci. and Engrg.           Dept. of Elect. Engrg. and Comp. Sci.
University of Washington FR-35           The George Washington University
Seattle, WA 98105                        Washington, DC 20052

mfeldman@cs.washington.edu               mfeldman@seas.gwu.edu
(206) 632-3794 (voice)                   (202) 994-5253 (voice)
(206) 543-2969 (fax)                     (202) 994-5296 (fax)
-------------------------------------------------------------------------------

Re: Red-faced professor gets bitten in search for portability

[Mike Murphy]

In article <1991Nov15.185959.5002@milton.u.washington.edu> mfeldman@milton.u.wa
shington.edu (Michael Feldman) writes:
>Asked by students whether Ada has an equivalent of the Pascal "Trunc" 
>operation, which just returns the integer part of its argument, 
>I put together a function based on the old trick we used
>in the Fortran days.
>
>  FUNCTION Trunc (X: Float) RETURN Integer IS
>  BEGIN
>    RETURN Integer(X - 0.5);
>  END Trunc;
>
>Since conversion of Float to Integer is, in Ada, a rounding operation,
>this looks like a good solution, right? WRONG! The trouble is that,
>according to the LRM, the result of the conversion is implementation-
>dependent if the fractional part of the float quantity lies just 
>between the two integers (that is, = 0.5).  (LRM sect. 4.6)
>
>To see the hidden nastiness here, suppose Y has an integral value.
>If Y = 10.0 (say), then Trunc(Y) returns 10 on compilers where the 0.5 
>case rounds _up_, and returns 9 on compilers where the 0.5 case rounds 
>_down_. Oops! A simple bit of code falls right into a portability trap.
>
>The only solution seems to lie in forcing an integer division, since integer
>division indeed truncates, portably. So our new function is
>
>  FUNCTION Trunc (X: Float) RETURN Integer IS
>  BEGIN
>    RETURN Integer(2.0 * X) / 2;
>  END Trunc;
>
>which, I suppose, can be optimized efficiently, but seems like a
>kludgy way to do a simple truncation. I think this is an example of
>an operation that's much more easily done as an intrinsic than as
>a programmer-defined function.
>
>Any thoughts in net-land?

You are right about your first version not being portable, but
your second version seems a bit convulated to me.  Why not just
check whether you indeed rounded down after doing the integer conversion?
For example:
	FUNCTION Trunc (X: Float) RETURN Integer IS
		itrunc : integer = integer(x - 0.5);
	BEGIN
		if x - float(itrunc) >= 1.0 then
			-- rounded down, so add back a 1
			return itrunc+1;
		else
			return itrunc;
		end if;
	END Trunc;

--mike
p.s. btw, some machines, e.g. MIPS, can either round down, or round up,
or round to the even number (the default setting).

Re: Red-faced professor gets bitten in search for portability

[Michael Feldman]

In article <11919@spim.mips.COM> murphy@mips.com (Mike Murphy) writes:
>
>You are right about your first version not being portable, but
>your second version seems a bit convulated to me.  Why not just
>check whether you indeed rounded down after doing the integer conversion?
>For example:
>	FUNCTION Trunc (X: Float) RETURN Integer IS
>		itrunc : integer = integer(x - 0.5);
>	BEGIN
>		if x - float(itrunc) >= 1.0 then
>			-- rounded down, so add back a 1
>			return itrunc+1;
>		else
>			return itrunc;
>		end if;
>	END Trunc;
Sure - this seems equivalent. I kinda like the halving and doubling,
because these are pretty fast to implement (as shifts) given a
good optimizer. Your method has the advantage, I suppose, of being immune
to a possible overflow from the doubling.

>p.s. btw, some machines, e.g. MIPS, can either round down, or round up,
>or round to the even number (the default setting).
Right. My argument about Trunc is that it's such a common (and easy to
implement) operation that it's better done as an intrinsic, so that
the compiler writer can use whatever hardware trickery is available.

Mike

Re: Red-faced professor gets bitten in search for portability

[Norman H. Cohen]

In article <1991Nov15.185959.5002@milton.u.washington.edu>
mfeldman@milton.u.washington.edu (Michael Feldman) describes the
difficulty in implementing a portable Trunc function in Ada.  He shows
a solution involing integer division and writes that this solution

>                                                ... seems like a
>kludgy way to do a simple truncation. I think this is an example of
>an operation that's much more easily done as an intrinsic than as
>a programmer-defined function.
>
>Any thoughts in net-land?

This is one of the primitive functions that the Numerics Rapporteur
Group of ISO-IEC/JTC1/SC22/WG9 is now working on, with the goal of
establishing an ISO standard Generic Primitive Functions package.
(Do not confuse this with the Generic Elementary Functions package,
which has, if I recall correctly, already been submitted to SC22 as a
draft standard.  The primitive functions are operations like truncation,
determining the exponent of a floating-point number, determining the
mantissa of a floating-point number, and so forth.  They are most
straightforwardly implemented by dealing directly with the target
machine's floating-point representation, so compilers are likely to
treat these functions as intrinsics, replacing calls to them with
inline code.  The elementary functions are higher-level operations like
the trig functions, logarithms, and square root, which can be implemented
portably in terms of the elementary functions.)

 > By the way - I am all the more red-faced because I let the dangerous
 > Trunc monster loose into a textbook. 'Course, I haven't run into
 > another Ada text that even discusses the question...

Don't lose sleep over it! I presume (pray) that no safety-critical
application will incorporate an algorithm from any textbook without
further careful scrutiny.  Compared to the flub in Dr. Ruth's book
(incorrectly describing the time of the month at which a woman can
conceive as the time at which the risk of pregnancy is LOWEST), this is
no big deal.

Norman H. Cohen

Re: Red-faced professor gets bitten in search for portability

[csus.edu!wupost!cs.utexas.edu!sun-barr!cronkite.Central.Sun.COM!newstop!s]

In <1991Nov16.000120.25606@milton.u.washington.edu> mfeldman@milton.u.washingto
n.edu (Michael Feldman) writes:

> In article <11919@spim.mips.COM> murphy@mips.com (Mike Murphy) writes:
> >
> >You are right about your first version not being portable, but
> >your second version seems a bit convulated to me.  Why not just
> >check whether you indeed rounded down after doing the integer conversion?
> >For example:
> >	FUNCTION Trunc (X: Float) RETURN Integer IS
> >		itrunc : integer = integer(x - 0.5);
> >	BEGIN
> >		if x - float(itrunc) >= 1.0 then
> >			-- rounded down, so add back a 1
> >			return itrunc+1;
> >		else
> >			return itrunc;
> >		end if;
> >	END Trunc;
> Sure - this seems equivalent. I kinda like the halving and doubling,
> because these are pretty fast to implement (as shifts) given a
> good optimizer. Your method has the advantage, I suppose, of being immune
> to a possible overflow from the doubling.

This code works (after fixing the syntax error).  Yours doesn't, Mike.

Your

  FUNCTION Trunc (X: Float) RETURN Integer IS
  BEGIN
    RETURN Integer(2.0 * X) / 2;
  END Trunc;

still takes the integer before the division.  In fact, as you've now doubled th
e
number, the horrible rounding beast is twice as active (i.e. 0.25 and 0.75
both fall victim).

I think Mike Murphy's code looks awkward because it's sort of upside down.  The
 
problem can be restated as: "if integer () rounds up, return one less than the
value of integer (), otherwise return the value of integer ()".

This translates to:

  function trunc (number: float) return integer is
  begin
    if float (integer (number)) > number then  -- does it round up?
      return integer (number) - 1;             -- yes - stop it!
    else
      return integer (number);                 -- no - ok as it is.
    end if;
  end trunc;

(I'm assuming that my compiler is clever enough to eliminate the redundant conv
ersions to integer.)

Dave

>From here down - my program which checks out all four trunc functions:

--------------------------------
with integer_io; use integer_io;
with text_io; use text_io;

procedure x is

package floating_io is new float_io (float); use floating_io;

r: float;
check: integer;
error,
first: boolean;

  FUNCTION Trunc_mf1 (X: Float) RETURN Integer IS
  BEGIN
    RETURN Integer(X - 0.5);
  END Trunc_mf1;

  FUNCTION Trunc_mf2 (X: Float) RETURN Integer IS
  BEGIN
    RETURN Integer(2.0 * X) / 2;
  END Trunc_mf2;

  FUNCTION Trunc_mm1 (X: Float) RETURN Integer IS
          itrunc : integer := integer(x - 0.5);
  BEGIN
    if x - float(itrunc) >= 1.0 then
      return itrunc+1; -- rounded down, so add back a 1
    else
      return itrunc;
    end if;
  END Trunc_mm1;
  
  function trunc_ds1 (number: float) return integer is
  begin
    if float (integer (number)) > number then  -- does it round up?
      return integer (number) - 1;             -- yes - stop it!
    else
      return integer (number);                 -- no - ok as it is.
    end if;
  end trunc_ds1;

  procedure try (message: string; value: integer) is
  begin
    put (" " & message & " => ");
    put (value, 2);
    if first then
      check := value;
      first := false;
    elsif check /= value then
      error := true;
    end if;
  end try;

begin -- x

  r := -5.0;
  loop
    put (r, 3, 2, 0);
    first := true;
    error := false;
    try ("mf1", trunc_mf1 (r));
    try ("mf2", trunc_mf2 (r));
    try ("mm1", trunc_mm1 (r));
    try ("ds1", trunc_ds1 (r));
    if error then
      put ("     ******");
    end if;
    new_line;
    r := r + 0.25;
    exit when r > 5.0;
  end loop;
  
end;
----------------------------

-- 
David Smart, Computer Sciences of Australia.
Net: daves@assip.csasyd.oz.au

Re: Red-faced professor gets bitten in search for portability

[micro-heart-of-gold.mit.edu!wupost!sdd.hp.com!uakari.primate.wisc.edu!use]

In article <9111191534.AA20295@ajpo.sei.cmu.edu> ncohen@WATSON.IBM.COM ("Norman
 H. Cohen") writes:
 [ stuff deleted ]
>
>This is one of the primitive functions that the Numerics Rapporteur
>Group of ISO-IEC/JTC1/SC22/WG9 is now working on, with the goal of
>establishing an ISO standard Generic Primitive Functions package.
>(Do not confuse this with the Generic Elementary Functions package,
>which has, if I recall correctly, already been submitted to SC22 as a
>draft standard.  The primitive functions are operations like truncation,
>determining the exponent of a floating-point number, determining the
>mantissa of a floating-point number, and so forth.  They are most
>straightforwardly implemented by dealing directly with the target
>machine's floating-point representation, so compilers are likely to
>treat these functions as intrinsics, replacing calls to them with
>inline code.  The elementary functions are higher-level operations like
>the trig functions, logarithms, and square root, which can be implemented
>portably in terms of the elementary functions.)

This is good. I think it's overdue and hope it happens soon.
>
> > By the way - I am all the more red-faced because I let the dangerous
> > Trunc monster loose into a textbook. 'Course, I haven't run into
> > another Ada text that even discusses the question...
>
>Don't lose sleep over it! I presume (pray) that no safety-critical
>application will incorporate an algorithm from any textbook without
>further careful scrutiny.  Compared to the flub in Dr. Ruth's book
>(incorrectly describing the time of the month at which a woman can
>conceive as the time at which the risk of pregnancy is LOWEST), this is
>no big deal.

Gosh! You mean safety-critical folks don't trust textbooks :-)

I'm red-faced but not insomniac about it. As a fellow author,
you'll understand the "Aarrgghh!" aspect of it. Dr. Ruth's bug was
clearly more serious.

Mike

Re: Red-faced professor gets bitten in search for portability

[micro-heart-of-gold.mit.edu!wupost!zaphod.mps.ohio-state.edu!uakari.prima]

In article <daves.690605701@condor> daves@assip.csasyd.oz.au writes:
 [ stuff deleted ]
>
>  FUNCTION Trunc (X: Float) RETURN Integer IS
>  BEGIN
>    RETURN Integer(2.0 * X) / 2;
>  END Trunc;
>
>still takes the integer before the division.  In fact, as you've now doubled t
he
>number, the horrible rounding beast is twice as active (i.e. 0.25 and 0.75
>both fall victim).

Indeed. Oops.
>
>I think Mike Murphy's code looks awkward because it's sort of upside down.  Th
e 
>problem can be restated as: "if integer () rounds up, return one less than the
>value of integer (), otherwise return the value of integer ()".
>
>This translates to:
>
>  function trunc (number: float) return integer is
>  begin
>    if float (integer (number)) > number then  -- does it round up?
>      return integer (number) - 1;             -- yes - stop it!
>    else
>      return integer (number);                 -- no - ok as it is.
>    end if;
>  end trunc;
>
>(I'm assuming that my compiler is clever enough to eliminate the redundant con
versions to integer.)
>

Better and simpler. Portable, too, I think.

Mike Feldman

Re: Red-faced professor gets bitten in search for portability

[psinntp!vitro.com!v7.vitro.com!vaxs09]

In article <11919@spim.mips.COM>, murphy@mips.com (Mike Murphy) writes:
> Why not just
> check whether you indeed rounded down after doing the integer conversion?
> For example:
> 	FUNCTION Trunc (X: Float) RETURN Integer IS
> 		itrunc : integer = integer(x - 0.5);
> 	BEGIN
> 		if x - float(itrunc) >= 1.0 then
> 			-- rounded down, so add back a 1
> 			return itrunc+1;
> 		else
> 			return itrunc;
> 		end if;
> 	END Trunc;
>
I think this version focuses more clearly on the original problem.

	FUNCTION Trunc ( X : Float ) RETURN Integer IS
		i : integer := integer ( x );
		if x < float(i) then
			return i - 1;
		else
			return i;
		endif;
	END Trunc;

Re: Red-faced professor gets bitten in search for portability

[The Sunset Kid]

In article <1991Nov22.090650.28@v7.vitro.com> vaxs09@v7.vitro.com writes:
>
>	FUNCTION Trunc ( X : Float ) RETURN Integer IS
>		i : integer := integer ( x );
>		if x < float(i) then
>			return i - 1;
>		else
>			return i;
>		endif;
>	END Trunc;

This version will not work for negative numbers.  If X is less than zero,
then a proper truncation would leave X < Float (I), subsequently subtracting 1,
and returning -2 for a parameter of -1.25.

Of course, I rarely use negative numbers, so change all your INTEGERs
to NATURALs and you'll have a pretty efficient truncator for non-negative
integers.  By the way, putting ABS calls in this thing would turn it into a
real mess... anyone else got an easier way of handling negatives?

C.A. Piepenbring

Re: Red-faced professor gets bitten in search for portability

[Jonathan Parker]

In article <1991Nov22.090650.28@v7.vitro.com> vaxs09@v7.vitro.com writes:
> In article <11919@spim.mips.COM>, murphy@mips.com (Mike Murphy) writes:
> > Why not just
> > check whether you indeed rounded down after doing the integer conversion?
> > For example:
> > 	FUNCTION Trunc (X: Float) RETURN Integer IS
> > 		itrunc : integer = integer(x - 0.5);
> > 	BEGIN
> > 		if x - float(itrunc) >= 1.0 then
> > 			-- rounded down, so add back a 1
> > 			return itrunc+1;
> > 		else
> > 			return itrunc;
> > 		end if;
> > 	END Trunc;
> >
> I think this version focuses more clearly on the original problem.
> 
> 	FUNCTION Trunc ( X : Float ) RETURN Integer IS
> 		i : integer := integer ( x );
> 		if x < float(i) then
> 			return i - 1;
> 		else
> 			return i;
> 		endif;
> 	END Trunc;
  There's still more to it, if I understand the Pascal Trunc function
correctly.  Trunc doesn't round, it truncates.  For example, Trunc(-3.999)
is -3, not -4.  Below is an excerpt from W. A. Whitaker's implementation
of the elementary math functions in Cody and Waite's famous book on the
subject.  I took this unmodified from files in the Ada repository on
SIMTEL-20. Floating is a generic Float type.   -- Jonathan  

    function Truncate (X : Floating) return Floating is
    --  Optimum code depends on how the system rounds at exact halves
        Z : Floating := Floating (Integer (X));
    begin
        if X = Z then
            return Z;
        else
            if X > Zero then
        --  For instance, you can get into trouble when 1.0E-20 is overwhelmed
        --  by HALF and X - HALF is -HALF and that is "truncated" to -1
        --  so a simple floating (Integer (X - Half)) does not work in all  
cases
                if X < Z then
                    return Z - One;
                else
                    return Z;
                end if;
            elsif X = Zero then
                return Zero;
            else
                if X < Z then
                    return Z;
                else
                    return Z + One;
                end if;
            end if;
        end if;
    end Truncate;

Re: Red-faced professor gets bitten in search for portability

[Dik T. Winter]

What is ignored in this thread is that floating-point arithmetic is a bit
indeterminate.  I will point to the problems in this version:

In article <1991Nov22.090650.28@v7.vitro.com> vaxs09@v7.vitro.com writes:
 > 	FUNCTION Trunc ( X : Float ) RETURN Integer IS
 > 		i : integer := integer ( x );
 > 		if x < float(i) then
What if x is the largest machine number smaller than i?  In that case it is
most probable that x is not a model number, so it is in the middle of a
model interval.  We find that the comparison is allowed to deliver both
true and false in this case.  On some machines you will indeed find false!
Even the comparison "x > float(i)" is allowed to return true (but I do not
know machines where that will happen).  The LRM clearly states what happens
when you compare numbers in different model intervals, or when you compare
different model numbers.  It also says that comparison of numbers in the
same model interval is indeterminate.

 > 			return i - 1;
 > 		else
 > 			return i;
 > 		endif;
 > 	END Trunc;


--
dik t. winter, cwi, kruislaan 413, 1098 sj  amsterdam, nederland
dik@cwi.nl (from bitnet also probably: dik@mcvax)

Re: Red-faced professor gets bitten in search for portability

[Dik T. Winter]

In article <1991Nov22.152426.25797@verdix.com> chrisp@verdix.com (The Sunset Ki
d) writes:
 > In article <1991Nov22.090650.28@v7.vitro.com> vaxs09@v7.vitro.com writes:
 > >
 > >	FUNCTION Trunc ( X : Float ) RETURN Integer IS
 > >		i : integer := integer ( x );
 > >		if x < float(i) then
 > >			return i - 1;
 > >		else
 > >			return i;
 > >		endif;
 > >	END Trunc;
 > 
 > This version will not work for negative numbers.
Actually it does not work at all.  I failed to notice in my followup, but
the two return statements must of course be interchanged.
--
dik t. winter, cwi, kruislaan 413, 1098 sj  amsterdam, nederland
dik@cwi.nl (from bitnet also probably: dik@mcvax)

Re: Red-faced professor gets bitten in search for portability

[csus.edu!wupost!spool.mu.edu!munnari.oz.au!metro!cluster!swift!sunaus!ass]

To follow up my own post, I've been reminded that my algorithm - and Mike 
Murphy's and both of Mike Feldman's - don't reproduce the Pascal trunc action 
for negative numbers.  E.g. -3.75 becomes -4 rather than -3.  In my own case, o
f
course, I was simply reproducing the action of Mike Murphy's algorithm :-).

There are tricky ways of getting the result right for negative values.  In 
Pascal, I'd probably use the abs () and sign () functions.  However, KISS:

  function trunc (number: float) return integer is
  begin
    if number >= 0.0 then                        -- positive cases
      if float (integer (number)) > number then    -- does it round up?
        return integer (number) - 1;               -- yes - stop it!
      else
        return integer (number);                   -- no - ok as it is.
      end if;
    else                                         -- negative cases
      if float (integer (number)) < number then    -- does it round down?
        return integer (number) + 1;               -- yes - stop it!
      else
        return integer (number);                   -- no - ok as it is.
      end if;
    end if;
  end trunc;

Dave
-- 
David Smart, Computer Sciences of Australia.
Net: daves@assip.csasyd.oz.au