수치 정밀도와 오차 - 결과를 믿어도 되나

이 장에서 배우는 것

앞 장에서 배열의 저장 공간과 값을 다루는 방법을 살펴보았다. 이제 배열 안에 들어 있는 실수가 계산 결과를 얼마나 충실하게 나타내는지 확인한다. 2차원 금속판의 온도를 배열에 저장했다고 해서 그 값이 언제나 믿을 만한 것은 아니다. 작은 온도 변화가 덧셈에서 사라질 수도 있고, 계산 중 생긴 유효하지 않은 값이 판 전체로 퍼질 수도 있다.

이 장에서는 열 방정식의 계산식을 아직 도입하지 않는다. 대신 작은 온도 변화량을 여러 번 더하는 실험으로 정밀도를 비교하고, 온도 배열의 상태를 검사하며, 기준값과 계산값을 오차 허용 범위 안에서 비교한다. 이후 수치 계산을 구현할 때 사용할 검증 도구를 먼저 갖추는 셈이다.

• 실수 종류(kind)를 선택하고 iso_fortran_env의 이름으로 계산 정밀도를 일관되게 지정한다.
• 작은 변화량을 반복해서 더할 때 반올림 오차가 어떻게 쌓이는지 설명한다.
• ieee_arithmetic으로 NaN과 무한대를 구분하고 온도 배열에서 비유한 값을 찾는다.
• 절대 허용 오차와 상대 허용 오차를 함께 사용하여 두 실수를 비교한다.

문제 상황

금속판의 초기 온도를 20도로 두고, 중앙 격자점에 매우 작은 온도 상승량을 여덟 번 더한다고 하자. 저장한 변화량은 양수이고 반복문도 정상적으로 실행된다. 그런데 결과를 확인하면 어떤 실수 종류에서는 온도가 여전히 20도다. 출력 자릿수가 부족한 것인지, 덧셈 자체에서 변화가 사라진 것인지 구분해야 한다.

다른 날에는 온도 배열의 일부에 NaN이 들어온다. 잘못된 입력값이나 유효하지 않은 연산이 원인일 수 있다. 이런 값을 이웃 격자점 계산에 그대로 사용하면 다음 계산에서도 NaN이 만들어질 수 있다. 최종 결과만 살펴보면 문제가 처음 발생한 위치를 찾기 어렵다.

검증 코드에도 문제가 생긴다. 기준 온도와 계산 온도가 화면에서는 똑같이 보이지만, == 비교는 거짓을 반환한다. 반대로 소수 둘째 자리까지만 출력하면 의미 있는 오차가 가려질 수 있다. 화면에 보이는 자릿수, 내부 표현의 정밀도, 결과를 받아들일 허용 오차는 각각 다른 기준이다.

필요한 것은 실수 종류를 한곳에서 정하는 규칙, 유효하지 않은 값을 일찍 찾는 검사, 수치의 크기에 맞는 비교 방법이다. 이 세 가지를 마련하면 결과가 다르게 나왔을 때 원인을 좁혀 갈 수 있다.

실수 종류와 표현 가능한 간격

Fortran의 kind 값은 자료형의 표현 방식을 식별하는 정수다. 이 정수를 바이트 수라고 해석하면 안 된다. 예를 들어 real(kind=8)이라는 선언은 종류 번호가 8인 실수라는 뜻이며, 언어 자체가 이를 8바이트 실수로 정의하지는 않는다. 컴파일러가 제공하는 종류 번호와 저장 크기를 구분해야 한다.

iso_fortran_env의 real32와 real64는 각각 저장 크기가 32비트와 64비트인 실수 종류를 가리킨다. 해당 종류를 지원하지 않는 처리계에서는 이름의 값이 음수다. 이 장의 GNU Fortran 실행 환경에서는 두 종류를 사용할 수 있다. 다만 저장 크기 이름만으로 모든 처리계의 기수와 유효 자릿수가 같다고 가정하지는 않는다.

금속판 계산에는 integer, parameter :: rk = real64처럼 공통 이름을 정한다. 온도 배열, 변화량, 허용 오차에 같은 이름을 사용하면 정밀도를 바꾸거나 선언을 검토하기 쉽다. 여러 모듈로 프로그램을 나눌 때에도 이 이름을 공통 모듈에서 제공하는 방식으로 확장할 수 있다.

