Re: Inline round for IA64

[Canqun Yang <canqun at nudt dot edu dot cn>]

Geert Bosch <bosch@adacore.com>:

> As far as I can seem from this patch, it rounds 
incorrectly.
> This is a problem with the library version as well, 
I believe.
>
> The issue is that one cannot round a positive float 
to int
> by adding 0.5 and truncating. (Same issues with 
negative values
> and subtracting 0.5, of course). This gives an error 
for the
> predecessor of 0.5. The between Pred (0.5) and 0.5 
is half that of
> pred (1.0) and 1.0. So the value of Pred (0.5) + 0.5 
lays exactly
> halfway Pred (1.0) and 1.0. The CPU rounds this 
halfway value to
> even, or 1.0 in this case.
>
> So try 
rounding .499999999999999944488848768742172978818416595
4589843750
> using IEEE double on non-x86 platform, and you'll 
see it gets rounded
> to 1.0.

Do you mean the correct value should be 0.0 ?

> A similar  problem exists with large odd integers 
between 2^52+1 and
> 2^53-1,
> where adding 0.5 results in a value exactly halfway 
two integers,
> rounding up to the nearest even integer. So, for 
IEEE double,
> 4503599627370497 would round to 4503599627370498.

Do you mean 4503599627370498 is a wrong result?

>
> These issues can be fixed by not adding/subtracting 
0.5, but Pred (0.5).
> As shown above, this rounds to 1.0 correctly for 
0.5. For larger values
> halfway two integers, the gap with the next higher 
representable number
> will
> only decrease so the result will always be rounded 
up to the next higher
> integer. For this technique to work, however, it is 
necessary that the
> addition will be rounded to the target precision 
according to IEEE
> round-to-even semantics. On platforms such as x86, 
where GCC implicitly
> widens intermediate results for IEEE double, the 
rounding to integer
> should be performed entirely in long double mode, 
using the long double
> predecessor of 0.5.
>
> See ada/trans.c around line 5340 for an example of 
how Ada does this.
>
>    -Geert
>
> On Apr 7, 2005, at 05:38, Canqun Yang wrote:
> > Gfortran translates the Fortran 95 intrinsic DNINT 
to
> > round operation with double precision type argument
> > and return value. Inline round operation will 
speed up
> > the SPEC CFP2000 benchmark 189.lucas which contains
> > function calls of intrinsic DNINT from 706 (SPEC
> > ratio) to 783 on IA64 1GHz system.
> >
> > I have implemented the double precison version of
> > inline round. If it is worth doing, I can go on to
> > finish the other precision mode versions.
>
> 

I attached an example for intrinsic DNINT with its 
output. Would you please check it, and tell me whether 
the result is correct.

Canqun Yang
Creative Compiler Research Group.
National University of Defense Technology, China.

! Test case for inline round
subroutine dnint_ex (a, b, n)
   real*8 a(n), b(n)
   integer n
   do i = 1, n
      b(i) = dnint (a(i))
   enddo
end

program round_test
   real*8 a(2), b(2)

   a(:) = (/.4999999999999999444888487687421729788184165954589843750_8,&
	    4503599627370497.0_8/)
   call dnint_ex (a, b, 2)
   write (*,*) b
end
  
The output is:
 
   0.00000000000000       4.503599627370497E+015