I’ve read (and cited above) this statement. Still, I consider it contradicts the previous. Left-to-right interpretation of adjacent ± or */ operators was in Fortran 77 (sect. 6.1.2.3-4) . Now if in the Modern Fortran the processor may choose any association, then the contents of the Note 10.9 has no sense. Because it states, e.g., that 2.1*3.4*4.9 has the same interpretation as (2.1*3.4)*4.9. And the latter form, having parentheses added, leave no choice to the processor.
If the standard committee decided to abandon this rule, it should have been clearly stated. Instead, we have contradicting clauses in the very same section of the standard.
Point the lecturer at this discussion!
And well done for digging
The left-to-right refers to how the expression is parsed (“interpretation”). The subsequent section (“evaluation”) states that it needs not honor the expression as parsed and can replace it.
1 Like
So, what is exactly the rule for the type of the iterator? Is this standard, or a coincidence that it works with three compilers?
program testit
use,intrinsic::iso_fortran_env,only:int8,int16,int32,int64
integer(kind=int64) I
do i=1,huge(0_int32)
! do stuff
enddo
write(*,*)i
Since everyone was just citing the standard, I was wondering what the type of the iterator was allowed to be and whether it is specified what type is used. I was actually expecting to see something saying it had to be of default integer kind but did not see it right away; thinking this might be non-standard and a side-effect of the compilers allowing other iterator types like reals as an extension.
Does not eliminate the same issue with huge(0_int64), but lets you iterator over all the int32 values (although infinity/huge(0_int32) is still a pretty big number 
)
PS: Found it, where it specifies M1, M2, M3 are the same type as the iterator, so that is standard, but could cause problems in the body of the loop if the type of the iterator matters and so on.
Why is this worrying or an operator error? He knew exactly how a do-loop is defined in the Fortran standard. If compilers would implement it exactly like that, do i=1,huge(i) wouldn’t be a problem. The loop would have to exit, because the iteration count is zero, despite that i overflows.
I suggested it is worrying to have an example “do i=1, huge(i)” .
For a loop to execute correctly as an indexed DO loop, i would = huge(i)+1 on exit; which is not possible.
I don’t know the context of a Fortran lecturer providing this example; perhaps an example of what to avoid ?
The example was the sum from i=1 to infinity over 1./real(i). He wanted to explain the difference between forward and backward summation. He then asked the students, how they would implement this in Fortran. So to be 100% precise, it was the suggestion from a student, but since he knew the loop definition with the iteration count he didn’t complain about it. Since all took place on a remote version of a whiteboard, there was neither compilation nor execution involved.
BTW: I assisted the lesson and tried to run the example. This is how I discovered this problem. Of course I corrected the example and later I even introduced Kahan to the students.
In my opinion it doesn’t matter if i overflows at the end of a loop, because the loop exits before it could use the last value. As long as you don’t expect i to be useful after the loop, this shouldn’t break your program.
If so, that being the harmonic series, which is known to be divergent, the students are being invited to explore different ways of computing a result that is known to be infinity, no calculations needed. Some students may resent being asked to spend time doing the impossible, and the experience may make them less willing to embark on an exploration in the next exercise.
I don’t understand why everyone tries to make a calamity out of this… 
The students experienced in programming and calculus, they know how to calculate a series and how to check if it is divergent or convergent. It was a simple example to show different ways of adding up numbers. The example was never meant to be a 100% accurate scientific calculation, but it shows very clearly how important it is to think about the order of summation (the difference with real32 is 18.8 - 15.4, whereas the sum over 1/i**2 only has a difference of 1.6449 - 1.6447).
The only thing I wanted to share is an edge-case, where compilers don’t behave like they should according to the standard. So please forget about the lecturer, the students and the example.
Next year I will be the lecturer. So if you’ve got a better idea, how to explain the difference between forward and backward summation, please let me know. I would simply ask the students to calculate the sum from 1 to 1e9 over 1/i. To show them they how to approximate a converging series (which they should already know) I’ll let them calculate the Leibniz-formular for pi.
1 Like
The example I have used in similar circumstances is the Logistic map, which, for very small changes in the parameter mu can display very different behaviour, and can be used to talk about sensitivity to initial conditions, and iterated maps more generally, for example the iterated maps we use to solve Odes numerically, just a thought
1 Like
I wrote a code for the sum of 1/i^2
program main
! demonstrate that summing positive floats in ascending order
! is more accurate
implicit none
integer :: i
integer, parameter :: n = huge(i) - 1, dp = kind(1.0d0)
real :: sum_ascend = 0.0, sum_descend = 0.0
real(dp), parameter :: pi = 3.141592653589793238462643_dp
do i=1,n
sum_descend = sum_descend + 1/(real(i)**2)
sum_ascend = sum_ascend + 1/(real(n-i+1)**2)
end do
print "(3a20)","limit","ascending_sum","descending_sum"
print "(3f20.8)",pi**2/6,sum_ascend,sum_descend
! sum of the infinite series 1/n^2 is pi^2/6
! https://www.quora.com/What-is-the-sum-of-the-series-sum-dfrac-1-n-2
end program main
! gfortran output:
! limit ascending_sum descending_sum
! 1.64493407 1.64493406 1.64472532
! Intel Fortran output:
! limit ascending_sum descending_sum
! 1.64493407 1.64493406 1.64482379
Knowing what the truncated sums are for various N, I believe there are ways to extrapolate to the N=Inf case.
Because the Fortran standard is “funny” that way. If one wants to use HUGE integer of one kind, say X, the wording in the standard effectively suggests using an integer kind that has a higher range than the one with kind X. So you can try something like this:
use, intrinsic :: iso_fortran_env, only : I8 => int64
integer(I8), parameter :: HUGEI = huge(1)
integer(I8), parameter :: INTVL = huge(1)/50
integer(I8) :: i
do i = 1, HUGEI
if ( (i==1).or.(mod(i,INTVL)==0) ) then
print *, "i = ", i
end if
end do
print "(a,i0)", "'do i=1,huge(i)' exits with i=", i
end
1 Like
So in other words, one should consider huge(i) to not belong to the range of numbers one would use for i. Or to be extra cautious when using it for i.
Do you know where in the standard I can find this wording?
The comparison of the two compilers is interesting by itself. Not only is the Intel’s result (descending) pretty closer to the real value but also, at least on my Ubuntu 20.04, i5-10505 CPU machine, more than 4 times faster (0.97 vs. 4.3 seconds) with no options to compiler, 2 times with “-O3” (0.96 vs. 2.04 sec)
Using a larger integer type does work and is standard in my example above as well and it acts “as expected” for the INT32 type by using an INT64 type with the caveats mentioned; but in the case of the series it would be interesting to show how a simple sum would quit changing when the delta becomes smaller than epsilon and how to make a more accurate sum, and how a simple-looking array syntax solution might use a massive amount of memory just as examples of the differences between mathematics and computing, as well as showing the index overflow issue in the example if that type of issue is part of the lecture (else if might scare them away from computation if more introductory!)
The descending sum is less accurate in this case, but gets as close as it will ever get in just 4096 iterations, which is lot less than huge(0_int32) which I think is interesting; and
calculating what N you should start with instead of 1 for the ascending sum is cool, and shows how many cycles you can crunch without much return ; and that this particular problem has a determinant answer you can compute with a single line and to always use parens around powers (if new to programming, what would most people think pi**2/6 should be?).
So:
- why this would be a bad idea
sum_ascend=sum([(1.0d0/(real(n-i+1)**2),i=1,n)])
- why for the descending sum does this give the same answer?
do i=1,4096
sum_descend = sum_descend + 1/(real(i)**2)
end do
print "(3f20.8)",pi**2/6.0d0,sum_ascend,sum_descend
as “do i=1,n”?
- and so on. So look at the timing differences (with default compiler switches):
module M_tictoc
use,intrinsic::iso_fortran_env,only:real32,real64
implicit none
private
public :: tic
public :: toc
real(kind=real64) :: t1,t2
contains
subroutine tic()
call cpu_time(t1)
end subroutine tic
subroutine toc()
call cpu_time(t2)
print*,"Time Taken -->", real(t2-t1)
end subroutine toc
end module M_tictoc
program main
! demonstrate that summing positive floats in ascending order
! is more accurate
use,intrinsic::iso_fortran_env,only:int8,int16,int32,int64
use M_tictoc
implicit none
integer(kind=int64) :: i ! WHY IS THIS INT64?
integer, parameter :: dp = kind(1.0d0)
integer, parameter :: n=huge(0)
real :: sum_ascend = 0.0, sum_descend = 0.0
real(dp), parameter :: pi = 3.141592653589793238462643_dp
call tic()
print*, 'IDEAL',pi**2/6.0d0
call toc()
print*, 'TAKING INTO ACCOUNT PRECISION '
sum_ascend=0.0
sum_descend=0.0
call tic()
do i=1,4096
sum_descend = sum_descend + 1/(real(i)**2)
end do
call toc()
call tic()
do i=n-31663061,n
sum_ascend = sum_ascend + 1/(real(n-i+1)**2)
end do
call toc()
print "(3a20)","limit","ascending_sum","descending_sum"
print "(3f20.8)",pi**2/6.0d0,sum_ascend,sum_descend
print*, 'MATHEMATICAL APPROACH WITH A FINITE LIMIT'
sum_ascend=0.0
sum_descend=0.0
call tic()
do i=1,n
sum_descend = sum_descend + 1/(real(i)**2)
end do
call toc()
call tic()
do i=1,n
sum_ascend = sum_ascend + 1/(real(n-i+1)**2)
end do
call toc()
print "(3a20)","limit","ascending_sum","descending_sum"
print "(3f20.8)",pi**2/6.0d0,sum_ascend,sum_descend
end program main
IDEAL 1.6449340668482264
Time Taken --> 3.99999990E-05
TAKING INTO ACCOUNT PRECISION
Time Taken --> 7.00000010E-06
Time Taken --> 5.66099994E-02
limit ascending_sum descending_sum
1.64493407 1.64493406 1.64472532
MATHEMATICAL APPROACH WITH A FINITE LIMIT
Time Taken --> 3.88258600
Time Taken --> 3.84209609
limit ascending_sum descending_sum
1.64493407 1.64493406 1.64472532
- why in some old code you see something like
function my_dp_eps()
! calculate the epsilon value of a machine the hard way
real(kind=dp) :: t
real(kind=dp) :: my_dp_eps
! starting with a value of 1, keep dividing the value
! by 2 until no change is detected. Note that with
! infinite precision this would be an infinite loop,
! but floating point values in Fortran have a defined
! and limited precision.
my_dp_eps = 1.0d0
SET_ST: do
my_dp_eps = my_dp_eps/2.0d0
t = 1.0d0 + my_dp_eps
if (t <= 1.0d0) exit
enddo SET_ST
my_dp_eps = 2.0d0*my_dp_eps
end function my_dp_eps
and why there are functions like EPSILON() and RANGE().
- how the standard does not define the SUM() command to a level that it is guaranteed to use a remainder method or binning or … so why you want to do
some low-level things yourself sometimes.
That should put the scare in them enough to read about working with floating-point values.
1 Like
@Carltoffel , you’ll find it is spread across the standard: an underlying theme you will notice is that the onus is on the program writer to ensure the evaluation of expressions is possible by the types and kinds of the objects involved in such expressions. Thus it’s a trivial consideration then that if a writer wants to use HUGE integer of one kind in a program where said limiting value can then enter into expressions that are not representable by the kinds chosen for the integer variables, then practically it will be up to the writer to make adjustments or contend with whatever processor-dependent and unexpected behavior might arise.
As mentioned upthread, the expression toward the iteration count can fail to be representable by the integer kind chosen for the loop iteration-variable. Hence my suggestion above where an integer with a higher range is chosen as the loop iteration-variable to ensure the iteration count, which the standard says has the same kind as the iteration-variable becomes representable as a result of the expression evaluation.
So summing to infinity using real numbers in a DO loop is not such a trivial example !
It is interesting that, from my approach to testing, the accumulated sum results are not affected by -ffast-math rounding.
The ascending sum converges after 2^25 itterations for the real32, where 1/(real(i)**2) = 1.e-15.
The decending sum converges (incorrectly) for 2^12 itterations, where 1/(real(i)**2) = 6.e-8 (close to epsilon)
This use of assending sum would be limited to this class of problem where many values are to be accumulated and each xi can be calculated to high accuracy proportional to each xi value.
For most sum of xi’s where the accuracy of each xi is similar, ascending sum is not as effective.
It is important to note that the accuracy of the accumulated sum is only as accurate as the least accurate value being accumulated, being about 1.e-8 for this real32 sum.