use iso_fortran_env, only : real64
integer, parameter :: rk = real64

real(rk) :: temperature
temperature = 20.0_rk

리터럴에도 종류 접미사를 붙인다. 0.1은 기본 실수 종류로 먼저 해석된다. 이를 나중에 real64 변수에 대입해도 처음 표현하면서 생긴 오차가 없어지지는 않는다. 0.1_real64처럼 쓰면 처음부터 선택한 종류로 상수를 표현한다. 반면 20처럼 이진 표현이 정확한 값을 골랐을 때에는 기본 종류를 거쳐도 같은 값이 될 수 있다. 이런 우연에 의존하지 않고 작성 규칙을 통일하는 편이 좋다.

실수 종류와 계산 특성을 확인하는 이름

이름	의미	해석할 때 주의할 점
real32, real64	저장 크기로 지정한 실수 종류	종류 번호 자체가 비트 수는 아니다.
precision(x)	십진 정밀도	출력할 소수 자리 수와 다르다.
epsilon(x)	1보다 큰 다음 표현값과 1의 차이	모든 크기의 값에서 간격이 같지는 않다.
spacing(x)	x 부근의 모델 수 간격	작은 변화량이 저장될 수 있는지 살펴보는 단서다.

이진 부동소수점(floating point)은 표현 가능한 값 사이에 간격이 있다. 이 장의 환경에서 real32의 20 부근 간격은 2의 −19승이다. 이에 비해 이번 실험의 변화량은 2의 −20승이다. 변화량이 간격의 절반이므로 덧셈 결과는 두 표현값의 가운데에 놓인다.

가장 가까운 값으로 반올림하고 가운데에서는 끝 비트가 짝수인 쪽을 선택하는 기본 반올림 방식에서는 이 덧셈이 다시 20으로 저장된다. 여덟 번 반복해도 매번 저장된 20에서 출발하므로 온도는 올라가지 않는다. real64에서는 같은 변화량을 표현하고 더할 수 있어 여덟 번의 증가가 남는다.

이 실험은 기본 반올림 모드를 사용하는 GNU Fortran 환경을 전제로 한다. 반올림 모드를 바꾸면 가운데 값의 처리 결과가 달라질 수 있다. 또한 real64로 바꾸었다고 반올림 오차가 사라지는 것은 아니다. 표현 간격이 더 작아져 이번 변화량을 보존할 수 있게 되었을 뿐이다.

반올림 오차와 결과 비교

반올림 오차(round-off error)는 정확한 연산 결과를 저장 가능한 값으로 바꾸면서 생긴 차이다. 반복 계산에서는 앞서 저장된 값이 다음 연산의 입력이 된다. 따라서 한 번의 오차가 다음 단계로 전달되고, 여러 오차가 같은 방향으로 쌓이거나 일부 상쇄될 수 있다.

이번 실험은 작은 증가량이 매번 사라지는 경우다. 일반적인 오차 누적이 언제나 반복 횟수에 정확히 비례하는 것은 아니다. 값의 크기, 연산 순서, 중간 결과의 부호가 영향을 준다. 실수 덧셈에서 (a + b) + c와 a + (b + c)가 다른 결과를 낼 수 있는 이유도 중간 단계의 반올림 때문이다.

예를 들어 크기가 큰 두 값이 거의 상쇄되는 계산에서는 작은 차이가 중요해진다. 큰 값 각각에 이미 포함된 오차가 마지막 차이에 비해 커질 수 있다. 정밀도를 높이는 것은 도움이 되지만, 계산식과 값의 규모를 검토하는 일도 필요하다. 온도 변화량을 합산할 때에는 큰 기준 온도와 작은 증가량을 어떤 순서로 더하는지 함께 살핀다.

반올림 오차와 수치 방법의 근사 오차도 구분한다. 공간이나 시간을 일정 간격으로 나누어 연속적인 현상을 근사하면, 실수를 더 정밀하게 저장해도 그 근사에서 생긴 오차는 남는다. 정밀도 변경과 격자 간격 변경은 서로 다른 실험이다. 이 장에서는 덧셈과 저장에서 생기는 오차에 집중한다.

절대 오차와 상대 오차를 함께 본다

