compute_real32 Subroutine

private 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

Warning

Calling this subroutine overwrites any previous computation stored in the object.

Type Bound

arithmetic_geometric_mean_real32_type

Arguments

Type IntentOptional Attributes Name
class(arithmetic_geometric_mean_real32_type), intent(inout) :: agm
real(kind=real32), intent(in) :: x
real(kind=real32), intent(in) :: y

Variables

Type Visibility Attributes Name Initial
real(kind=real32), private :: xy

Source Code

    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