! Chapter 4: three experiments with floating-point arithmetic.
! Compile: gfortran ch04_roundoff.f95 -o ch04_roundoff
! Run: ./ch04_roundoff
! Use ordinary floating-point settings; do not enable fast-math.
! The extra printed digits expose rounding; they do not add accuracy.

program ch04_roundoff
    implicit none
    integer, parameter :: sp = selected_real_kind(6, 30)
    integer, parameter :: dp = selected_real_kind(15, 300)
    integer :: k
    real(sp) :: one_single, delta_single, single_sum
    real(dp) :: a, b, c, s, double_sum, x, direct, stable
    real(dp) :: temp, first, second

    ! 1. Representation and precision.
    a = 0.1_dp
    b = 0.2_dp
    c = 0.3_dp
    s = a + b
    one_single = 1.0_sp
    delta_single = 1.0e-8_sp
    single_sum = one_single + delta_single
    double_sum = 1.0_dp + 1.0e-8_dp

    write(*,'(a)') '1. Representation and precision'
    write(*,'(a,es24.16)') '0.1 + 0.2:                    ', s
    write(*,'(a,es24.16)') '0.3:                          ', c
    write(*,'(a,es24.16)') 'Difference:                   ', s-c
    write(*,'(a,es24.16)') 'Gap above 1 (single):         ', real(epsilon(1.0_sp),dp)
    write(*,'(a,es24.16)') 'Gap above 1 (double):         ', epsilon(1.0_dp)
    write(*,'(a,es24.16)') '1 + 1e-8 (single):            ', real(single_sum,dp)
    write(*,'(a,es24.16)') '1 + 1e-8 (double):            ', double_sum

    ! 2. Cancellation. The stable formula is a reference, not exact arithmetic.
    write(*,*)
    write(*,'(a)') '2. Cancellation: sqrt(1+x) - 1'
    write(*,'(a12,2a24)') 'x', 'Direct formula', 'Stable formula'
    do k = 2, 18, 2
        x = 10.0_dp**(-k)
        direct = sqrt(1.0_dp+x) - 1.0_dp
        stable = x / (sqrt(1.0_dp+x) + 1.0_dp)
        write(*,'(es12.2,2es24.16)') x, direct, stable
    end do

    ! 3. Order of summation. Store each intermediate result explicitly.
    a = 1.0e16_dp
    b = 1.0_dp
    c = -1.0e16_dp
    temp = a + b
    first = temp + c
    temp = a + c
    second = temp + b

    write(*,*)
    write(*,'(a)') '3. Order of summation'
    write(*,'(a,es24.16)') '(a + b) + c:                  ', first
    write(*,'(a,es24.16)') '(a + c) + b:                  ', second
    write(*,'(a,es24.16)') 'Exact sum:                    ', 1.0_dp
end program ch04_roundoff