절대 오차(absolute error)는 두 값의 차이의 크기인 abs(a - b)다. 온도를 비교하면 그 단위도 온도다. 상대 오차(relative error)는 이 차이를 비교 대상의 규모로 나눈 값이다. 이 장에서는 두 값 중 더 큰 절댓값을 규모로 삼아 대칭적인 비교를 만든다.

abs(a - b) <= max(atol, rtol * max(abs(a), abs(b)))

atol은 절대 허용 오차이고, rtol은 상대 허용 오차다. 위 조건은 절대 기준과 상대 기준 중 더 넓은 범위를 적용한다. 다른 프로그램에서는 두 허용량을 더하는 규칙을 사용하기도 한다. 어느 규칙을 택했는지 명시하고 검증 전반에 일관되게 사용해야 한다.

0 부근에서는 상대 기준만으로 비교하기 어렵다. 한 값이 0이고 다른 값이 아주 작아도, 두 값 중 큰 절댓값으로 나누면 상대 차이는 1이 된다. 따라서 0 부근에서는 의미 있는 절대 허용 오차가 필요하다. 반대로 크기가 큰 값에는 같은 절대 허용 오차가 지나치게 엄격할 수 있어 상대 기준이 유용하다.

완성 코드의 두 온도 차이는 약 0.0000076294도다. 상대 허용 오차를 1.0e-6_rk으로 두면 20도 부근에서 약 0.000020도의 차이를 허용하므로 비교가 참이다. 1.0e-8_rk으로 줄이면 상대 허용량은 약 0.0000002도가 되어 거짓이다. 같은 결과도 요구 정확도에 따라 판정이 달라진다.

허용 오차는 epsilon에서 자동으로 얻는 정답이 아니다. 측정 자료의 정확도, 계산 방법의 오차, 결과를 사용하는 목적을 보고 정한다. epsilon은 자료형의 표현 특성이고, 허용 오차는 사용자가 받아들이는 차이의 기준이다. 온도 단위나 기준점을 바꾸면 상대 비교의 의미도 달라질 수 있으므로 물리적 해석을 확인해야 한다.

위 비교식은 이 장처럼 값의 범위가 제한된 온도 검사에 사용한다. 매우 큰 유한 값에서는 a - b나 rtol * max(abs(a), abs(b))가 넘칠 수 있다. 넓은 지수 범위를 처리하는 범용 비교 함수라면 크기를 조정한 식과 경계 사례 검토가 추가로 필요하다.

NaN과 무한대를 계산에서 분리한다

NaN은 수로 표현할 수 없는 결과를 나타내는 값이다. 무한대(infinity)는 양의 무한대와 음의 무한대로 구분된다. IEEE 연산을 지원하는 환경에서는 유효하지 않은 연산이나 범위 초과 등으로 이러한 값이 생길 수 있다. 다만 예외를 중단하도록 실행 환경을 설정했다면 값을 받아 계속 진행하기 전에 프로그램이 멈출 수도 있다.

ieee_arithmetic은 이 상태를 검사하는 표준 모듈이다. ieee_is_nan(x)는 NaN인지 검사하고, ieee_is_finite(x)는 유한한 값인지 검사한다. 무한대만 찾고 싶다면 유한하지 않으면서 NaN도 아닌 값을 찾는다. 두 검사 함수는 배열에도 원소별로 적용할 수 있다.

완성 코드에서는 ieee_value로 조용한 NaN과 양의 무한대를 직접 만든다. 검사 동작을 확인하기 위한 의도적인 시험값이다. 0.0_rk / 0.0_rk 같은 식으로 NaN을 만들면 컴파일 단계에서 진단되거나 예외 설정의 영향을 받을 수 있다. 표준 기능으로 시험값을 만들면 무엇을 검사하는 코드인지 더 분명해진다.

ieee_support_nan과 ieee_support_inf로 선택한 실수 종류가 필요한 값을 지원하는지 먼저 확인한다. 지원하지 않으면 이 실험은 설명을 붙여 중단한다. 지원 여부를 확인하는 것과 예외 발생 이력을 확인하는 것은 다른 작업이다. 여기서는 배열에 현재 저장된 값의 상태를 검사한다.

