43 lines
1.4 KiB
Fortran
43 lines
1.4 KiB
Fortran
program main
|
|
implicit none
|
|
|
|
integer :: i
|
|
double precision, allocatable :: vals(:)
|
|
|
|
vals = [ 100000000000000.01d0, 100000000000000.011d0, &
|
|
& 100.01d0, 100.011d0, &
|
|
& 10000000000000.001d0/10000d0, 1000000000.0000001000d0, &
|
|
& 0.001d0, 0.0010000001d0, &
|
|
& 0.000000000000000000000101d0, 0d0, &
|
|
& sqrt(2d0)*sqrt(2d0), 2d0, &
|
|
& -sqrt(2d0)*sqrt(2d0), -2d0, &
|
|
& 3.14159265358979323846d0, 3.14159265358979324d0 ]
|
|
|
|
do i = 1, size(vals)/2
|
|
print '(ES30.18, A, ES30.18, A, L)', vals(2*i-1), ' == ', vals(2*i), ' ? ', eq_approx(vals(2*i-1), vals(2*i))
|
|
end do
|
|
|
|
contains
|
|
|
|
logical function eq_approx(a, b, reltol, abstol)
|
|
!! is a approximately equal b?
|
|
|
|
double precision, intent(in) :: a, b
|
|
!! values to compare
|
|
double precision, intent(in), optional :: reltol, abstol
|
|
!! relative and absolute error thresholds.
|
|
!! defaults: epsilon, smallest non-denormal number
|
|
|
|
double precision :: rt, at
|
|
|
|
rt = epsilon(1d0)
|
|
at = tiny(1d0)
|
|
if (present(reltol)) rt = reltol
|
|
if (present(abstol)) at = abstol
|
|
|
|
eq_approx = abs(a - b) .le. max(rt * max(abs(a), abs(b)), at)
|
|
return
|
|
end function
|
|
|
|
end program
|