module arithmetic_geometric_mean_fortran
    !! Compute arithmetic–geometric mean.

    use, intrinsic :: iso_fortran_env, only: real32
    use, intrinsic :: iso_fortran_env, only: real64
    use, intrinsic :: iso_fortran_env, only: real128

    use, intrinsic :: ieee_arithmetic, only: ieee_quiet_nan
    use, intrinsic :: ieee_arithmetic, only: ieee_unordered
    use, intrinsic :: ieee_arithmetic, only: ieee_value



    implicit none



    private

    public :: arithmetic_geometric_mean
    public :: arithmetic_geometric_mean_kernel
    public :: gap
    public :: max
    public :: min
    public :: n_iter

    public :: arithmetic_geometric_mean_real32_type
    public :: arithmetic_geometric_mean_real64_type
    public :: arithmetic_geometric_mean_real128_type



    integer, parameter :: initial_n_iter = -1

    integer, parameter :: max_n_iter_real32 = digits(0.0_real32)
    integer, parameter :: max_n_iter_real64 = digits(0.0_real64)
    integer, parameter :: max_n_iter_real128 = digits(0.0_real128)




    interface arithmetic_geometric_mean
        !! Safe wrapper for the arithmetic-geometric mean (AGM) computation.
        !!
        !! This interface performs lightweight AGM computation
        !! with input validation, automatic ordering.
        !! Unlike the type-bound subroutine instead,
        !! this interface does not retain intermediate calculation results,
        !! so they cannot be referenced later.
        !!
        !! @note
        !! - If either input was NaN: returns NaN
        !! - If `x` and `y` had opposite signs (`x * y .lt. 0`): returns NaN
        !! - If either `x` or `y` is zero (`x * y .eq. 0`): returns 0
        !! - Otherwise: computes AGM using the iterative kernel
        !! @endnote

        module procedure :: arithmetic_geometric_mean_real32
        module procedure :: arithmetic_geometric_mean_real64
        module procedure :: arithmetic_geometric_mean_real128

    end interface arithmetic_geometric_mean



    interface arithmetic_geometric_mean_kernel
        !! Compute arithmetic–geometric mean (AGM)
        !! using the given arithmetic mean `a` and geometric mean `g`.
        !!
        !! This interface provides a lightweight AGM computation
        !! that returns only the final converged value
        !! without storing iteration history.  
        !! For applications that need to analyze the convergence process,
        !! use the type-bound subroutine instead,
        !! which preserves the full iteration history.
        !!
        !! @warning
        !! - This interface assumes both inputs are positive.
        !! - No validation is performed on inputs.
        !! @endwarning
        !!
        !! @note
        !! **Convergence criterion**  
        !! See [[is_not_converged]]
        !! @endnote

        module procedure :: arithmetic_geometric_mean_kernel_real32
        module procedure :: arithmetic_geometric_mean_kernel_real64
        module procedure :: arithmetic_geometric_mean_kernel_real128

    end interface arithmetic_geometric_mean_kernel



    interface compute_step
        !! Compute arithmetic and geometric mean using the given `prev_a` and `prev_g`.  
        !!
        !! @warning
        !! - This interface assumes both inputs are positive.
        !! - No validation is performed on inputs.
        !! @endwarning


        module procedure :: compute_step_real32
        module procedure :: compute_step_real64
        module procedure :: compute_step_real128

    end interface compute_step



    interface initialize
        !! Initialize components: `n_iter`, `list_a` and `list_g`

        module procedure :: initialize_real32
        module procedure :: initialize_real64
        module procedure :: initialize_real128

    end interface initialize



    interface is_not_converged
        !! Check if the arithmetic-geometric mean iteration has not yet converged.
        !!
        !! @note
        !! **Convergence criterion**  
        !! Returns `.true.` if the absolute difference between the arithmetic mean `a`
        !! and geometric mean `g` exceeds the machine epsilon relative to the smaller value,
        !! as determined by `spacing(min(a, g))`.
        !!
        !! Mathematically: `|a - g| > spacing(min(a, g))`
        !!
        !! **Appendix**  
        !! This interface is designed for internal use within the AGM iteration where
        !! both values are guaranteed to be positive and converging.
        !! @endnote

        module procedure :: is_not_converged_real32
        module procedure :: is_not_converged_real64
        module procedure :: is_not_converged_real128

    end interface is_not_converged



    interface gap

        module procedure :: gap_final_real32
        module procedure :: gap_final_real64
        module procedure :: gap_final_real128

        module procedure :: gap_selectable_real32
        module procedure :: gap_selectable_real64
        module procedure :: gap_selectable_real128

    end interface gap



    interface max

        module procedure :: max_final_real32
        module procedure :: max_final_real64
        module procedure :: max_final_real128

        module procedure :: max_selectable_real32
        module procedure :: max_selectable_real64
        module procedure :: max_selectable_real128

    end interface max



    interface min

        module procedure :: min_final_real32
        module procedure :: min_final_real64
        module procedure :: min_final_real128

        module procedure :: min_selectable_real32
        module procedure :: min_selectable_real64
        module procedure :: min_selectable_real128

    end interface min



    type, abstract :: arithmetic_geometric_mean_base_type

        integer, private :: n_iter_ = initial_n_iter
        !! the number of iterations performed during AGM calculation

        contains

        procedure, pass, public :: n_iter

    end type arithmetic_geometric_mean_base_type



    type, extends(arithmetic_geometric_mean_base_type) :: arithmetic_geometric_mean_real32_type

        real(real32), private :: list_a(0:max_n_iter_real32)
        !! history of the arithmetic mean

        real(real32), private :: list_g(0:max_n_iter_real32)
        !! history of the geometric mean

        contains

        procedure, pass, private :: compute_real32
        procedure, pass, private :: compute_kernel_real32

        generic, public  :: compute        => compute_real32
        generic, private :: compute_kernel => compute_kernel_real32

    end type arithmetic_geometric_mean_real32_type



    type, extends(arithmetic_geometric_mean_base_type) :: arithmetic_geometric_mean_real64_type

        real(real64), private :: list_a(0:max_n_iter_real64)
        !! history of the arithmetic mean

        real(real64), private :: list_g(0:max_n_iter_real64)
        !! history of the geometric mean

        contains

        procedure, pass, private :: compute_real64
        procedure, pass, private :: compute_kernel_real64

        generic, public  :: compute        => compute_real64
        generic, private :: compute_kernel => compute_kernel_real64

    end type arithmetic_geometric_mean_real64_type



    type, extends(arithmetic_geometric_mean_base_type) :: arithmetic_geometric_mean_real128_type

        real(real128), private :: list_a(0:max_n_iter_real128)
        !! history of the arithmetic mean

        real(real128), private :: list_g(0:max_n_iter_real128)
        !! history of the geometric mean

        contains

        procedure, pass, private :: compute_real128
        procedure, pass, private :: compute_kernel_real128

        generic, public  :: compute        => compute_real128
        generic, private :: compute_kernel => compute_kernel_real128

    end type arithmetic_geometric_mean_real128_type



    contains



    elemental function arithmetic_geometric_mean_real32(x, y) result(agm)
        !! Safe wrapper for the arithmetic-geometric mean (AGM) computation.
        !!
        !! This function performs lightweight AGM computation
        !! with input validation, automatic ordering.
        !! Unlike the type-bound subroutine [[compute_real32]] instead,
        !! this function does not retain intermediate calculation results,
        !! so they cannot be referenced later.
        !!
        !! @note
        !! - If either input was NaN: returns NaN
        !! - If `x` and `y` had opposite signs (`x * y .lt. 0`): returns NaN
        !! - If either `x` or `y` is zero (`x * y .eq. 0`): returns 0
        !! - Otherwise: computes AGM using the iterative kernel
        !! @endnote

        real(real32), intent(in) :: x, y



        real(real32) :: agm ! return value



        real(real32) :: xy !! x * y



        if ( ieee_unordered(x, y) ) then

            agm = ieee_value(agm, ieee_quiet_nan); return

        end if



        xy = x * y



        if ( xy .lt. 0.0_real32 ) then

            agm = ieee_value(agm, ieee_quiet_nan)

        else if ( xy .gt. 0.0_real32 ) then

            if (x .lt. y) then
                agm = arithmetic_geometric_mean_kernel( a = y, g = x )
            else
                agm = arithmetic_geometric_mean_kernel( a = x, g = y )
            end if

        else

            agm = 0.0_real32

        end if

    end function arithmetic_geometric_mean_real32



    elemental function arithmetic_geometric_mean_real64(x, y) result(agm)
        !! Safe wrapper for the arithmetic-geometric mean (AGM) computation.
        !!
        !! This function performs lightweight AGM computation
        !! with input validation, automatic ordering.
        !! Unlike the type-bound subroutine [[compute_real64]] instead,
        !! this function does not retain intermediate calculation results,
        !! so they cannot be referenced later.
        !!
        !! @note
        !! - If either input was NaN: returns NaN
        !! - If `x` and `y` had opposite signs (`x * y .lt. 0`): returns NaN
        !! - If either `x` or `y` is zero (`x * y .eq. 0`): returns 0
        !! - Otherwise: computes AGM using the iterative kernel
        !! @endnote

        real(real64), intent(in) :: x, y



        real(real64) :: agm ! return value



        real(real64) :: xy !! x * y



        if ( ieee_unordered(x, y) ) then

            agm = ieee_value(agm, ieee_quiet_nan); return

        end if



        xy = x * y



        if ( xy .lt. 0.0_real64 ) then

            agm = ieee_value(agm, ieee_quiet_nan)

        else if ( xy .gt. 0.0_real64 ) then

            if (x .lt. y) then
                agm = arithmetic_geometric_mean_kernel( a = y, g = x )
            else
                agm = arithmetic_geometric_mean_kernel( a = x, g = y )
            end if

        else

            agm = 0.0_real64

        end if

    end function arithmetic_geometric_mean_real64



    elemental function arithmetic_geometric_mean_real128(x, y) result(agm)
        !! Safe wrapper for the arithmetic-geometric mean (AGM) computation.
        !!
        !! This function performs lightweight AGM computation
        !! with input validation, automatic ordering.
        !! Unlike the type-bound subroutine [[compute_real128]] instead,
        !! this function does not retain intermediate calculation results,
        !! so they cannot be referenced later.
        !!
        !! @note
        !! - If either input was NaN: returns NaN
        !! - If `x` and `y` had opposite signs (`x * y .lt. 0`): returns NaN
        !! - If either `x` or `y` is zero (`x * y .eq. 0`): returns 0
        !! - Otherwise: computes AGM using the iterative kernel
        !! @endnote

        real(real128), intent(in) :: x, y



        real(real128) :: agm ! return value



        real(real128) :: xy !! x * y



        if ( ieee_unordered(x, y) ) then

            agm = ieee_value(agm, ieee_quiet_nan); return

        end if



        xy = x * y



        if ( xy .lt. 0.0_real128 ) then

            agm = ieee_value(agm, ieee_quiet_nan)

        else if ( xy .gt. 0.0_real128 ) then

            if (x .lt. y) then
                agm = arithmetic_geometric_mean_kernel( a = y, g = x )
            else
                agm = arithmetic_geometric_mean_kernel( a = x, g = y )
            end if

        else

            agm = 0.0_real128

        end if

    end function arithmetic_geometric_mean_real128



    elemental function arithmetic_geometric_mean_kernel_real32(a, g) result(agm)
        !! Compute arithmetic–geometric mean (AGM)
        !! using the given arithmetic mean `a` and geometric mean `g`.
        !!
        !! This function provides a lightweight AGM computation
        !! that returns only the final converged value
        !! without storing iteration history.  
        !! For applications that need to analyze the convergence process,
        !! use the type-bound subroutine [[compute_kernel_real32]] instead,
        !! which preserves the full iteration history.
        !!
        !! @warning
        !! - This function assumes both inputs are positive.
        !! - No validation is performed on inputs.
        !! @endwarning
        !!
        !! @note
        !! **Convergence criterion**  
        !! See [[is_not_converged]]
        !! @endnote

        real(real32), intent(in) :: a !! arithmetic mean

        real(real32), intent(in) :: g !! geometric mean



        real(real32) :: agm ! return value



        real(real32) :: prev_a !! previous arithmetic mean

        real(real32) :: prev_g !! previous geometric mean

        real(real32) :: next_a !! next arithmetic mean

        real(real32) :: next_g !! next geometric mean



        prev_a = a
        prev_g = g



        do

            call compute_step( &!
                prev_a = prev_a , &!
                prev_g = prev_g , &!
                next_a = next_a , &!
                next_g = next_g   &!
            )

            if ( is_not_converged(next_a, next_g) ) then

                prev_a = next_a
                prev_g = next_g

                cycle

            else

                agm = max(next_a, next_g)

                return

            end if

        end do

    end function arithmetic_geometric_mean_kernel_real32



    elemental function arithmetic_geometric_mean_kernel_real64(a, g) result(agm)
        !! Compute arithmetic–geometric mean (AGM)
        !! using the given arithmetic mean `a` and geometric mean `g`.
        !!
        !! This function provides a lightweight AGM computation
        !! that returns only the final converged value
        !! without storing iteration history.  
        !! For applications that need to analyze the convergence process,
        !! use the type-bound subroutine [[compute_kernel_real64]] instead,
        !! which preserves the full iteration history.
        !!
        !! @warning
        !! - This function assumes both inputs are positive.
        !! - No validation is performed on inputs.
        !! @endwarning
        !!
        !! @note
        !! **Convergence criterion**  
        !! See [[is_not_converged]]
        !! @endnote

        real(real64), intent(in) :: a !! arithmetic mean

        real(real64), intent(in) :: g !! geometric mean



        real(real64) :: agm ! return value



        real(real64) :: prev_a !! previous arithmetic mean

        real(real64) :: prev_g !! previous geometric mean

        real(real64) :: next_a !! next arithmetic mean

        real(real64) :: next_g !! next geometric mean



        prev_a = a
        prev_g = g



        do

            call compute_step( &!
                prev_a = prev_a , &!
                prev_g = prev_g , &!
                next_a = next_a , &!
                next_g = next_g   &!
            )

            if ( is_not_converged(next_a, next_g) ) then

                prev_a = next_a
                prev_g = next_g

                cycle

            else

                agm = max(next_a, next_g)

                return

            end if

        end do

    end function arithmetic_geometric_mean_kernel_real64



    elemental function arithmetic_geometric_mean_kernel_real128(a, g) result(agm)
        !! Compute arithmetic–geometric mean (AGM)
        !! using the given arithmetic mean `a` and geometric mean `g`.
        !!
        !! This function provides a lightweight AGM computation
        !! that returns only the final converged value
        !! without storing iteration history.  
        !! For applications that need to analyze the convergence process,
        !! use the type-bound subroutine [[compute_kernel_real128]] instead,
        !! which preserves the full iteration history.
        !!
        !! @warning
        !! - This function assumes both inputs are positive.
        !! - No validation is performed on inputs.
        !! @endwarning
        !!
        !! @note
        !! **Convergence criterion**  
        !! See [[is_not_converged]]
        !! @endnote

        real(real128), intent(in) :: a !! arithmetic mean

        real(real128), intent(in) :: g !! geometric mean



        real(real128) :: agm ! return value



        real(real128) :: prev_a !! previous arithmetic mean

        real(real128) :: prev_g !! previous geometric mean

        real(real128) :: next_a !! next arithmetic mean

        real(real128) :: next_g !! next geometric mean



        prev_a = a
        prev_g = g



        do

            call compute_step( &!
                prev_a = prev_a , &!
                prev_g = prev_g , &!
                next_a = next_a , &!
                next_g = next_g   &!
            )

            if ( is_not_converged(next_a, next_g) ) then

                prev_a = next_a
                prev_g = next_g

                cycle

            else

                agm = max(next_a, next_g)

                return

            end if

        end do

    end function arithmetic_geometric_mean_kernel_real128



    elemental function gap_final_real32(agm) result(gap_final)
        !! Extract the final gap (arithmetic mean - geometric mean).
        !!
        !! @warning
        !! This function assumes the AGM computation has been performed via the `compute` method.  
        !! If called on an uninitialized or improperly computed AGM object,
        !! the result may be NaN or undefined.
        !! @endwarning

        type(arithmetic_geometric_mean_real32_type), intent(in) :: agm



        real(real32) :: gap_final



        gap_final = gap(agm, agm%n_iter_)

    end function gap_final_real32



    elemental function gap_final_real64(agm) result(gap_final)
        !! Extract the final gap (arithmetic mean - geometric mean).
        !!
        !! @warning
        !! This function assumes the AGM computation has been performed via the `compute` method.  
        !! If called on an uninitialized or improperly computed AGM object,
        !! the result may be NaN or undefined.
        !! @endwarning

        type(arithmetic_geometric_mean_real64_type), intent(in) :: agm



        real(real64) :: gap_final



        gap_final = gap(agm, agm%n_iter_)

    end function gap_final_real64



    elemental function gap_final_real128(agm) result(gap_final)
        !! Extract the final gap (arithmetic mean - geometric mean).
        !!
        !! @warning
        !! This function assumes the AGM computation has been performed via the `compute` method.  
        !! If called on an uninitialized or improperly computed AGM object,
        !! the result may be NaN or undefined.
        !! @endwarning

        type(arithmetic_geometric_mean_real128_type), intent(in) :: agm



        real(real128) :: gap_final



        gap_final = gap(agm, agm%n_iter_)

    end function gap_final_real128



    elemental function gap_selectable_real32(agm, i) result(gap_selectable)
        !! Extract the gap (arithmetic mean - geometric mean) at a specific iteration.
        !!
        !! @warning
        !! This function assumes the AGM computation has been performed via the `compute` method.  
        !! If called on an uninitialized or improperly computed AGM object,
        !! the result may be NaN or undefined.
        !!
        !! The iteration index `i` must be valid: `0 <= i <= n_iter`.  
        !! No bounds checking is performed; invalid indices may cause undefined behavior.
        !! @endwarning

        type(arithmetic_geometric_mean_real32_type), intent(in) :: agm

        integer, intent(in) :: i



        real(real32) :: gap_selectable



        gap_selectable = agm%list_a(i) - agm%list_g(i)

    end function gap_selectable_real32



    elemental function gap_selectable_real64(agm, i) result(gap_selectable)
        !! Extract the gap (arithmetic mean - geometric mean) at a specific iteration.
        !!
        !! @warning
        !! This function assumes the AGM computation has been performed via the `compute` method.  
        !! If called on an uninitialized or improperly computed AGM object,
        !! the result may be NaN or undefined.
        !!
        !! The iteration index `i` must be valid: `0 <= i <= n_iter`.  
        !! No bounds checking is performed; invalid indices may cause undefined behavior.
        !! @endwarning

        type(arithmetic_geometric_mean_real64_type), intent(in) :: agm

        integer, intent(in) :: i



        real(real64) :: gap_selectable



        gap_selectable = agm%list_a(i) - agm%list_g(i)

    end function gap_selectable_real64



    elemental function gap_selectable_real128(agm, i) result(gap_selectable)
        !! Extract the gap (arithmetic mean - geometric mean) at a specific iteration.
        !!
        !! @warning
        !! This function assumes the AGM computation has been performed via the `compute` method.  
        !! If called on an uninitialized or improperly computed AGM object,
        !! the result may be NaN or undefined.
        !!
        !! The iteration index `i` must be valid: `0 <= i <= n_iter`.  
        !! No bounds checking is performed; invalid indices may cause undefined behavior.
        !! @endwarning

        type(arithmetic_geometric_mean_real128_type), intent(in) :: agm

        integer, intent(in) :: i



        real(real128) :: gap_selectable



        gap_selectable = agm%list_a(i) - agm%list_g(i)

    end function gap_selectable_real128



    elemental function is_not_converged_real32(a, g) result(stat)
        !! Check if the arithmetic-geometric mean iteration has not yet converged.
        !!
        !! @note
        !! **Convergence criterion**  
        !! Returns `.true.` if the absolute difference between the arithmetic mean `a`
        !! and geometric mean `g` exceeds the machine epsilon relative to the smaller value,
        !! as determined by `spacing(min(a, g))`.
        !!
        !! Mathematically: `|a - g| > spacing(min(a, g))`
        !!
        !! **Appendix**  
        !! This function is designed for internal use within the AGM iteration where
        !! both values are guaranteed to be positive and converging.
        !! @endnote

        real(real32), intent(in) :: a !! arithmetic mean

        real(real32), intent(in) :: g !! geometric mean



        logical :: stat



        stat = abs(a - g) .gt. spacing( min(a, g) )

    end function is_not_converged_real32



    elemental function is_not_converged_real64(a, g) result(stat)
        !! Check if the arithmetic-geometric mean iteration has not yet converged.
        !!
        !! @note
        !! **Convergence criterion**  
        !! Returns `.true.` if the absolute difference between the arithmetic mean `a`
        !! and geometric mean `g` exceeds the machine epsilon relative to the smaller value,
        !! as determined by `spacing(min(a, g))`.
        !!
        !! Mathematically: `|a - g| > spacing(min(a, g))`
        !!
        !! **Appendix**  
        !! This function is designed for internal use within the AGM iteration where
        !! both values are guaranteed to be positive and converging.
        !! @endnote

        real(real64), intent(in) :: a !! arithmetic mean

        real(real64), intent(in) :: g !! geometric mean



        logical :: stat



        stat = abs(a - g) .gt. spacing( min(a, g) )

    end function is_not_converged_real64



    elemental function is_not_converged_real128(a, g) result(stat)
        !! Check if the arithmetic-geometric mean iteration has not yet converged.
        !!
        !! @note
        !! **Convergence criterion**  
        !! Returns `.true.` if the absolute difference between the arithmetic mean `a`
        !! and geometric mean `g` exceeds the machine epsilon relative to the smaller value,
        !! as determined by `spacing(min(a, g))`.
        !!
        !! Mathematically: `|a - g| > spacing(min(a, g))`
        !!
        !! **Appendix**  
        !! This function is designed for internal use within the AGM iteration where
        !! both values are guaranteed to be positive and converging.
        !! @endnote

        real(real128), intent(in) :: a !! arithmetic mean

        real(real128), intent(in) :: g !! geometric mean



        logical :: stat



        stat = abs(a - g) .gt. spacing( min(a, g) )

    end function is_not_converged_real128



    elemental function max_final_real32(agm) result(max_final)
        !! Extract the final arithmetic-geometric mean from completed AGM calculations.  
        !! Specifically, this function returns the larger of the two final values.
        !!
        !! @warning
        !! This function assumes the AGM computation has been performed via the `compute` method.  
        !! If called on an uninitialized or improperly computed AGM object,
        !! the result may be NaN or undefined.
        !! @endwarning

        type(arithmetic_geometric_mean_real32_type), intent(in) :: agm



        real(real32) :: max_final



        max_final = max(agm, agm%n_iter_)

    end function max_final_real32



    elemental function max_final_real64(agm) result(max_final)
        !! Extract the final arithmetic-geometric mean from completed AGM calculations.  
        !! Specifically, this function returns the larger of the two final values.
        !!
        !! @warning
        !! This function assumes the AGM computation has been performed via the `compute` method.  
        !! If called on an uninitialized or improperly computed AGM object,
        !! the result may be NaN or undefined.
        !! @endwarning

        type(arithmetic_geometric_mean_real64_type), intent(in) :: agm



        real(real64) :: max_final



        max_final = max(agm, agm%n_iter_)

    end function max_final_real64



    elemental function max_final_real128(agm) result(max_final)
        !! Extract the final arithmetic-geometric mean from completed AGM calculations.  
        !! Specifically, this function returns the larger of the two final values.
        !!
        !! @warning
        !! This function assumes the AGM computation has been performed via the `compute` method.  
        !! If called on an uninitialized or improperly computed AGM object,
        !! the result may be NaN or undefined.
        !! @endwarning

        type(arithmetic_geometric_mean_real128_type), intent(in) :: agm



        real(real128) :: max_final



        max_final = max(agm, agm%n_iter_)

    end function max_final_real128



    elemental function max_selectable_real32(agm, i) result(max_selectable)
        !! Extract the arithmetic-geometric mean value at a specific iteration.  
        !! Specifically, this function returns the larger of the two values.
        !!
        !! @warning
        !! This function assumes the AGM computation has been performed via the `compute` method.  
        !! If called on an uninitialized or improperly computed AGM object,
        !! the result may be NaN or undefined.
        !!
        !! The iteration index `i` must be valid: `0 <= i <= n_iter`.  
        !! No bounds checking is performed; invalid indices may cause undefined behavior.
        !! @endwarning

        type(arithmetic_geometric_mean_real32_type), intent(in) :: agm

        integer, intent(in) :: i



        real(real32) :: max_selectable



        max_selectable = max( agm%list_a(i), agm%list_g(i) )

    end function max_selectable_real32



    elemental function max_selectable_real64(agm, i) result(max_selectable)
        !! Extract the arithmetic-geometric mean value at a specific iteration.  
        !! Specifically, this function returns the larger of the two values.
        !!
        !! @warning
        !! This function assumes the AGM computation has been performed via the `compute` method.  
        !! If called on an uninitialized or improperly computed AGM object,
        !! the result may be NaN or undefined.
        !!
        !! The iteration index `i` must be valid: `0 <= i <= n_iter`.  
        !! No bounds checking is performed; invalid indices may cause undefined behavior.
        !! @endwarning

        type(arithmetic_geometric_mean_real64_type), intent(in) :: agm

        integer, intent(in) :: i



        real(real64) :: max_selectable



        max_selectable = max( agm%list_a(i), agm%list_g(i) )

    end function max_selectable_real64



    elemental function max_selectable_real128(agm, i) result(max_selectable)
        !! Extract the arithmetic-geometric mean value at a specific iteration.  
        !! Specifically, this function returns the larger of the two values.
        !!
        !! @warning
        !! This function assumes the AGM computation has been performed via the `compute` method.  
        !! If called on an uninitialized or improperly computed AGM object,
        !! the result may be NaN or undefined.
        !!
        !! The iteration index `i` must be valid: `0 <= i <= n_iter`.  
        !! No bounds checking is performed; invalid indices may cause undefined behavior.
        !! @endwarning

        type(arithmetic_geometric_mean_real128_type), intent(in) :: agm

        integer, intent(in) :: i



        real(real128) :: max_selectable



        max_selectable = max( agm%list_a(i), agm%list_g(i) )

    end function max_selectable_real128



    elemental function min_final_real32(agm) result(min_final)
        !! Extract the final arithmetic-geometric mean from completed AGM calculations.  
        !! Specifically, this function returns the smaller of the two final values.
        !!
        !! @warning
        !! This function assumes the AGM computation has been performed via the `compute` method.  
        !! If called on an uninitialized or improperly computed AGM object,
        !! the result may be NaN or undefined.
        !! @endwarning

        type(arithmetic_geometric_mean_real32_type), intent(in) :: agm



        real(real32) :: min_final



        min_final = min(agm, agm%n_iter_)

    end function min_final_real32



    elemental function min_final_real64(agm) result(min_final)
        !! Extract the final arithmetic-geometric mean from completed AGM calculations.  
        !! Specifically, this function returns the smaller of the two final values.
        !!
        !! @warning
        !! This function assumes the AGM computation has been performed via the `compute` method.  
        !! If called on an uninitialized or improperly computed AGM object,
        !! the result may be NaN or undefined.
        !! @endwarning

        type(arithmetic_geometric_mean_real64_type), intent(in) :: agm



        real(real64) :: min_final



        min_final = min(agm, agm%n_iter_)

    end function min_final_real64



    elemental function min_final_real128(agm) result(min_final)
        !! Extract the final arithmetic-geometric mean from completed AGM calculations.  
        !! Specifically, this function returns the smaller of the two final values.
        !!
        !! @warning
        !! This function assumes the AGM computation has been performed via the `compute` method.  
        !! If called on an uninitialized or improperly computed AGM object,
        !! the result may be NaN or undefined.
        !! @endwarning

        type(arithmetic_geometric_mean_real128_type), intent(in) :: agm



        real(real128) :: min_final



        min_final = min(agm, agm%n_iter_)

    end function min_final_real128



    elemental function min_selectable_real32(agm, i) result(min_selectable)
        !! Extract the arithmetic-geometric mean value at a specific iteration.  
        !! Specifically, this function returns the smaller of the two values.
        !!
        !! @warning
        !! This function assumes the AGM computation has been performed via the `compute` method.  
        !! If called on an uninitialized or improperly computed AGM object,
        !! the result may be NaN or undefined.
        !!
        !! The iteration index `i` must be valid: `0 <= i <= n_iter`.  
        !! No bounds checking is performed; invalid indices may cause undefined behavior.
        !! @endwarning

        type(arithmetic_geometric_mean_real32_type), intent(in) :: agm

        integer, intent(in) :: i



        real(real32) :: min_selectable



        min_selectable = min( agm%list_a(i), agm%list_g(i) )

    end function min_selectable_real32



    elemental function min_selectable_real64(agm, i) result(min_selectable)
        !! Extract the arithmetic-geometric mean value at a specific iteration.  
        !! Specifically, this function returns the smaller of the two values.
        !!
        !! @warning
        !! This function assumes the AGM computation has been performed via the `compute` method.  
        !! If called on an uninitialized or improperly computed AGM object,
        !! the result may be NaN or undefined.
        !!
        !! The iteration index `i` must be valid: `0 <= i <= n_iter`.  
        !! No bounds checking is performed; invalid indices may cause undefined behavior.
        !! @endwarning

        type(arithmetic_geometric_mean_real64_type), intent(in) :: agm

        integer, intent(in) :: i



        real(real64) :: min_selectable



        min_selectable = min( agm%list_a(i), agm%list_g(i) )

    end function min_selectable_real64



    elemental function min_selectable_real128(agm, i) result(min_selectable)
        !! Extract the arithmetic-geometric mean value at a specific iteration.  
        !! Specifically, this function returns the smaller of the two values.
        !!
        !! @warning
        !! This function assumes the AGM computation has been performed via the `compute` method.  
        !! If called on an uninitialized or improperly computed AGM object,
        !! the result may be NaN or undefined.
        !!
        !! The iteration index `i` must be valid: `0 <= i <= n_iter`.  
        !! No bounds checking is performed; invalid indices may cause undefined behavior.
        !! @endwarning

        type(arithmetic_geometric_mean_real128_type), intent(in) :: agm

        integer, intent(in) :: i



        real(real128) :: min_selectable



        min_selectable = min( agm%list_a(i), agm%list_g(i) )

    end function min_selectable_real128


    elemental function n_iter(agm)
        !! return the number of iterations performed during AGM calculation

        class(arithmetic_geometric_mean_base_type), intent(in) :: agm



        integer :: n_iter



         n_iter = agm%n_iter_

    end function n_iter



    elemental subroutine compute_real32(agm, x, y)
        !! Safe wrapper for the arithmetic-geometric mean (AGM) computation.
        !!
        !! This type-bound subroutine performs AGM computation
        !! with input validation, automatic ordering,
        !! and complete convergence history storage.
        !! Unlike [[arithmetic_geometric_mean_kernel]],
        !! this subroutine retains all intermediate calculation results
        !! and can be retrieved later.
        !!
        !! @note
        !! - If either input was NaN: returns NaN
        !! - If `x` and `y` had opposite signs (`x * y .lt. 0`): returns NaN
        !! - If either `x` or `y` is zero (`x * y .eq. 0`): returns 0
        !! - Otherwise: full AGM iteration using [[compute_kernel_real32]]
        !! @endnote
        !!
        !! @warning
        !! Calling this subroutine overwrites any previous computation
        !! stored in the object.
        !! @endwarning

        class(arithmetic_geometric_mean_real32_type), intent(inout) :: agm

        real(real32), intent(in) :: x, y



        real(real32) :: xy



        if ( ieee_unordered(x, y) ) then

            call initialize(agm)

            agm%n_iter_ = agm%n_iter_ + 1

            return

        end if



        xy = x * y



        if ( xy .lt. 0.0_real32 ) then

            call initialize(agm)

            agm%n_iter_ = agm%n_iter_ + 1

        else if ( xy .gt. 0.0_real32 ) then

            if (x .lt. y) then
                call agm%compute_kernel(init_a = y, init_g = x)
            else
                call agm%compute_kernel(init_a = x, init_g = y)
            end if

        else

            call initialize(agm)

            agm%n_iter_ = 0

            if (x .lt. y) then

                agm%list_a(agm%n_iter_) = y
                agm%list_g(agm%n_iter_) = x

            else

                agm%list_a(agm%n_iter_) = x
                agm%list_g(agm%n_iter_) = y
        
            end if

            agm%n_iter_ = agm%n_iter_ + 1

            agm%list_a(agm%n_iter_) = 0.0_real32
            agm%list_g(agm%n_iter_) = 0.0_real32

        end if

    end subroutine compute_real32



    elemental subroutine compute_real64(agm, x, y)
        !! Safe wrapper for the arithmetic-geometric mean (AGM) computation.
        !!
        !! This type-bound subroutine performs AGM computation
        !! with input validation, automatic ordering,
        !! and complete convergence history storage.
        !! Unlike [[arithmetic_geometric_mean_kernel]],
        !! this subroutine retains all intermediate calculation results
        !! and can be retrieved later.
        !!
        !! @note
        !! - If either input was NaN: returns NaN
        !! - If `x` and `y` had opposite signs (`x * y .lt. 0`): returns NaN
        !! - If either `x` or `y` is zero (`x * y .eq. 0`): returns 0
        !! - Otherwise: full AGM iteration using [[compute_kernel_real64]]
        !! @endnote
        !!
        !! @warning
        !! Calling this subroutine overwrites any previous computation
        !! stored in the object.
        !! @endwarning

        class(arithmetic_geometric_mean_real64_type), intent(inout) :: agm

        real(real64), intent(in) :: x, y



        real(real64) :: xy



        if ( ieee_unordered(x, y) ) then

            call initialize(agm)

            agm%n_iter_ = agm%n_iter_ + 1

            return

        end if



        xy = x * y



        if ( xy .lt. 0.0_real64 ) then

            call initialize(agm)

            agm%n_iter_ = agm%n_iter_ + 1

        else if ( xy .gt. 0.0_real64 ) then

            if (x .lt. y) then
                call agm%compute_kernel(init_a = y, init_g = x)
            else
                call agm%compute_kernel(init_a = x, init_g = y)
            end if

        else

            call initialize(agm)

            agm%n_iter_ = 0

            if (x .lt. y) then

                agm%list_a(agm%n_iter_) = y
                agm%list_g(agm%n_iter_) = x

            else

                agm%list_a(agm%n_iter_) = x
                agm%list_g(agm%n_iter_) = y
        
            end if

            agm%n_iter_ = agm%n_iter_ + 1

            agm%list_a(agm%n_iter_) = 0.0_real64
            agm%list_g(agm%n_iter_) = 0.0_real64

        end if

    end subroutine compute_real64



    elemental subroutine compute_real128(agm, x, y)
        !! Safe wrapper for the arithmetic-geometric mean (AGM) computation.
        !!
        !! This type-bound subroutine performs AGM computation
        !! with input validation, automatic ordering,
        !! and complete convergence history storage.
        !! Unlike [[arithmetic_geometric_mean_kernel]],
        !! this subroutine retains all intermediate calculation results
        !! and can be retrieved later.
        !!
        !! @note
        !! - If either input was NaN: returns NaN
        !! - If `x` and `y` had opposite signs (`x * y .lt. 0`): returns NaN
        !! - If either `x` or `y` is zero (`x * y .eq. 0`): returns 0
        !! - Otherwise: full AGM iteration using [[compute_kernel_real128]]
        !! @endnote
        !!
        !! @warning
        !! Calling this subroutine overwrites any previous computation
        !! stored in the object.
        !! @endwarning

        class(arithmetic_geometric_mean_real128_type), intent(inout) :: agm

        real(real128), intent(in) :: x, y



        real(real128) :: xy



        if ( ieee_unordered(x, y) ) then

            call initialize(agm)

            agm%n_iter_ = agm%n_iter_ + 1

            return

        end if



        xy = x * y



        if ( xy .lt. 0.0_real128 ) then

            call initialize(agm)

            agm%n_iter_ = agm%n_iter_ + 1

        else if ( xy .gt. 0.0_real128 ) then

            if (x .lt. y) then
                call agm%compute_kernel(init_a = y, init_g = x)
            else
                call agm%compute_kernel(init_a = x, init_g = y)
            end if

        else

            call initialize(agm)

            agm%n_iter_ = 0

            if (x .lt. y) then

                agm%list_a(agm%n_iter_) = y
                agm%list_g(agm%n_iter_) = x

            else

                agm%list_a(agm%n_iter_) = x
                agm%list_g(agm%n_iter_) = y
        
            end if

            agm%n_iter_ = agm%n_iter_ + 1

            agm%list_a(agm%n_iter_) = 0.0_real128
            agm%list_g(agm%n_iter_) = 0.0_real128

        end if

    end subroutine compute_real128



    elemental subroutine compute_kernel_real32(agm, init_a, init_g)
        !! Compute arithmetic–geometric mean (AGM)
        !! using the given arithmetic mean `init_a` and geometric mean `init_g`
        !! and store the complete iteration history.
        !!
        !! This type-bound subroutine performs the AGM iteration starting from
        !! given initial values, storing each intermediate arithmetic and geometric mean
        !! in `list_a` and `list_g` respectively.
        !! The final iteration count is stored in `n_iter`.
        !!
        !! Unlike the [[arithmetic_geometric_mean_kernel]],
        !! this subroutine preserves the full convergence history for analysis.
        !!
        !! @warning
        !! - This subroutine/interface assumes both inputs are positive.
        !! - No validation is performed on inputs.
        !! - Calling this subroutine overwrites any previous computation stored in the object.
        !! @endwarning
        !!
        !! @note
        !! **Convergence criterion**  
        !! See [[is_not_converged]]
        !! @endnote

        class(arithmetic_geometric_mean_real32_type), intent(inout) :: agm

        real(real32), intent(in) :: init_a
        !! initial value: arithmetic mean

        real(real32), intent(in) :: init_g
        !! initial value: geometric mean



        call initialize(agm)

        agm%list_a(0) = init_a
        agm%list_g(0) = init_g
        agm%n_iter_    = 0

        do

            associate(prev_iter => agm%n_iter_, next_iter => agm%n_iter_ + 1)

                associate( &!
                    prev_a => agm%list_a(prev_iter) , &!
                    prev_g => agm%list_g(prev_iter) , &!
                    next_a => agm%list_a(next_iter) , &!
                    next_g => agm%list_g(next_iter)   &!
                )

                    agm%n_iter_ = agm%n_iter_ + 1

                    call compute_step( &!
                        prev_a = prev_a , &!
                        prev_g = prev_g , &!
                        next_a = next_a , &!
                        next_g = next_g   &!
                    )

                    if ( is_not_converged(next_a, next_g) ) then

                        prev_iter = next_iter

                        cycle

                    else

                        prev_iter = next_iter

                        return

                    end if

                end associate

            end associate

        end do

    end subroutine compute_kernel_real32



    elemental subroutine compute_kernel_real64(agm, init_a, init_g)
        !! Compute arithmetic–geometric mean (AGM)
        !! using the given arithmetic mean `init_a` and geometric mean `init_g`
        !! and store the complete iteration history.
        !!
        !! This type-bound subroutine performs the AGM iteration starting from
        !! given initial values, storing each intermediate arithmetic and geometric mean
        !! in `list_a` and `list_g` respectively.
        !! The final iteration count is stored in `n_iter`.
        !!
        !! Unlike the [[arithmetic_geometric_mean_kernel]],
        !! this subroutine preserves the full convergence history for analysis.
        !!
        !! @warning
        !! - This subroutine/interface assumes both inputs are positive.
        !! - No validation is performed on inputs.
        !! - Calling this subroutine overwrites any previous computation stored in the object.
        !! @endwarning
        !!
        !! @note
        !! **Convergence criterion**  
        !! See [[is_not_converged]]
        !! @endnote

        class(arithmetic_geometric_mean_real64_type), intent(inout) :: agm

        real(real64), intent(in) :: init_a
        !! initial value: arithmetic mean

        real(real64), intent(in) :: init_g
        !! initial value: geometric mean



        call initialize(agm)

        agm%list_a(0) = init_a
        agm%list_g(0) = init_g
        agm%n_iter_    = 0

        do

            associate(prev_iter => agm%n_iter_, next_iter => agm%n_iter_ + 1)

                associate( &!
                    prev_a => agm%list_a(prev_iter) , &!
                    prev_g => agm%list_g(prev_iter) , &!
                    next_a => agm%list_a(next_iter) , &!
                    next_g => agm%list_g(next_iter)   &!
                )

                    agm%n_iter_ = agm%n_iter_ + 1

                    call compute_step( &!
                        prev_a = prev_a , &!
                        prev_g = prev_g , &!
                        next_a = next_a , &!
                        next_g = next_g   &!
                    )

                    if ( is_not_converged(next_a, next_g) ) then

                        prev_iter = next_iter

                        cycle

                    else

                        prev_iter = next_iter

                        return

                    end if

                end associate

            end associate

        end do

    end subroutine compute_kernel_real64



    elemental subroutine compute_kernel_real128(agm, init_a, init_g)
        !! Compute arithmetic–geometric mean (AGM)
        !! using the given arithmetic mean `init_a` and geometric mean `init_g`
        !! and store the complete iteration history.
        !!
        !! This type-bound subroutine performs the AGM iteration starting from
        !! given initial values, storing each intermediate arithmetic and geometric mean
        !! in `list_a` and `list_g` respectively.
        !! The final iteration count is stored in `n_iter`.
        !!
        !! Unlike the [[arithmetic_geometric_mean_kernel]],
        !! this subroutine preserves the full convergence history for analysis.
        !!
        !! @warning
        !! - This subroutine/interface assumes both inputs are positive.
        !! - No validation is performed on inputs.
        !! - Calling this subroutine overwrites any previous computation stored in the object.
        !! @endwarning
        !!
        !! @note
        !! **Convergence criterion**  
        !! See [[is_not_converged]]
        !! @endnote

        class(arithmetic_geometric_mean_real128_type), intent(inout) :: agm

        real(real128), intent(in) :: init_a
        !! initial value: arithmetic mean

        real(real128), intent(in) :: init_g
        !! initial value: geometric mean



        call initialize(agm)

        agm%list_a(0) = init_a
        agm%list_g(0) = init_g
        agm%n_iter_    = 0

        do

            associate(prev_iter => agm%n_iter_, next_iter => agm%n_iter_ + 1)

                associate( &!
                    prev_a => agm%list_a(prev_iter) , &!
                    prev_g => agm%list_g(prev_iter) , &!
                    next_a => agm%list_a(next_iter) , &!
                    next_g => agm%list_g(next_iter)   &!
                )

                    agm%n_iter_ = agm%n_iter_ + 1

                    call compute_step( &!
                        prev_a = prev_a , &!
                        prev_g = prev_g , &!
                        next_a = next_a , &!
                        next_g = next_g   &!
                    )

                    if ( is_not_converged(next_a, next_g) ) then

                        prev_iter = next_iter

                        cycle

                    else

                        prev_iter = next_iter

                        return

                    end if

                end associate

            end associate

        end do

    end subroutine compute_kernel_real128



    elemental subroutine compute_step_real32(prev_a, prev_g, next_a, next_g)
        !! Compute arithmetic and geometric mean using the given `prev_a` and `prev_g`.  
        !!
        !! @warning
        !! - This subroutine assumes both inputs are positive.
        !! - No validation is performed on inputs.
        !! @endwarning


        real(real32), intent(in) :: prev_a
        !! previous arithmetic mean

        real(real32), intent(in) :: prev_g
        !! previous geometric mean

        real(real32), intent(out) :: next_a
        !! next arithmetic mean

        real(real32), intent(out) :: next_g
        !! next geometric mean



        next_a =     (prev_a + prev_g) * 0.5_real32
        next_g = sqrt(prev_a * prev_g)

    end subroutine compute_step_real32



    elemental subroutine compute_step_real64(prev_a, prev_g, next_a, next_g)
        !! Compute arithmetic and geometric mean using the given `prev_a` and `prev_g`.  
        !!
        !! @warning
        !! - This subroutine assumes both inputs are positive.
        !! - No validation is performed on inputs.
        !! @endwarning


        real(real64), intent(in) :: prev_a
        !! previous arithmetic mean

        real(real64), intent(in) :: prev_g
        !! previous geometric mean

        real(real64), intent(out) :: next_a
        !! next arithmetic mean

        real(real64), intent(out) :: next_g
        !! next geometric mean



        next_a =     (prev_a + prev_g) * 0.5_real64
        next_g = sqrt(prev_a * prev_g)

    end subroutine compute_step_real64



    elemental subroutine compute_step_real128(prev_a, prev_g, next_a, next_g)
        !! Compute arithmetic and geometric mean using the given `prev_a` and `prev_g`.  
        !!
        !! @warning
        !! - This subroutine assumes both inputs are positive.
        !! - No validation is performed on inputs.
        !! @endwarning


        real(real128), intent(in) :: prev_a
        !! previous arithmetic mean

        real(real128), intent(in) :: prev_g
        !! previous geometric mean

        real(real128), intent(out) :: next_a
        !! next arithmetic mean

        real(real128), intent(out) :: next_g
        !! next geometric mean



        next_a =     (prev_a + prev_g) * 0.5_real128
        next_g = sqrt(prev_a * prev_g)

    end subroutine compute_step_real128



    elemental subroutine initialize_real32(agm)
        !! Initialize components: `n_iter`, `list_a` and `list_g`

        type(arithmetic_geometric_mean_real32_type), intent(inout) :: agm



        agm%n_iter_ = initial_n_iter

        agm%list_a(:) = ieee_value( x = 0.0_real32, class = ieee_quiet_nan )
        agm%list_g(:) = agm%list_a(:)

    end subroutine initialize_real32



    elemental subroutine initialize_real64(agm)
        !! Initialize components: `n_iter`, `list_a` and `list_g`

        type(arithmetic_geometric_mean_real64_type), intent(inout) :: agm



        agm%n_iter_ = initial_n_iter

        agm%list_a(:) = ieee_value( x = 0.0_real64, class = ieee_quiet_nan )
        agm%list_g(:) = agm%list_a(:)

    end subroutine initialize_real64



    elemental subroutine initialize_real128(agm)
        !! Initialize components: `n_iter`, `list_a` and `list_g`

        type(arithmetic_geometric_mean_real128_type), intent(inout) :: agm



        agm%n_iter_ = initial_n_iter

        agm%list_a(:) = ieee_value( x = 0.0_real128, class = ieee_quiet_nan )
        agm%list_g(:) = agm%list_a(:)

    end subroutine initialize_real128

end module arithmetic_geometric_mean_fortran