NaN이 포함된 값을 일반적인 대소 비교로 처리하려 해서는 안 된다. 예를 들어 abs(a - b) > limit일 때만 오류로 판정하는 코드는, 차이가 NaN이면 조건이 거짓이 되어 오류를 놓칠 수 있다. 비교 함수는 먼저 두 입력이 유한한지 검사하고, 그다음 오차를 계산해야 한다.

무한대도 정상 온도 비교에서는 제외한다. 같은 부호의 무한대 두 개를 빼면 유효한 온도 차이가 나오지 않는다. 이 장의 비교 함수는 입력값이나 허용 오차가 비유한 값이면 거짓을 반환한다. 음수 허용 오차도 거부한다. 따라서 거짓은 차이가 너무 크다는 뜻뿐 아니라 비교 조건이 유효하지 않다는 뜻일 수도 있다.

실제 시뮬레이터에서는 비유한 값을 찾았을 때 위치와 계산 단계를 기록한 뒤 중단하는 정책을 둘 수 있다. 이 장에서는 검사 결과를 보여 주기 위해 개수만 출력한다. 비유한 값을 0으로 바꾸어 계산을 계속하는 처리는 하지 않는다. 그렇게 바꾸면 열량이 갑자기 달라지고 원래 오류가 가려질 수 있다.

표준 모듈과 종류 이름의 세부 사항은 GNU Fortran의 ISO_FORTRAN_ENV 설명과 IEEE 모듈 설명에서 확인할 수 있다. 이 장의 문장과 실험 코드는 온도 검사 예제를 위해 새로 구성한 것이다.

완성 코드

프로그램은 외부 입력 없이 실행된다. 첫 번째 실험에서는 두 종류의 실수로 같은 온도 증가를 계산한다. 이어서 상대 허용 오차를 바꾸어 비교하고, 3×3 온도 배열에 시험용 비유한 값을 넣어 상태를 집계한다. NaN과 무한대 자체는 출력하지 않으므로 해당 값의 출력 문자열에 의존하지 않는다.

main.f90

program precision_check
  use iso_fortran_env, only : real32, real64
  use ieee_arithmetic, only : ieee_is_finite, ieee_is_nan, &
       ieee_value, ieee_quiet_nan, ieee_positive_inf, &
       ieee_support_nan, ieee_support_inf
  implicit none

  integer, parameter :: rk = real64
  integer, parameter :: nsteps = 8
  real(real32) :: temp32, delta32
  real(rk) :: temp64, delta64, expected
  real(rk) :: plate(3, 3), bad_nan, bad_inf
  integer :: step

  if (.not. ieee_support_nan(0.0_rk)) then
    error stop 'NaN support is required'
  end if
  if (.not. ieee_support_inf(0.0_rk)) then
    error stop 'Infinity support is required'
  end if

  temp32 = 20.0_real32
  temp64 = 20.0_rk
  delta32 = 2.0_real32 ** (-20)
  delta64 = 2.0_rk ** (-20)
  expected = 20.0_rk + real(nsteps, rk) * delta64

  do step = 1, nsteps
    temp32 = temp32 + delta32
    temp64 = temp64 + delta64
  end do

  write (*, '(A,I0)') 'steps = ', nsteps
  write (*, '(A,F16.10)') 'real32 temperature = ', temp32
  write (*, '(A,F16.10)') 'real64 temperature = ', temp64
  write (*, '(A,F16.10)') 'expected temperature = ', expected
  write (*, '(A,F16.10)') 'absolute difference = ', &
       abs(real(temp32, rk) - temp64)

  write (*, '(A,L1)') 'close, rtol=1e-6 = ', &
       is_close(real(temp32, rk), temp64, 1.0e-8_rk, 1.0e-6_rk)
  write (*, '(A,L1)') 'close, rtol=1e-8 = ', &
       is_close(real(temp32, rk), temp64, 1.0e-8_rk, 1.0e-8_rk)
  write (*, '(A,L1)') 'close near zero = ', &
       is_close(0.0_rk, 5.0e-10_rk, 1.0e-9_rk, 1.0e-6_rk)

  bad_nan = ieee_value(0.0_rk, ieee_quiet_nan)
  bad_inf = ieee_value(0.0_rk, ieee_positive_inf)
  plate = 20.0_rk
  plate(2, 2) = temp64
  plate(1, 3) = bad_nan
  plate(3, 1) = bad_inf

  write (*, '(A,I0)') 'finite cells = ', &
       count(ieee_is_finite(plate))
  write (*, '(A,I0)') 'NaN cells = ', count(ieee_is_nan(plate))
  write (*, '(A,I0)') 'infinite cells = ', &
       count((.not. ieee_is_finite(plate)) .and. &
             (.not. ieee_is_nan(plate)))
  write (*, '(A,L1)') 'NaN accepted = ', &
       is_close(bad_nan, temp64, 1.0e-8_rk, 1.0e-6_rk)
  write (*, '(A,L1)') 'infinity accepted = ', &
       is_close(bad_inf, temp64, 1.0e-8_rk, 1.0e-6_rk)

