Additive underflow occurs when two numbers of different magnitudes are added or subtracted. The processor adjusts the exponent of the number with the smaller magnitude so that the two exponents are the same. To do this, the mantissa of the number with the smaller magnitude is shifted to the right. The least significant bits of this number are therefore rounded in the result. If the magnitudes are very different, the number with the smaller magnitude may disappear completely from the result.
The program shown below describes a space ship. The ship starts from a starting position of 0 metres and moves at 1/3 metre per second. Every second, we add 1/3 of a metre to its position. It is moving in empty space, and it should go on moving forever. However, in our simulation it doesn't. The distance travelled per second becomes so small in comparison with the distance travelled that the result of the addition stops changing, and the spaceship stops moving.
! ***************************************************************************** ! spaceship.f90 8-Oct-15 John Collins ! ! A simple demonstration of additive underflow ! ! ***************************************************************************** PROGRAM spaceship IMPLICIT NONE INTEGER,PARAMETER :: kr4 = KIND(1.0E0) REAL(KIND=kr4) :: dist,step,prev_dist ! metres INTEGER(KIND=8):: time ! seconds REAL(KIND=kr4) :: buffer(20) ! Last 20 positions INTEGER buf_index,i,j time = 0 dist = 0.0 buf_index = 0 ! We travel at one third of a metre per second step = 1.0/3.0 prev_dist = -step DO WHILE (dist > prev_dist) ! Move the spaceship time = time+1 prev_dist = dist dist = dist + step ! Capture the last 20 positions buf_index = buf_index+1 IF (buf_index > 20) buf_index=1 buffer(buf_index) = dist ENDDO WRITE(*,'(/,"Distance travelled per second: ",E16.10)')step WRITE(*,'(/," Time (secs) Distance Change distance Time*step",/)') IF (buf_index == 20) THEN i = 1 ELSE i = buf_index+1 ENDIF DO j = 1,20 IF (j==1) THEN WRITE(*,'(I12,F16.3,16X,F16.3)')time-19,buffer(i),(time-19)*step ELSE WRITE(*,'(I12,3F16.3)')time-20+j,buffer(i),buffer(i)-prev_dist,(time-20+j)*step ENDIF prev_dist = buffer(i) i = i+1 IF (i > 20) i = 1 ENDDO END PROGRAM spaceship
We capture the position of the ship in a circular buffer, and the program reports the last 20 positions:
Distance travelled per second: 0.3333333433E+00 Time (secs) Distance Change distance Time*step 22805181 8388599.000 7601727.000 22805182 8388599.500 0.500 7601727.500 22805183 8388600.000 0.500 7601728.000 22805184 8388600.500 0.500 7601728.000 22805185 8388601.000 0.500 7601728.000 22805186 8388601.500 0.500 7601729.000 22805187 8388602.000 0.500 7601729.500 22805188 8388602.500 0.500 7601729.500 22805189 8388603.000 0.500 7601729.500 22805190 8388603.500 0.500 7601730.000 22805191 8388604.000 0.500 7601731.000 22805192 8388604.500 0.500 7601731.000 22805193 8388605.000 0.500 7601731.000 22805194 8388605.500 0.500 7601731.500 22805195 8388606.000 0.500 7601732.000 22805196 8388606.500 0.500 7601732.000 22805197 8388607.000 0.500 7601732.000 22805198 8388607.500 0.500 7601733.000 22805199 8388608.000 0.500 7601733.500 22805200 8388608.000 0.000 7601733.500
The change in distance at the last time step is zero. Our space ship has stopped moving. Note also that the change in position in the preceding 19 time steps is half a metre, even though we are adding 1/3 metre every second. Half a metre represents the least significant bit of the mantissa of the distance travelled. It is the smallest change that the program can make.
The last column shows the distance the ship should have travelled. The result is accurate to better than 6 decimal places.
The program uses 32-bit real numbers with 23 mantissa bits. If we had used 64-bit real numbers, with 52 mantissa bits, our space ship would still stop, but it would travel 2**29 times further before it did so.