I think this is a bug in the F2J code generation:
As I recall, the fortran rule for "mixed mode arithmetic" is that a calculation is to be performed in double precision when the widest of the operants is double. I'm not sure if this also applies when this widest operant is the result of a call to a double precision function (?)
- if so, then the following output of F2J is incorrect:
(Input from Slatec, drd.f)
(1) TUPLIM = D1MACH(1)**(1.0E0/3.0E0)
(2) TUPLIM = (0.10D0*ERRTOL)**(1.0E0/3.0E0)/TUPLIM
UPLIM = TUPLIM**2.0D0
(3) IF (MAX(X,Y,Z).GT.UPLIM) THEN
IER = 3
// Line (1)
(double)(((float)Math.pow(D1mach.d1mach(1), ((1.0E0f / 3.0E0f)))));
// (Here D1MACH(1) yields Double.MIN_VALUE.) tumplim in line 1
// evaluates to about 1.70D-108 which becomes 0.0 when casted to float.
// Now tuplim in line (2) and later uplim evaluate to +Infinity:
tuplim = ((Math.pow(((0.10e0 * errtol)), ((1.0E0f / 3.0E0f)))) / tuplim);
// ... so the check of the arguments range in line (3) is rendered
// ineffective as it now reads "if (some value > +Infinity)". Later this
// lets DRD loop infinitely.
I'm aware that the true reason for the problem is due to the use "1.0E0/3.0E0" in the original fortran code where "1.0D0/3.0D0" should have been used.
Nevertheless the intermediate result should not been casted to float in line (1) (note that for line (2) the mixed mode rule has been applied correctly).