To handle integer overflow in a multiplication, you could do the multiplication in floating point arithmetic, check if it exceeds huge(i), and return a sentinel value if it does, as done in the code below. The user should test if the result equals the sentinel. Is that a good approach? A related 2017 thread from the Intel Fortran forum was How best to handle integer overflow situations?
module product_mod
implicit none
private
public :: prod, dp
integer, parameter :: dp = kind(1.0d0)
contains
elemental function prod(i, j, bad) result(ij)
! return the product of i and j, unless it overflows, in which case return bad
integer, intent(in) :: i, j, bad
integer :: ij
real(kind=dp) :: xij
xij = real(i, kind=dp) * real(j, kind=dp)
if (abs(xij) > huge(i)) then
ij = bad
else
ij = i*j
end if
end function prod
end module product_mod
!
program xproduct
use product_mod, only: prod, dp
implicit none
integer :: i, isq, imin
imin = int(sqrt(real(huge(i), kind=dp)))
print*,"huge(i) =", huge(i)
print "(/,*(a20))", "i", "prod(i, i, -1)", "i*i"
do i=imin, imin+2
isq = i*i
print "(*(i20))", i, prod(i, i, -1), isq
end do
end program xproduct
output:
huge(i) = 2147483647
i prod(i, i, -1) i*i
46340 2147395600 2147395600
46341 -1 -2147479015
46342 -1 -2147386332
Nice question I would start by asking for an integer kind that can store large numbers, with selected_int_kind(). Then, we could do the multiplication and check if the result is finite with ieee_is_finite(). Now, what if the largest kind available was not enough? Well, then I would probably start translating Java’s BigInteger in Fortran, or look for a library that can handle arbitrary precision arithmetic.
I think there are several approaches that can be taken, depending on how much effort the programmer is willing to expend. If the default integer kind is int32 and if int64 is available, then as mentioned above the result can be computed and tested with int64, and then converted to int32 if successful. If only one integer kind is available, or if the calculations are being done in the longest kind already, then the following approach can be used.
integer(ik) :: i, j, ij, bigj
bigj = huge(ij) / abs(i)
if ( abs(j) > bigj ) then
istat = overflow_code
else
ij = i * j
istat = normal_code
endif
This obviously costs an integer division in addition to the multiplication, so it is expensive. I also would recommend using a subroutine interface rather than a function interface for this operation so that the istat value can be tested in a straightforward way in the calling program. There is also the potential problem of false negative when the product is -huge(ij)-1. That is a representable number, but it would be falsely flagged as an overflow in the above test (it does overflow the fortran model number range). Those special cases could also be tested, but with even more effort.
We have used procedures essentially identical to what you describe but we did use them as functions but they stopped the program if overflow was detected so there was no flag being returned.
We did something different for the largest kind on the Cray because (for reasons I do not remember at the moment) the cost of the division was onerous but can not access the archives to see if notes remain on what the alternative was.
But as stated, just considering commonly used kinds, our favorite little prime huge(0) has the condition huge(0_int32)**2 < huge(0_int64) , so all products of two 32-bit values are safe.
The alternative to division involved checking if both values had an absolute value < sqrt(huge(0)) and just proceeding and then went ahead and did the multiplication and did some checks on the sign bit I think, but
I think it was only safe for a particular programming environment.
@FortranFan comments in that thread that when multiplying two large enough integer numbers a wrap around the integer(4) representable number line can occur. My way to go is similar to what is proposed in the post and OP: since all integer(4) representable numbers can be represented as real(8) variables exactly, I keep a double precision “control” variable and compare the results between integer and real multiplications. For example, this program calculates 15!=1307674368000, which wraps around the number line several times,
program test_integer_overflow
implicit none
integer, parameter :: dp = 8
integer :: i
integer :: a(15) = [(i, i = 1, size(a))]
real(dp) :: tst = 1.0_dp
do i = 1, size(a)
tst = tst*real(a(i), dp)
enddo
print*, "Multiplication with reals: ", tst, "."
print*, "Integer conversion: ", int(tst), "."
print*, "Multiplication with integers, ", product(a), "."
print*, "Absolute difference:", abs(tst - real(product(a), dp)), "."
print*, "Epsilon (double precision):", epsilon(1.0_dp), "."
print*, "Overflow?", (abs(tst - real(product(a), dp))>epsilon(1.0_dp)), "."
end program test_integer_overflow
Output:
Multiplication with reals: 1307674368000.0000 .
Integer conversion: -2147483648 .
Multiplication with integers, 2004310016 .
Absolute difference: 1305670057984.0000 .
Epsilon (double precision): 2.2204460492503131E-016 .
Overflow? T .
There are some situations where the programmer cannot recover from the overflow and that is the only option. In other cases, the return status code allows the programmer to use alternative algorithms, or to modify the algorithm somehow, and still obtain the final solution despite the initial integer overflow.
One way to handle both situations is to let the istat dummy variable have the optional attribute. If it is present, then the overflow_code value is returned and execution continues. If it is not present, then the program can abort at that time (e.g. with error stop). Thus both situations are handled in a straightforward way with a single subroutine.
Regarding whether to use int64 or real64 to compute the results in extended precision, modern fortran (since f2008 I think) has required an extended integer kind with at least 18 decimal digits. All compilers I use provide int64 to satisfy this requirement, although a longer integer type would also be allowed. Thus the integer product can be computed exactly with this extended integer kind and the programmer can depend on its support by the compiler in a portable way. Before f2008, only the default integer kind was required to exist, typically int32, so relying on compiler support of any extended integer kind was risky. In contrast, double precision was required going back to f77 (at least in the full implementation), so within these earlier standards that was the most reliable and portable way to compute the extended precision intermediate; even if it was typically only 53-bits of mantissa, the overflow of 31 bits could still be detected reliably.
I would start by asking for an integer kind that can store large numbers, with