/*=========================================================================
 *
 *  Copyright NumFOCUS
 *
 *  Licensed under the Apache License, Version 2.0 (the "License");
 *  you may not use this file except in compliance with the License.
 *  You may obtain a copy of the License at
 *
 *         https://www.apache.org/licenses/LICENSE-2.0.txt
 *
 *  Unless required by applicable law or agreed to in writing, software
 *  distributed under the License is distributed on an "AS IS" BASIS,
 *  WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 *  See the License for the specific language governing permissions and
 *  limitations under the License.
 *
 *=========================================================================*/
#ifndef itkCompensatedSummation_hxx
#define itkCompensatedSummation_hxx


namespace itk
{

void ITKCommon_EXPORT
     CompensatedSummationAddElement(float & compensation, float & sum, const float element);
void ITKCommon_EXPORT
     CompensatedSummationAddElement(double & compensation, double & sum, const double element);

#ifndef itkCompensatedSummation_cxx
// We try the looser pragma guards if we don't have an explicit instantiation.
#  ifdef __INTEL_COMPILER
#    pragma optimize("", off)
#  endif // __INTEL_COMPILER
#  ifdef _MSC_VER
#    pragma float_control(push)
#    pragma float_control(precise, on)
#  endif // _MSC_VER i.e. Microsoft Visual Studio
#endif   // not itkCompensatedSummation_cxx
template <typename TFloat>
/** A helper for the CompensatedSummation class. */
void
CompensatedSummationAddElement(TFloat & compensation, TFloat & sum, const TFloat & element, int = 0)
{
  using AccumulateType = typename NumericTraits<TFloat>::AccumulateType;
  const auto           compensatedInput = static_cast<AccumulateType>(element - compensation);
  const AccumulateType tempSum = sum + compensatedInput;
  // Warning: watch out for the compiler optimizing this out!
  compensation = static_cast<TFloat>((tempSum - sum) - compensatedInput);
  sum = static_cast<TFloat>(tempSum);
}
#ifndef itkCompensatedSummation_cxx
#  ifdef __INTEL_COMPILER
#    pragma optimize("", on)
#  endif // __INTEL_COMPILER
#  ifdef _MSC_VER
#    pragma float_control(pop)
#  endif // _MSC_VER
#endif   // not itkCompensatedSummation_cxx

template <typename TFloat>
CompensatedSummation<TFloat>::CompensatedSummation(const TFloat value)
  : m_Sum(value)
  , m_Compensation(AccumulateType{})
{}

template <typename TFloat>
CompensatedSummation<TFloat>::CompensatedSummation(const Self & rhs)
{
  this->m_Sum = rhs.m_Sum;
  this->m_Compensation = rhs.m_Compensation;
}

template <typename TFloat>
auto
CompensatedSummation<TFloat>::operator=(const Self & rhs) -> Self &
{
  if (this != &rhs)
  {
    this->m_Sum = rhs.m_Sum;
    this->m_Compensation = rhs.m_Compensation;
  }
  return *this;
}

template <typename TFloat>
void
CompensatedSummation<TFloat>::AddElement(const FloatType & element)
{
  CompensatedSummationAddElement(this->m_Compensation, this->m_Sum, element);
}

template <typename TFloat>
CompensatedSummation<TFloat> &
CompensatedSummation<TFloat>::operator+=(const FloatType & rhs)
{
  this->AddElement(rhs);
  return *this;
}


template <typename TFloat>
CompensatedSummation<TFloat> &
CompensatedSummation<TFloat>::operator+=(const Self & rhs)
{
  this->AddElement(rhs.m_Compensation);
  this->AddElement(rhs.m_Sum);
  return *this;
}

template <typename TFloat>
CompensatedSummation<TFloat> &
CompensatedSummation<TFloat>::operator-=(const FloatType & rhs)
{
  this->AddElement(-rhs);
  return *this;
}

template <typename TFloat>
CompensatedSummation<TFloat> &
CompensatedSummation<TFloat>::operator*=(const FloatType & rhs)
{
  this->m_Sum *= rhs;
  this->m_Compensation *= rhs;
  return *this;
}

template <typename TFloat>
CompensatedSummation<TFloat> &
CompensatedSummation<TFloat>::operator/=(const FloatType & rhs)
{
  this->m_Sum /= rhs;
  this->m_Compensation /= rhs;
  return *this;
}

template <typename TFloat>
void
CompensatedSummation<TFloat>::ResetToZero()
{
  this->m_Sum = AccumulateType{};
  this->m_Compensation = AccumulateType{};
}

template <typename TFloat>
CompensatedSummation<TFloat> &
CompensatedSummation<TFloat>::operator=(const FloatType & rhs)
{
  this->m_Sum = rhs;
  this->m_Compensation = AccumulateType{};

  return *this;
}

template <typename TFloat>
auto
CompensatedSummation<TFloat>::GetSum() const -> const AccumulateType &
{
  return this->m_Sum;
}


template <typename TFloat>
CompensatedSummation<TFloat>::operator typename CompensatedSummation<TFloat>::FloatType() const
{
  return this->m_Sum;
}

} // end namespace itk

#endif