Bug (?) in Intel Math library

[WSinc]

This issue came up before regarding computations involving literals.
Apparently the same library is being used for variables as well.


integer(4) ip3
integer(8) x3
do ip3=1,1000
x3=ip3**4
print *,ip3,x3
if(x3.lt.0)pause
end do

This loop should go to completion, since an 8 byte integer should easily contain
1000**4 = 10^12. Its highest number is about 9.22 * 10^18 in fact.

However the loop quits at ip3=216, giving a NEGATIVE result for x3.
This is because whoever coded the library routine did not understand that the
result does not necessarily have the same number of bytes as the input
quantity.

In other words, the internal calculations should have more bytes than the input
quantity. This seems to me to be a very obvious thing they overlooked.

Of course, if x3 were the same size, you would get an integer overflow.
We don't see these problems with floating point, as far as I know.

[mecej4]

Your comments, however reasonable they may be from a mathematical point of view, are in conflict with the Fortran standard.

In an assignment statement, the expression has a type. This type is determined by the contents of the expression, completely independent of the type of the assignee on the left of the = sign.

If the expression is built up from a number of operands of numeric types and operators, there is a table that prescribes what the type of the expression is. For example, see Table 3.1 of Metcalf, Reid and Cohen, Fortran 95/2003 Simplified, OUP, 2006.

The statement

x3=ip3**4

tries to assign to x3, an 8-byte integer variable, the value of an expression of type default integer, i.e., a 4-byte integer. The evaluation of this expression causes overflow to occur before the necessary promotion of the result to an 8-byte integer.

Replacing the statement in question by

x3=ip3**4_8

makes the code work as you desired.

And, finally, the same qualitative behavior occurs even with floating point variables, as this short program will show you:

[fortran]program r8
real(4) p3
real(8) x3
integer ip3
do ip3=1,1000
    p3=ip3*1e9
    x3=p3**4
    print *,ip3,p3,x3
    if(x3.gt.huge(p3))stop 'Overflow'
end do
end program r8
[/fortran]

which prints out the following:

[bash]           1  1.0000000E+09  9.999999616903162E+035
           2  2.0000000E+09  1.599999938704506E+037
           3  3.0000000E+09  8.100000339362494E+037
           4  4.0000000E+09  2.559999901927210E+038
           5  5.0000000E+09 Infinity
Overflow
[/bash]

[WSinc]

Well, the problem seems to be how we approach this question.

When we are doing simulations of a complex physical phenomenon,
we usually approach things from "a mathematical point of view."
Otherwise, why bother?

Unfortunately, Intel does not catch integer overflows, so when the internal calculations
are not done with enough byte allocation, the answers can be completely wrong, and
we have no way of knowing it. Floating point is no problem, it DOES trap.

For example, this simple multiply would give a wrong answer:

integer(4) n1,n2
integer(8) n3
data n1,n2/100000, 120000/
n3=n1*n2

A good general rule is: If either the LEFT or the RIGHT side is 8 bytes,
all internal calculations should be also. There would be some loss of perfermance,
but at least its safe. That was the rule at JPL anyway.

I hope in the future, Intel will put an integer overflow trap in their CPU.