finite_charge_sum Function

  pure real(dp) function finite_charge_sum(values, context) result(total)
    real(dp), intent(in) :: values(:)
    character(len=*), intent(in) :: context
    real(dp) :: compensation, correction, scaled_sum, scaled_value, scale_value, updated_sum
    integer :: i

    if (.not. all(ieee_is_finite(values))) error stop trim(context)//' contains a non-finite charge.'
    if (size(values) == 0) then
      total = 0.0_dp
      return
    end if
    scale_value = maxval(abs(values))
    if (scale_value == 0.0_dp) then
      total = 0.0_dp
      return
    end if
    scaled_sum = 0.0_dp
    compensation = 0.0_dp
    do i = 1, size(values)
      scaled_value = values(i)/scale_value
      updated_sum = scaled_sum + scaled_value
      ! 局所補正を先に確定し、既存の compensation が式中間で丸め落ちるのを防ぐ。
      if (abs(scaled_sum) >= abs(scaled_value)) then
        correction = (scaled_sum - updated_sum) + scaled_value
      else
        correction = (scaled_value - updated_sum) + scaled_sum
      end if
      compensation = compensation + correction
      scaled_sum = updated_sum
    end do
    scaled_sum = scaled_sum + compensation
    if (scale_value >= 1.0_dp) then
      if (abs(scaled_sum) > huge(total)/scale_value) error stop trim(context)//' charge sum overflowed.'
    end if
    total = scale_value*scaled_sum
    if (.not. ieee_is_finite(total)) error stop trim(context)//' charge sum is not finite.'
  end function finite_charge_sum