experchange > fortran

steve kargl (02-28-19, 05:41 PM)
john.collins.simcon wrote:

[..]
> If we can't split an expression into triples because the evaluation depends on the context, we need a more sophisticated pre-processor.
> Do we?
> Best wishes,


The expression is invalid Fortran for Fortran processors where the default integer kind
corresponds to a signed 32-bit integer. As it is invalid Fortran, a Fortran processor
can do anything with expression. This includes giving a result that may look
like what one was hoping to achieve.
gah4 (02-28-19, 09:47 PM)
On Thursday, February 28, 2019 at 7:41:37 AM UTC-8, steve kargl wrote:

(snip)

> The expression is invalid Fortran for Fortran processors where
> the default integer kind corresponds to a signed 32-bit integer.
> As it is invalid Fortran, a Fortran processor can do anything with
> expression. This includes giving a result that may look
> like what one was hoping to achieve.


Reminds me of a program years ago that accidentally worked when
an uninitialized variable was not zero. Later, on a system that
initialized variables to zero, it failed.

While the above is convenient, it could surprise someone later
porting to another system.

Traditionally, Fortran systems were good at detecting floating
point overflow, but ignoring fixed point overflow. Some programs
(especially random number generators) depend on that.

It would be nice to be able to turn this feature off when debugging,
so as to detect the bug. (Even if the overflow isn't detected,
the very wrong result might be.)
steve kargl (03-01-19, 07:44 PM)
steve kargl wrote:

> john.collins.simcon wrote:
> The expression is invalid Fortran for Fortran processors where the default integer kind
> corresponds to a signed 32-bit integer. As it is invalid Fortran, a Fortran processor
> can do anything with expression. This includes giving a result that may look
> like what one was hoping to achieve.


Consider,

% cat a.f90
i = 10**23
x = 6.02 * 10**23
print *, i, x
end

My version of gfortran gives

% gfcx -c a.f90
a.f90:1:9:

1 | i = 10**23
| 1
Error: Result of exponentiation at (1) exceeds range of INTEGER(4)
a.f90:2:16:

2 | x = 6.02 * 10**23
| 1
Error: Result of exponentiation at (1) exceeds range of INTEGER(4)
robin.vowels (03-02-19, 09:30 AM)
On Thursday, February 28, 2019 at 11:45:23 PM UTC+11, john.coll...@gmail.com wrote:
> On Friday, 15 February 2019 00:30:00 UTC, spectrum wrote:
> I find Spectrum's result worrying. The rules of operator precedence imply that
> 6.02 * 10 ** 23
> will be evaluated as
> 6.02 * ( 10 ** 23 )
> The behaviour of gfortran implies that the evaluation of 10**23 depends on its context in an expression. Does the standard imply that?


No. Previous discussions here already gave that answer.

> This is important to us because we created a pre-processor to convert Fortran for a fine-grain parallel processor, which split all expressions into triples. Thus,
> r = 6.02 * 10 ** 23
> would be split into
> itemp_1 = 10 ** 23
> r = 6.02 ** itemp_1
> Where itemp_1 is a default integer (And would overflow)


Naturally.

> If we can't split an expression into triples because the evaluation depends on the context, we need a more sophisticated pre-processor.
> Do we?


No.
robin.vowels (03-02-19, 09:45 AM)
On Friday, March 1, 2019 at 12:50:56 AM UTC+11, Wolfgang Kilian wrote:
> On 28.02.2019 13:45, john.....@gmail.com wrote:
> According to the Fortran 2018 standard, 10.1.5.2.4:
> I think that gfortran takes advantage of that rule here,


The rule is meant to cover components of an expression
(such as a repeated sub-expression) that may be factored out
and evaluated prior to evaluating the expression, or which
may be evaluated in a different order (but paying attention to
operator precedence, of course).

The fact that exponentiation and multiplication are the only
operators here, and that exponentiation has higher priority,
means that 10**23 must be evaluated first.
The only way to evaluate that is with integer arithmetic.

And multiplication follows.
[..]
robin.vowels (03-02-19, 09:49 AM)
On Friday, March 1, 2019 at 6:47:04 AM UTC+11, ga...@u.washington.edu wrote:
[..]
> It would be nice to be able to turn this feature off when debugging,
> so as to detect the bug. (Even if the overflow isn't detected,
> the very wrong result might be.)


