Блог

Я пытался найти совместимый со стандартами способ проверки бесконечности и значений NaN в Fortran 90/95, но это оказалось сложнее, чем я думал.

• Я попытался вручную создать переменные Inf и NaN, используя двоичное представление, описанное в IEEE 754, но я не нашел такой функциональности.
• Я знаю о встроенном ieee_arithmetic модуле в Fortran 2003 с ieee_is_nan() и ieee_is_finite() внутренними функциями. Однако он поддерживается не всеми компиляторами (особенно gfortran начиная с версии 4.9).

Определение infinity и NaN в начале, как pinf = 1. / 0 и nan = 0. / 0 , кажется мне хакерским, и ИМХО может вызвать некоторые проблемы со сборкой — например, если некоторые компиляторы проверяют это во время компиляции, нужно было бы предоставить специальный флаг.

Есть ли способ, который я могу реализовать в стандартном Fortran 90/95?

 function isinf(x)
! Returns .true. if x is infinity, .false. otherwise
...
end function isinf
  

и isnan() ?

Простой способ без использования ieee_arithmatic — сделать следующее.

Бесконечность: определите свою переменную infinity = HUGE(dbl_prec_var) (или, если она у вас есть, переменную с четырехкратной точностью). Затем вы можете просто проверить, равна ли ваша переменная infinity by if(my_var > infinity) .

NAN: Это еще проще. По определению, NAN не равен ничему, даже самому себе. Просто сравните переменную с самой собой: if(my_var /= my_var) .

У меня недостаточно репутации для комментариев, поэтому я «отвечу» на предложение Рика Томпсона по тестированию infinity.

 if (A-1 .eq. A) 
  

Это также будет верно, если A является очень большим числом с плавающей запятой и 1 ниже точности A.

Простой тест:

 subroutine test_inf_1(A)
    real, intent(in) :: A
    print*, "Test (A-1 == A)"
    if (A-1 .eq. A) then
        print*, "    INFINITY!!!"
    else
        print*, "    NOT infinite"
    endif
end subroutine

subroutine test_inf_2(A)
    real, intent(in) :: A
    print*, "Test (A > HUGE(A))"
    if (A > HUGE(A)) then
        print*, "    INFINITY!!!"
    else
        print*, "    NOT infinite"
    endif
end subroutine


program test
    real :: A,B

    A=10
    print*, "A = ",A
    call test_inf_1(A)
    call test_inf_2(A)
    print*, ""

    A=1e20
    print*, "A = ",A
    call test_inf_1(A)
    call test_inf_2(A)
    print*, ""

    B=0.0 ! B is necessary to trick gfortran into compiling this
    A=1/B
    print*, "A = ",A
    call test_inf_1(A)
    call test_inf_2(A)
    print*, ""

end program test
  

выходы:

 A =    10.0000000    
Test (A-1 == A)
    NOT infinite
Test (A > HUGE(A))
    NOT infinite

A =    1.00000002E 20
Test (A-1 == A)
    INFINITY!!!
Test (A > HUGE(A))
    NOT infinite

A =          Infinity
Test (A-1 == A)
    INFINITY!!!
Test (A > HUGE(A))
    INFINITY!!!
  

Нет.

Основные части IEEE_ARITHMETIC для генерации / проверки NaN достаточно легко написать для gfortran для конкретной архитектуры.

Я использовал:

   PROGRAM MYTEST
  USE, INTRINSIC :: IEEE_ARITHMETIC, ONLY: IEEE_IS_FINITE      
  DOUBLE PRECISION :: number, test
  number = 'the expression to test'
  test = number/number
  IF (IEEE_IS_FINITE(test)) THEN
     WRITE(*,*) 'We are OK'
  ELSE
     WRITE(*,*) 'Got a problem'
  END IF         
     WRITE(*,*) number, test
  END PROGRAM MYTEST
  

При этом будет выведено ‘Возникла проблема’ для number = 0.0D0, 1.0D0 /0.0D0, 0.0D0 /0.0D0, SQRT(-2.0D0), а также для переполнений и недопотоков, таких как number = EXP(1.0D800) или number = EXP(-1.0D800). Обратите внимание, что обычно такие вещи, как number = EXP (1.0D-800), просто устанавливают number = 1.0 и выдают предупреждение во время компиляции, но программа напечатает «У нас все в порядке», что я нахожу приемлемым.

OL.

для тестирования NaN ничего не сработало, например. если тестировать реальный s2p, чтобы увидеть, является ли это NaN, тогда

 if(isnan(s2p)) 
  

не работает в gfortran и не

 if(s2p.ne.s2p). 
  

Единственное, что сработало, это

 if(.not.s2p<1.and..not.s2p>1)  
  

хотя, чтобы убедиться, что вы можете добавить

 if(.not.s2p<1.and..not.s2p>1.and..not.s2p==1)    
  

Нет.

В Fortran 90/95 также нет способа проверки бесконечности или NAN, соответствующего стандартам, и не может быть способа, соответствующего стандартам. В Fortran 90/95 нет совместимого со стандартами способа определения любого из этих квазичислов.

Для Fortran 1 / infinity = 0 таким образом, разделите вашу переменную на ноль, т. Е.

 program test
implicit none
real :: res
integer :: i

do i=1,1000000
    res=-log((1. (10**(-real(i))))-1.)
    print*,i,res
    if ((1./res)==0.) then
        exit
    end if
end do

end program
  

вот ваша проверка на бесконечность. Никаких усложнений не требуется.

Для Inf, похоже, работает, что если (A-1 .eq. A) верно, тогда A является Inf