contains

  logical function is_close(a, b, atol, rtol) result(ok)
    real(rk), intent(in) :: a, b, atol, rtol
    real(rk) :: scale

    ok = .false.
    if (.not. ieee_is_finite(a)) return
    if (.not. ieee_is_finite(b)) return
    if (.not. ieee_is_finite(atol)) return
    if (.not. ieee_is_finite(rtol)) return
    if (atol < 0.0_rk .or. rtol < 0.0_rk) return

    scale = max(abs(a), abs(b))
    ok = abs(a - b) <= max(atol, rtol * scale)
  end function is_close

end program precision_check

줄별 해설

use iso_fortran_env는 실수 종류 이름을 가져온다. only로 필요한 이름만 가져오면 선언의 출처를 확인하기 쉽다. 다음 use ieee_arithmetic에는 검사 함수, 시험값 생성 함수, 시험값의 분류 이름, 지원 여부를 확인하는 함수가 들어 있다.

implicit none 뒤에서 rk를 real64로 정한다. 실제 온도 배열과 비교 함수는 이 종류를 사용한다. temp32와 delta32만 별도로 real32로 선언하여 정밀도 차이를 드러낸다. nsteps는 반복 횟수이며 실행 중 바뀌지 않는다.

두 지원 검사는 시험값을 만들기 전에 수행한다. 0.0_rk는 검사할 실수 종류를 전달하는 역할을 한다. 지원하지 않을 때의 error stop은 프로그램이 이 실험을 수행할 수 없다는 사실을 알린다. 제시한 실행 결과는 두 기능이 지원되는 환경의 정상 경로다.

초기 온도와 변화량은 각각의 종류 접미사를 붙여 지정한다. 2의 음의 정수승을 사용한 이유는 변화량 자체를 이진수로 정확히 표현할 수 있게 하기 위해서다. 따라서 이번 결과를 설명할 때 십진 상수의 변환 오차와 덧셈의 반올림을 섞지 않아도 된다.

expected는 여덟 변화량을 먼저 합한 뒤 기준 온도에 더한다. 이번 값들은 real64에서 정확히 표현되는 이진 값이므로 기대 온도를 명확하게 만들 수 있다. 일반적인 수치 문제에서는 다른 계산식으로 구한 값이 곧 정확한 해라는 보장은 없다. 이 변수는 이번 통제된 실험의 기준값이다.

반복문 안의 두 대입문은 모두 기존 온도에 같은 크기의 변화량을 더한다. temp32는 매번 20으로 반올림되고, temp64에는 증가가 남는다. 반복문 뒤의 형 변환 real(temp32, rk)는 비교할 두 인수의 종류를 맞춘다. 앞서 사라진 증가량을 되살리는 변환은 아니다.

온도와 차이는 F16.10으로 출력한다. 전체 필드 폭은 16칸이고 소수점 아래는 10자리다. 이번 값의 작은 차이를 관찰할 수 있는 폭을 택했다. 출력 형식은 내부 저장 정밀도를 바꾸지 않는다. 폭이 부족할 때 생기는 별표 출력도 값 자체가 NaN이라는 뜻은 아니다.

세 번의 정상 비교는 각각 느슨한 상대 기준, 엄격한 상대 기준, 0 부근의 절대 기준을 확인한다. 논리값은 L1으로 한 글자씩 출력한다. 이후 ieee_value로 만든 두 시험값은 서로 다른 격자점에 저장한다. 중앙 격자점에는 계산한 temp64를 저장한다.