In PL/I, integer overflow is detetecd as part of the language.
Integer over overflow can be ignored when required, with a particular
language feature.
Dick Hendrickson (03-02-19, 07:23 PM)
On 3/2/19 1:30 AM, robin.vowels wrote:
> On Thursday, February 28, 2019 at 11:45:23 PM UTC+11, john.coll...@gmail.com wrote:
> No. Previous discussions here already gave that answer.
> Naturally.

I'm not sure. Think about
X*Y * I/J
versus
X*Y + I/J
I'm sure in the first one the division is a real division and in the
second the division is an integer division. Depends on how you split
things and whether or not you figure out intermediate types for the
triples.

Dick Hendrickson
gah4 (03-03-19, 12:45 AM)
On Saturday, March 2, 2019 at 9:24:02 AM UTC-8, Dick Hendrickson wrote:
> On 3/2/19 1:30 AM, robin.vowels wrote:
> I'm not sure. Think about
> X*Y * I/J
> versus
> X*Y + I/J
> I'm sure in the first one the division is a real division and in the
> second the division is an integer division. Depends on how you split
> things and whether or not you figure out intermediate types for the
> triples.


In the first, it is (X*Y*I)/J so real division, or at least
mathematically, if not numerically, equal.

Reminds me of a story from a book on how not to write programs:

DO 1 I=1,N
DO 1 J=1,N
1 A(I,J) = (I/J)*(J/I)

(Fortran 66 days.)

Note that this won't work for C programmers, or Fortran
programmers with 0 origin arrays. But most matrix arithmetic
is done in mathematics with 1 origin arrays.
john.collins.simcon (03-04-19, 01:33 PM)
> I'm not sure. Think about
> X*Y * I/J
> versus
> X*Y + I/J
> I'm sure in the first one the division is a real division and in the
> second the division is an integer division. Depends on how you split
> things and whether or not you figure out intermediate types for the
> triples.
> Dick Hendrickson


In the first case I would expect left to right evaluation, so:
X*Y * I/J
would evaluate as
(X*Y*I) / J
and the division would be real.

Consider:

PROGRAM t
INTEGER :: i,j
REAL :: x,y,r1,r2,r3
DO
WRITE(*,'("i,j: ",$)')
READ(*,*)i,j
WRITE(*,'("x,y: ",$)')
READ(*,*)x,y
r1 = x*y*i/j
r2 = i/j*x*y
r3 = i/(j*x*y)
WRITE(*,'("x*y*i/j: ",E12.3)')r1
WRITE(*,'("i/j*x*y: ",E12.3)')r2
WRITE(*,'("i/(j*x*y): ",E12.3,/)')r3
ENDDO
END PROGRAM t

Under gfortran (GNU Fortran (Ubuntu 7.1.0-10ubuntu1~16.04.york0) 7.1.0) a typical behaviour is:

i,j: 4 2
x,y: 0.1 0.2
x*y*i/j: 0.400E-01
i/j*x*y: 0.400E-01
i/(j*x*y): 0.100E+03

i,j: 2 4
x,y: 0.1 0.2
x*y*i/j: 0.100E-01
i/j*x*y: 0.000E+00
i/(j*x*y): 0.250E+02

i/j*x*y is evaluated as (i/j)*x*y, apparently with left to right ordering of operations at the same level of precedance. (i/j) is an integer division.
steve kargl (03-04-19, 05:07 PM)
john.collins.simcon wrote:

> In the first case I would expect left to right evaluation, so:


The Fortran standard does require left-to-right evaluation.
It requiers integriity of parentheses.

> X*Y * I/J
> would evaluate as
> (X*Y*I) / J
> and the division would be real.


That is one possible choice. A conforming processor can
evaluate the original expression as if it were written as
(x*y)*(i/j) or x*(y*i)/j or x * i * (y / j) or ...

[..]
> i/j*x*y: 0.000E+00
> i/(j*x*y): 0.250E+02
> i/j*x*y is evaluated as (i/j)*x*y, apparently with left to right ordering of operations at the same level of precedance. (i/j) is an integer division.


