Peryapse/macheps.f at master · drbitboy/Peryapse · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C Purpose:  approximate [machine+Fortran compiler] epsilon
C
C Build:  gfortran -static macheps.f -o macheps64lebeau
C Build:  g77 -static macheps.f -o macheps32hogan (hogan)
C
C Run:  ./macheps32hogan
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      program macheps
      implicitnone
      doubleprecision dble_epsi
      doubleprecision dble_eps

      doubleprecision dble_dble
      doubleprecision dble_one / 1d0 /

      integer*8 i8_dble
      integer*8 i8_one

      equivalence (dble_dble,i8_dble)
      equivalence (dble_one,i8_one)

      real*4 real_eps
      real*4 real_epsi

      real*4 real_real
      real*4 real_one / 1.0 /

      integer*4 i4_real
      integer*4 i4_one

      equivalence (real_real,i4_real)
      equivalence (real_one,i4_one)

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

      dble_eps = 1d0
      dble_dble = 1.5d0
      do while (dble_dble .ne. dble_one)
        dble_eps = dble_eps * 0.5d0
        dble_dble = 1d0+(0.5d0*dble_eps)
      enddo

      dble_epsi = 1d0
      dble_dble = 1.5d0
      do while (i8_dble .ne. i8_one)
        dble_epsi = dble_epsi * 0.5d0
        dble_dble = 1d0+(0.5d0*dble_epsi)
      enddo

      real_eps = 1.0
      real_real = 1.5
      do while (real_real .ne. real_one)
        real_eps = real_eps * 0.5
        real_real = 1.0+(0.5*real_eps)
      enddo

      real_epsi = 1.0
      real_real = 1.5
      do while (i4_real .ne. i4_one)
        real_epsi = real_epsi * 0.5
        real_real = 1.0+(0.5*real_epsi)
      enddo

      print'(1x,e24.16,a)',dble_eps ,'~Double epsilon (DOUBLE compare)'
     &                    ,dble_epsi,'~Double epsilon (INT*8 compare)'
     &                    ,real_eps ,'~Single epsilon (REAL compare)'
     &                    ,real_epsi,'~Single epsilon (INT*4 compare)'

      dble_eps = 1d0
      do while ((1.0d0+(0.5d0*dble_eps)) .ne. dble_one)
        dble_eps = dble_eps * 0.5d0
      enddo

      real_eps = 1.0
      do while ((1.0+(0.5*real_eps)) .ne. real_one)
        real_eps = real_eps * 0.5
      enddo

      print'(1x,e24.16,a)',dble_eps ,'~Double epsilon (Expr. compare)'
     &                    ,real_eps ,'~Single epsilon (Expr. compare)'

      end