ieee_is_finite(plate)는 3×3 논리 배열을 만든다. count는 그중 참인 원소 수를 센다. 무한대 개수를 계산하는 식은 비유한 원소에서 NaN을 제외한다. 이 식의 두 피연산자는 모두 모든 원소에 안전하게 적용할 수 있는 검사다.

내부 함수 is_close는 결과를 거짓으로 시작한다. 입력과 허용 오차를 순서대로 검사하고 유효하지 않으면 즉시 반환한다. Fortran의 .and.가 뒤쪽 식의 평가를 생략한다고 가정하지 않기 위해, 유한성 검사와 차이 계산을 서로 다른 문장으로 분리했다. 마지막 두 줄에서만 비교 규모와 실제 허용 범위를 계산한다.

실행 결과

파일을 main.f90으로 저장하고 다음 명령을 실행한다. 이 실험은 기본 반올림 설정을 사용하며 -ffast-math 같은 부동소수점 의미를 바꾸는 옵션을 추가하지 않는다. 제시한 코드는 GNU Fortran 16에서 아래 옵션으로 경고 없이 컴파일되도록 구성되어 있다.

gfortran -std=f2018 -Wall main.f90 -o precision_check
./precision_check

예상 출력은 다음과 같다. 실수 필드 앞의 공백도 지정한 서식의 일부다.

steps = 8
real32 temperature =    20.0000000000
real64 temperature =    20.0000076294
expected temperature =    20.0000076294
absolute difference =     0.0000076294
close, rtol=1e-6 = T
close, rtol=1e-8 = F
close near zero = T
finite cells = 7
NaN cells = 1
infinite cells = 1
NaN accepted = F
infinity accepted = F

배열에는 아홉 원소가 있고, 그중 일곱 원소가 유한하다. 중앙의 증가한 온도도 유한 원소에 포함된다. 시험용 NaN과 무한대는 각각 하나씩 검출되고, 둘 다 온도 비교에서 받아들여지지 않는다. 이 출력은 정밀도 실험과 상태 검사가 서로 다른 질문에 답한다는 점을 보여 준다.

실무에서 자주 틀리는 것

넓은 변수에 대입하면 상수도 정밀해진다고 생각한다

다음 코드는 상수를 기본 실수 종류로 해석한 뒤 변환한다. 경고가 없더라도 의도한 정밀도로 상수가 만들어진 것은 아니다.

real(real64) :: increment
increment = 0.1

처음부터 사용할 종류로 상수를 작성한다. 계산식 안의 다른 실수 리터럴에도 같은 규칙을 적용한다.

real(real64) :: increment
increment = 0.1_real64

계산한 실수를 등호로 검증한다

다음 비교는 두 내부 값이 정확히 같은지를 묻는다. 반올림이 포함된 계산 결과가 요구 정확도 안에 있는지 묻는 검사와는 다르다.

if (computed == reference) then
  write (*, '(A)') 'accepted'
end if

완성 코드의 함수를 사용하여 허용 오차를 명시한다. 허용 오차의 선택 이유도 검증 코드나 작업 기록에 남긴다.

if (is_close(computed, reference, 1.0e-8_rk, 1.0e-6_rk)) then
  write (*, '(A)') 'accepted'
end if

실수에 대한 등호 자체가 금지되는 것은 아니다. 정확히 표현되는 상태값을 의도적으로 구별하는 등의 용도가 있다. 문제는 근사 계산의 일치 여부를 등호 하나로 판단하는 데 있다.

차이가 큰 경우만 거부하면 NaN도 잡힌다고 생각한다

다음 코드에서 value가 NaN이면 대소 비교가 거짓이 되어 정상으로 남을 수 있다.

accepted = .true.
if (abs(value - reference) > limit) accepted = .false.

유한성을 먼저 검사한다. 아래 코드는 reference가 유한하고 limit이 유한한 음수 아닌 값이라는 전제에서 단순화한 예다. 완성 코드의 함수는 그 전제도 직접 검사한다.

accepted = .false.
if (ieee_is_finite(value)) then
  accepted = abs(value - reference) <= limit
end if

상대 오차를 기준값 하나로 나누어 계산한다

다음 식은 기준값이 0이면 나눗셈이 유효하지 않으며, 기준값이 매우 작으면 판정이 민감해진다.

