The ran_init subrroutine of the Numerical Recipes contains this lines:
INTEGER(K4B) :: new,j,hgt
...
hgt=hg
...
if (hgt+1 /= hgng) call nrerror('ran_init: arith assump 3 fails')
Where K4B, hgng and hg are globally declared in the module via:
INTEGER, PARAMETER :: K4B=selected_int_kind(9)
INTEGER(K4B), PARAMETER :: hg=huge(1_K4B), hgm=-hg, hgng=hgm-1
The problem is that in one particular computer (but not in the others) I get the error ran_init: arith assump 3 fails. The only thing that I get from the documentation about this error is:
Bit of dirty laundry here! We are testing whether the most positive integer hg wraps around to the most negative integer hgng when 1 is added to it. We can’t just write hg+1 , since some compilers will evaluate this at compile time and return an overflow error message. If your compiler sees through the charade of the temporary variable hgt , you’ll have to find another way to trick it!
How can I trick it?
the compiler during the optimization step can easily prove that hgt+1 and hgng cannot be the same and therefore it optimizes the whole condition to .false.. Adding two positive integers cannot create a negative integer.
You might be able to avoid it by using smaller level of optimization or by some masquerade with temporary variables which would be more elaborate then what they did in the old Numerical Recipes. You may try to use the -fwrapv or -fno-strict-overflow flag. But it is highly unsure and will work only in a given version of the compiler with certain compiler flags.
The easiest "fix" is probably just to delete the offending check. It is quite likely to work in many circumstances (and fail terribly in others).
What Numerical Recipes do is against the standard of Fortran. They assume that if you add 1 to the greatest integer, you get the lowest one.
This is not allowed for Fortran integers and neither is it allowed for signed integers in C. It is only allowed for unsigned integers in C and Fortran doesn't have them.
So the compiler can safely assume that adding two positive integers can never result in a negative integer. See https://gcc.gnu.org/bugzilla/show_bug.cgi?id=30475 for thorough discussion of these optimizations.
One option is to rewrite the operations which may lead to overflow in C and use conversion to unsigned integers. I use this in the random number generator I use (based on http://www.cmiss.org/openCMISS/wiki/RandomNumberGenerationWithOpenMP):
#include <stdint.h>
#include <memory.h>
int32_t sum_and_overflow (int32_t a, int32_t b)
{
uint32_t au, bu, su;
int32_t s;
memcpy(&au, &a, sizeof(a));
memcpy(&bu, &b, sizeof(b));
su = au + bu;
memcpy(&s, &su, sizeof(s));
return s;
}
with Fortran interface
interface
function sum_and_overflow(a, b) result(res) bind(C, name="sum_and_overflow")
use, intrinsic :: iso_c_binding
integer(c_int32_t) :: res
integer(c_int32_t), value :: a, b
end function
end interface
and instead of
c = a + b
I call
c = sum_and_overflow(a, b)
So in your case
if (sum_and_overflow(hgt,1) /= hgng) ...
but not only there, you would have to at least one more place where is this assumption used in the generator. In the generator I use there is only one such line.
There are many other hacks which may lead for temporary success, but will fail later with some other compiler options. For example, the undefined behaviour sanitizations in GCC don't like overflow in signed integers in Fortran and C and will terminate your program if you do that.
This is why I try to show a solution which is more complicated, but is not working around the standard, but following it.