Yes, gfortran ueses left-to-right evaulation. It is not required. The gfortran frontend
leverages the middle and back end of GCC, which started life as a C compiler. Guess
what language requires left-to-right evaulation?
steve kargl (03-04-19, 05:41 PM)
steve kargl wrote:

> john.collins.simcon wrote:
> The Fortran standard does require left-to-right evaluation.
> It requiers integriity of parentheses.


Dang (never post at 7:00 am).

does *NOT* require.
Dick Hendrickson (03-04-19, 08:08 PM)
On 3/4/19 9:41 AM, steve kargl wrote:
> steve kargl wrote:
> Dang (never post at 7:00 am).
> does *NOT* require.

True about evaluation, but the syntax rules that define the
interpretation of an expression are left-to-right (except for **).

From 7.1.8.3 (F2003?)
The rules given in 7.2.1 specify the interpretation of a numeric
intrinsic operation. Once the interpretation has been established in
accordance with those rules, the processor may evaluate any
mathematically equivalent expression, provided that the integrity of
parentheses is not violated.

I think this means that in X*Y*I/J the divide is a real divide. In
fact, now that I'm actually reading the standard instead of going on
memory, note 7.18 expressly forbids transforming X*I/J into X*(I/J).

The [pre]processor can break up the expression in any way it likes, but
it must convert I and J into reals before using them.

Dick Hendrickson
john.collins.simcon (03-04-19, 08:41 PM)
> Dang (never post at 7:00 am).
> does *NOT* require.
> --
> steve


No problem - that's what I thought you meant.

The issue here is that the results of the expressions

r1 = x*y*i/j
r2 = i/j*x*y

or any multiplication - division chain where there is at least one real number and there are integers both sides of any division, are undefined by thestandard. A compiler may always choose to evaluate (i/j) first, causing an integer division, or not to evaluate (i/j) avoiding an integer division.

I will try the test program on other compilers, but even if a compiler behaves as though it uses left-to-right ordering on a simple test program it doesn't mean that it will always do so. I will also write a test for this situation and test our code library to see how often it occurs (But this willtake some time).

Best wishes,

John
gah4 (03-04-19, 09:08 PM)
On Monday, March 4, 2019 at 10:41:45 AM UTC-8, john.coll...@gmail.com wrote:

(snip)

> No problem - that's what I thought you meant.
> The issue here is that the results of the expressions
> r1 = x*y*i/j
> r2 = i/j*x*y


> or any multiplication - division chain where there is at least
> one real number and there are integers both sides of any division,
> are undefined by the standard.


> A compiler may always choose to evaluate (i/j) first, causing
> an integer division, or not to evaluate (i/j) avoiding an
> integer division.


No.

The result has to be mathematically equivalent to following
the precedence and association rules, which are left to right
in most cases.

The result doesn't have to be numerically equal, though.
(Especially floating point rounding might be different.)

Reordering might generate overflow, especially floating point
overflow, that a different order wouldn't generate.

But all the same integer division will happen, or not, as the
case may be.

I suspect that compilers won't reorder such that precision
would be reduced, though I am not sure that the standard
specifies that. As we know, especially from systems using
the x87 floating point processor with extra precision,
expressions can be evaluated with more precision.

Mostly this should matter when common subexpression elimination
is done, where the subexpressions are not the ones that would
happen with the obvious evaluation order.

X=B*C
Y=A*B*C*D

a compiler might notice that reorder could allow it to use
the already computed B*C. (Older compilers required common
subexpressions to be easier for the compiler to find. People
were encouraged to help the compiler with these.)

Since integer expressions don't have precision roundoff,
it is easier to rearrange for the exact result.

For example, with the usual twos complement overflow,
addition with intermediate (ignored) overflow gives the same
result for any order. That isn't true in floating point.
jfh (03-04-19, 11:00 PM)
On Saturday, March 2, 2019 at 6:45:00 AM UTC+13, steve kargl wrote:
[..]
> Error: Result of exponentiation at (1) exceeds range of INTEGER(4)
> --
> steve


I don't think anyone mentioned that if your compiler allows integers with 24 digits, and default integers may have 9 digits, both OK in gfortran, then you can get an integer value of Avogadro's number like this:

implicit none
integer,parameter:: bigi=selected_int_kind(24)
integer(bigi):: avo=602214076*10_bigi**15
print "(I0)",avo
end program

Similar Threads