accepted = abs((computed - reference) / reference) <= rtol

나눗셈을 없애고 절대 기준을 함께 둔다. 입력과 허용 오차의 검사는 완성 코드처럼 먼저 수행한다.

scale = max(abs(computed), abs(reference))
accepted = abs(computed - reference) <= max(atol, rtol * scale)

한눈에 보기

수치 결과를 점검할 때 구분해야 할 기준

점검 대상	사용할 수단	알 수 있는 것	알 수 없는 것
저장 정밀도	real64, 종류 접미사	값을 어떤 종류로 표현하는지	계산 방법이 충분히 정확한지
표현 간격	epsilon, spacing	작은 변화량의 저장 가능성	전체 계산의 오차 상한
배열 상태	ieee_is_finite, ieee_is_nan	현재 비유한 값이 있는지	그 값이 처음 발생한 원인
결과의 일치	절대·상대 허용 오차	정한 기준 안에 드는지	기준값 자체가 정확한지
화면 표시	F16.10 등의 서식	선택한 자릿수의 결과	내부 값의 정확한 일치

유한하다는 사실은 계산이 정확하다는 뜻이 아니다. 허용 오차 안에 든다는 사실도 물리 모델이나 기준값의 타당성을 대신하지 않는다. 종류 선택, 상태 검사, 수치 비교를 각각 수행해야 결과를 해석할 근거가 생긴다. 다음에 선형 연립방정식을 풀 때에도 해의 저장 정밀도와 계산 결과의 검증을 이 기준으로 나누어 살필 수 있다.

연습 문제

1. 완성 코드에서 두 변화량의 지수를 −20에서 −19로 바꾼다. real32 온도도 증가하는 이유를 설명하고, 최종 온도를 소수점 아래 10자리까지 예상한다.
2. temp32에 변화량을 여덟 번 더하는 대신, 반복문 뒤에서 20.0_real32 + real(nsteps, real32) * delta32를 계산한다. 원래 반복 결과와 다른 이유를 설명한다.
3. is_close에 두 값 0과 5.0e-10_rk를 넣되, atol을 0으로 바꾼다. rtol은 1.0e-6_rk으로 유지한다. 결과를 예상하고 상대 기준만으로 통과하지 못하는 이유를 계산한다.
4. 온도 배열의 plate(1, 1)에도 NaN을 넣는다. 세 가지 개수 출력을 예상한다. 이어서 허용 오차 중 하나가 음수인 경우 완성 코드의 비교 함수가 어떤 값을 반환하는지 설명한다.

정답과 해설

1. 2의 −19승은 이 환경에서 20 부근의 real32 표현 간격과 같다. 한 번 더할 때마다 다음 표현값으로 이동하므로 여덟 증가가 남는다. 전체 증가량은 2의 −16승이며, 두 종류의 최종 온도는 20.0000152588로 출력된다. F16.10에 따른 반올림을 포함한 표시값이다.

2. 원래 지수 −20을 유지하면 여덟 변화량의 합은 2의 −17승이다. 이는 20 부근 간격 네 개에 해당한다. 작은 값을 먼저 합산한 뒤 한 번 더하면 최종 온도는 20.0000076294로 출력된다. 반복 대입에서는 매 단계의 반올림으로 변화가 사라졌지만, 바꾼 식에서는 기준 온도에 더하기 전에 증가량이 충분히 커졌다. 이 사례는 연산 순서가 결과에 영향을 준다는 것을 보여 준다.

3. 결과는 거짓이다. 비교 규모는 5.0e-10_rk이고 상대 허용량은 약 5.0e-16_rk이다. 실제 차이 5.0e-10_rk가 훨씬 크므로 통과하지 못한다. 기존 코드에서는 절대 허용 오차 1.0e-9_rk가 이 차이를 받아들였다.

4. 유한 원소는 6개, NaN은 2개, 무한대는 1개다. 전체 개수는 여전히 9개다. atol이나 rtol 중 하나가 음수이면 함수는 차이 계산 전에 거짓을 반환한다. 호출하는 프로그램에서 비교 실패와 잘못된 허용 오차를 별도로 보고해야 한다면, 허용 오차를 설정할 때 검증하거나 함수에 상태 반환값을 추가할 수 있다.