fast and accurate in mixed mode operations
in [Fortran]
|
Prev: a simple overflow example
Next: Overload assignment operator for an user defined type array with ?vector scripts
From: Thomas Koenig on 12 Sep 2009 08:41 On 2009-09-12, Dan Nagle <dannagle(a)verizon.net> wrote: > Division by a constant, on modern compilers, it is almost always > replaced by multiplication by the reciprocal. gcc/gfortran doesn't do that, even with -ffast-math. Do you think it should? Here's an example, BTW: $ cat divide.f90 subroutine divide(a, b) real, dimension(2), intent(inout) :: a real, intent(in) :: b a(1) = a(1) / b a(2) = a(2) / b $ gfortran -O3 -ffast-math -march=native -mfpmath=sse -msse -S divide.f90 $ cat divide.s .file "divide.f90" .text .p2align 4,,15 ..globl divide_ .type divide_, @function divide_: pushl %ebp movl %esp, %ebp movl 8(%ebp), %eax movl 12(%ebp), %edx movss (%eax), %xmm0 movss (%edx), %xmm1 divss %xmm1, %xmm0 movss %xmm0, (%eax) movss 4(%eax), %xmm0 divss %xmm1, %xmm0 movss %xmm0, 4(%eax) leave ret
From: nmm1 on 12 Sep 2009 08:49 In article <slrnhan5od.4p3.tkoenig(a)meiner.onlinehome.de>, Thomas Koenig <tkoenig(a)netcologne.de> wrote: >On 2009-09-12, Dan Nagle <dannagle(a)verizon.net> wrote: > >> Division by a constant, on modern compilers, it is almost always >> replaced by multiplication by the reciprocal. > >gcc/gfortran doesn't do that, even with -ffast-math. > >Do you think it should? Yes, I do. Strongly. Inter alia, it's arguably a bug not to so so. -ffast-math sets -funsafe-math-optimizations, the latter sets -freciprocal-math, and the last is specified to do precisely this. Regards, Nick Maclaren.
From: Tim Prince on 12 Sep 2009 09:55 robin wrote: > "Tim Prince" <tprince(a)nospamcomputer.org> wrote in message news:7gsp7pF2qhnrlU1(a)mid.individual.net... > | Dan Nagle wrote: > | > Hello, > | > > | > On 2009-09-10 06:51:39 -0400, "robin" <robin_v(a)bigpond.com> said: > | > > | >> There's no requirement to do any conversion at compile time. > | >> x/2 can be done at execution time on certain computers > | >> not by division but by scaling (and at considerable saving in time). > | > > | > Masking out the exponent, shifting, subtracting one, > | > and merging the new exponent back into the original number > | > may well take longer than one multiplication > | > on modern hardware. > > | Not to mention checking for over/underflow. > > A test for underflow has to be performed regardless of whether > actual division or simple halving is performed. That observation > therefore is irrelevant. A divide operation takes care of possible underflow. Subtracting from the exponent field without checking against limits may wrap the result around past HUGE. Since formats such as infinities and sub-normals were introduced (over 25 years ago), more cases of failure for such shortcuts have been present. > > | The compiler I learned on > | generated code for /2. which jumped over the subtraction from the > | exponent in the case of a 0. operand but didn't take care of all corner > | cases. > > What machine was that? > > On the Honeywell 6000 series, x/2. would not fail, but x/4. would produce HUGE(x) if x==TINY(x). Of course, there were no TINY or HUGE; they gave up after a partial implementation of f77.
From: steve on 12 Sep 2009 11:57 On Sep 12, 5:41 am, Thomas Koenig <tkoe...(a)netcologne.de> wrote: > On 2009-09-12, Dan Nagle <danna...(a)verizon.net> wrote: > > > Division by a constant, on modern compilers, it is almost always > > replaced by multiplication by the reciprocal. > > gcc/gfortran doesn't do that, even with -ffast-math. Read what Dan wrote. "Division by a constant". REMOVE:kargl[79] cat jo.f90 subroutine divide(a) real, dimension(2), intent(inout) :: a a(1) = a(1) / 2 a(2) = a(2) / 2 fc4x -c -ffast-math -O2 -save-temps jo.f90 divide_: pushl %ebp movl %esp, %ebp movl 8(%ebp), %eax flds .LC0 flds (%eax) fmul %st(1), %st fstps (%eax) fmuls 4(%eax) fstps 4(%eax) leave ret Note, you can change the 2 to any floating point number you like. You still get to fmul. -- steve
From: glen herrmannsfeldt on 12 Sep 2009 12:49
Dan Nagle <dannagle(a)verizon.net> wrote: (snip, someone wrote) <> Multipliocation? whio said anything about multiplication? <> The operation is DIVISION here. < Division by a constant, on modern compilers, it is almost always < replaced by multiplication by the reciprocal. < I thought you knew that. Floating point multiplication by a constant might round differently than division. Many will consider that reasonable. In fixed point, you can multiply be a reciprocal and take the high half of the product. The result is properly truncated for many but not all divisors. As people depend on exact results for fixed point divide, compilers better get that one right. -- glen |