# Difference between revisions of "Floating point math issues"

From Geos-chem

(→The pitfalls of floating point mathematics) |
(→The pitfalls of floating-point mathematics) |
||

Line 4: | Line 4: | ||

Here are some excellent references about floating-point math: | Here are some excellent references about floating-point math: | ||

− | * [http://www.lahey.com/float.htm | + | * [http://www.lahey.com/float.htm ''The Perils of Floating Point'' by Bruce M. Bush] |

* [http://blogs.msdn.com/ericlippert/archive/2005/01/20/fun-with-floating-point-arithmetic-part-five.aspx Fun with Floating Point Arithmetic] | * [http://blogs.msdn.com/ericlippert/archive/2005/01/20/fun-with-floating-point-arithmetic-part-five.aspx Fun with Floating Point Arithmetic] | ||

## Revision as of 14:39, 9 April 2008

## The pitfalls of floating-point mathematics

One important thing to remember when writing computer programs is that floating-point mathematics is not exact but is an approximation. In most programming languages, floating-point real numbers are composed of groups of 4 or 8 bytes. This means that floating-point numbers are not infinitely precise, but have a maximum precision. As a consequence, floating-point math operations (especially multiplication and division) can often lead to different results than anticipated.

Here are some excellent references about floating-point math:

## Safe floating-point division

From http://groups.google.com/group/comp.lang.fortran/browse_thread/thread/8b367f44c419fa1d/

- Thank you to all that offered their suggestions for a "safe division" routine to prevent overflow. For those that are curious about the solution to the problem, I found useful to adopt a subroutine along the lines suggested by William Long:

if(exponent(a) - exponent (b) >= maxexponent(a) .OR. b==0)Then q=altv else q=a/b endif

- The = in the >= is to take into account the case when the fractional part of a is 1.111... and that of b is 1.

- It works very well.

- Thank you again,
- Julio.

We have implemented the above algorithm into GEOS-Chem. See routine SAFE_DIV in error_mod.f:

FUNCTION SAFE_DIV( N, D, ALTV ) RESULT( Q ) ! !****************************************************************************** ! Subroutine SAFE_DIV performs "safe division", that is to prevent overflow, ! underflow, NaN, or infinity errors. An alternate value is returned if the ! division cannot be performed. (bmy, 2/26/08) ! ! For more information, see the discussion on: http://groups.google.com/group/comp.lang.fortran/browse_thread/thread/8b367f44c419fa1d/ ! ! Arguments as Input: ! ============================================================================ ! (1 ) N : Numerator for the division ! (2 ) D : Divisor for the division ! (3 ) ALTV : Alternate value to be returned if the division can't be done ! ! NOTES: !****************************************************************************** ! ! Arguments REAL*8, INTENT(IN) :: N, D, ALTV ! Function value REAL*8 :: Q !================================================================== ! SAFE_DIV begins here! !================================================================== IF ( EXPONENT(N) - EXPONENT(D) >= MAXEXPONENT(N) .or. D==0 ) THEN Q = ALTV ELSE Q = N / D ENDIF ! Return to calling program END FUNCTION SAFE_DIV