Work at SourceForge, help us to make it a better place! We have an immediate need for a Support Technician in our San Francisco or Denver office.

Close

#5 bug: casting to float where double precision should be used

open
nobody
None
5
2013-04-12
2013-04-12
Stefan Reiser
No

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
[...]

F2J produces:

// Line (1)

tuplim =
(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).

Regards
Stefan Reiser

Discussion