PROGRAM MachEps ! ------------------------------------------------------------------------------ ! ! Copyright 2006 Matthew Norton ! ! ------------------------------------------------------------------------------ ! This program is free software; you can redistribute it and/or modify ! it under the terms of the GNU General Public License as published by ! the Free Software Foundation; either version 2 of the License, or ! (at your option) any later version. ! ! This program is distributed in the hope that it will be useful, ! but WITHOUT ANY WARRANTY; without even the implied warranty of ! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the ! GNU General Public License for more details. ! ! You should have received a copy of the GNU General Public License ! along with this program; if not, write to the Free Software ! Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA ! ! NOTE: This is in FORTRAN 90 and you will need a FORTRAN 90 complier to run ! this code ! ------------------------------------------------------------------------------ ! ! This program calculates the machine epislon, the smallest number that, when ! added to 1, gives a number greater than 1. This is useful in finding out the ! precision of the computer you run your FORTRAN program on. ! ! ------------------------------------------------------------------------------ ! ! PROGRAM HISTORY ! ! AUTHOR: Matthew S. Norton ! CREATION DATE: January 18, 2007 ! ! ------------------------------------------------------------------------------ IMPLICIT NONE REAL*8 :: Eps, EpsA REAL*8, PARAMETER :: A = 1 CHARACTER(1) :: ForS, FileName*20 INTEGER :: OutUnit, Man 10 PRINT *, "Do you want to output to a file (f) or the screen (s)?" READ *, ForS SELECT CASE(ForS) CASE("F","f") OutUnit = 20 PRINT *, "What do you want the output file to be (include extension)?" READ *, FileName OPEN (OutUnit, File = TRIM(FileName)) CASE("S","s") OutUnit = 6 Case Default GOTO 10 END SELECT Eps = 1 EpsA = 0 Epislon: DO EpsA = Eps + A IF (EpsA > 1.) THEN WRITE(OutUnit, '(1X, I1, 1x, A1, 1X, E15.10E2, 1X, A1, 1X, E15.10E2)') 1,"+",Eps, "=", 1+Eps Eps = Eps / 2. ELSE IF (EpsA == 1.)THEN Eps = Eps * 2. EXIT Epislon END IF END DO Epislon WRITE(OutUnit, '(1X, A39, 1X, E20.14E2)') "The machine epislon of your machine is:", Eps !Man = 1. - NINT(ALOG10(2.*Eps)) !WRITE(OutUnit, '(1X, A32, 1X, I3)') "The mantissa of your machine is:", Man END PROGRAM