The Fortran function which rounds a 4-byte real number is shown below. Note that the function is elemental, which indicates that it can operate on whole arrays. The function rounds the 4-byte real number to the number of mantissa bits specified by the Fortran parameter kr4_mantissa_bits. This parameter is specified in the include file precision_experiments_dec.i90 and both this file and the code of the rounding function are included in the Fortran module module_emulate_real_arithmetic. The parameter kr4_mantissa_bits is used to set up additional parameters used to mask off the unused bits of the operand and to round the result. Only the parameter kr4_mantissa_bit needs to be modified to adjust the precision of 4-byte real numbers.
The declarations of the Fortran parameters (constants) are:
! ********************************************************************************************************************* ! precision_experiments_dec.i90 7-Sep-15 John Collins ! ********************************************************************************************************************* INTEGER,PARAMETER :: kr4_mantissa_bits = 23 INTEGER,PARAMETER :: kr8_mantissa_bits = 52 INTEGER(KIND=ki4),PARAMETER :: kr4_test_bit = ISHFT(1_ki4,(22-kr4_mantissa_bits)) INTEGER(KIND=ki4),PARAMETER :: kr4_mask = IEOR(-1_ki4,ISHFT(1_ki4,(23-kr4_mantissa_bits))-1) INTEGER(KIND=ki4),PARAMETER :: kr4_exponent_lsb = Z'00800000' INTEGER(KIND=ki4),PARAMETER :: kr4_exponent_mask = Z'FF800000' INTEGER(KIND=ki4),PARAMETER :: kr4_max_mantissa = IAND(Z'004FFFFF',(kr4_mask+kr4_test_bit)) INTEGER(KIND=ki8),PARAMETER :: kr8_test_bit = ISHFT(1_ki8,(51-kr8_mantissa_bits)) INTEGER(KIND=ki8),PARAMETER :: kr8_mask = IEOR(-1_ki4,ISHFT(1_ki8,(52-kr8_mantissa_bits))-1) INTEGER(KIND=ki8),PARAMETER :: kr8_exponent_lsb = Z'0010000000000000' INTEGER(KIND=ki8),PARAMETER :: kr8_exponent_mask = Z'FFF0000000000000' INTEGER(KIND=ki8),PARAMETER :: kr8_max_mantissa = IAND(Z'000FFFFFFFFFFFFF',(kr8_mask+kr8_test_bit)) PUBLIC show_parameters ! End of precision_experiments_dec.i90 ********************************************************************************
The function which rounds the 4-byte emulated real values is:
ELEMENTAL FUNCTION experiment_4(r) IMPLICIT NONE TYPE(em_real_k4),INTENT(IN) :: r REAL(KIND=kr4) :: experiment_4 INTEGER(KIND=ki4)ir ir = TRANSFER(r,ir) IF (IAND(ir,kr4_test_bit) /= 0) THEN IF (IAND(ir,kr4_test_field) /= 0) THEN IF (IAND(ir,kr4_max_mantissa) == kr4_max_mantissa) THEN ir = IAND(ir,kr4_exponent_mask) + kr4_exponent_lsb ELSE ir = ir+kr4_test_bit ENDIF ENDIF ENDIF ir = IAND(ir,kr4_mask) experiment_4 = TRANSFER(ir,experiment_4) END FUNCTION experiment_4
Note that when no bits in the mantissa beyond the test bit are set the mantissa is not rounded up. This is the IEEE behaviour. Thanks are due to Andrew Dawson at Atmospheric, Oceanic and Planetary Physics, Oxford, for pointing out this detail and correcting an error in the original code.