Странное поведение Fortran 90

program main
call findbracket(x0, a, b)
end program

function f(x)
double precision x,f
f = x
end function

subroutine findbracket(x0,a,b)
double precision x0, a, b
double precision fa, fb
double precision dx
dx = 0.001d0
x0 = 1.0d0
a = x0
b = x0
print*, a, b
print*, f(a)
print*, f(b)

do
fa = f(a)
fb = f(b)

print*, "what is fa", fa
print*, "what is fb", fb
a = a - dx

if (fa*fb < 0) then
        exit
end if
print*, b, dx
b = b + dx
if (fa*fb < 0) then
        exit
end if

dx = dx*2

end do

end subroutine 

Я пишу программу, которая решает корень заданной функции f (x). На этом этапе я просто хочу проверить, получаю ли я правильное значение для каждого шага, и заметил некоторые ошибки. Чтобы упростить задачу, у меня есть функция f (x) = x, и я ожидал, что программа напечатает f (a) = 1.0 и f (b) = 1.0, но результат будет f (a) = 1.8750000000000000 и f (b ) = 0,0000000000000000. Плюс, хотя я установил как a, так и b равными x0, кажется, что это a = 1.0000002381857485 и b = 1.0000000000000000. Может ли кто-нибудь объяснить, почему это происходит? Возможно, я упускаю какую-то глупость, но я не могу этого найти. Я ценю вашу помощь.