#pragma once

#ifndef _GEAR_TESTING__MATH__MULTIPRECISION_H_
#define _GEAR_TESTING__MATH__MULTIPRECISION_H_

#include <gear_testing/math/floatingpoint.h>

#if defined(min)
#   undef min
#endif

#if defined(max)
#   undef max
#endif

namespace G4
{

    class MPNumeric
    {
    public:

        MPNumeric(double value = 0.0)
        : m_value(value)
        {
            // Note that we somewhat arbitrarily choose to give 0 the same error as 1
            m_error = AsFloatingPoint(value!=0.0?value:1.0).ULPToAbsolute(2);
        }

        MPNumeric(float value)
        : m_value(value)
        {
            // Note that we somewhat arbitrarily choose to give 0 the same error as 1
            m_error = AsFloatingPoint(value!=0.0f?value:1.0f).ULPToAbsolute(2);
        }

        MPNumeric operator-() const
        {
            return MPNumeric(-m_value, m_error);
        }

        MPNumeric& operator+=(const MPNumeric& value)
        {
            m_value += value.m_value;
            m_error += value.m_error;
            return *this;
        }

        MPNumeric operator+(const MPNumeric& value) const
        {
            MPNumeric t(*this);
            t += value;
            return t;
        }

        MPNumeric& operator-=(const MPNumeric& value)
        {
            m_value -= value.m_value;
            m_error += value.m_error;
            return *this;
        }

        MPNumeric operator-(const MPNumeric& value) const
        {
            MPNumeric t(*this);
            t -= value;
            return t;
        }

        MPNumeric& operator*=(const MPNumeric& value)
        {
            m_error = AddRelativeError(value) * 1.2;
            m_value *= value.m_value;
            m_error *= std::abs(m_value);
            return *this;
        }

        MPNumeric operator*(const MPNumeric& value) const
        {
            MPNumeric t(*this);
            t *= value;
            return t;
        }

        MPNumeric& operator/=(const MPNumeric& value)
        {
            if ( AsFloatingPoint( value.m_value ).IsInfinite() )
            {
                m_value = 0.0;
                // Note that we somewhat arbitrarily choose to give 0 the same error as 1
                m_error = AsFloatingPoint( 1.0 ).ULPToAbsolute(2);
                return( *this );
            }

            // This is basically doing the same thing as interval arithmetic in an attempt to find good bounds
            double lower[2];
            lower[0] = m_value - m_error;
            lower[1] = value.m_value - value.m_error;
            double upper[2];
            upper[0] = m_value + m_error;
            upper[1] = value.m_value + value.m_error;
            double interval[2];
            double val = lower[0] / lower[ 1 ];
            interval[ 0 ] = val;
            interval[ 1 ] = val;
            val = lower[0] / upper[ 1 ];
            interval[ 0 ] = std::min( interval[ 0 ], val );
            interval[ 1 ] = std::max( interval[ 1 ], val );
            val = upper[0] / lower[ 1 ];
            interval[ 0 ] = std::min( interval[ 0 ], val );
            interval[ 1 ] = std::max( interval[ 1 ], val );
            val = upper[0] / upper[ 1 ] ;
            interval[ 0 ] = std::min( interval[ 0 ], val );
            interval[ 1 ] = std::max( interval[ 1 ], val );

            m_error = AddRelativeError(value);
            m_value /= value.m_value;
            m_error *= std::abs(m_value);

            double intervalError = std::max( std::abs( m_value - interval[0] ), std::abs( m_value - interval[1] ) );
            m_error = std::max( m_error, intervalError ); // Use nominal error value if the error determined by interval arithmetic is not tight enough
            return *this;
        }

        MPNumeric operator/(const MPNumeric& value) const
        {
            MPNumeric t(*this);
            t /= value;
            return t;
        }

        MPNumeric Abs() const
        {
            return( MPNumeric( std::abs( m_value ), m_error ) );
        }

        MPNumeric Cosine() const
        {
            double value = std::cos(m_value);
            double err1 = std::abs( value - std::cos( m_value + m_error ) );
            double err2 = std::abs( value - std::cos( m_value - m_error ) );
            double error = std::max( err1, err2 ) * 2.0;
            double err3 = AsFloatingPoint(1.0f).ULPToAbsolute(2);
            error = std::max( error, err3 ) * 2.0; // Use nominal error value if the error determined by interval arithmetic is too tight
            return( MPNumeric( value, error ) );
        }

        MPNumeric Sine() const
        {
            double value = std::sin(m_value);
            double err1 = std::abs( value - std::sin( m_value + m_error ) );
            double err2 = std::abs( value - std::sin( m_value - m_error ) );
            double error = std::max( err1, err2 ) * 2.0;
            double err3 = AsFloatingPoint(1.0f).ULPToAbsolute(2);
            error = std::max( error, err3 ) * 2.0; // Use nominal error value if the error determined by interval arithmetic is too tight
            return( MPNumeric( value, error ) );
        }

        MPNumeric ArcCos() const
        {
            double value = std::acos(m_value);
            double err1 = std::abs( value - std::acos( std::min( m_value + m_error, 1.0 ) ) );
            double err2 = std::abs( value - std::acos( std::max( m_value - m_error, -1.0 ) ) );
            double error = std::max( err1, err2 ) * 2.0;
            double err3 = AsFloatingPoint(G4PI).ULPToAbsolute(2);
            error = std::max( error, err3 ); // Use nominal error value if the error determined by interval arithmetic is too tight
            return( MPNumeric( value, error ) );
        }

        MPNumeric ArcTan()
        {
            double value = std::atan( m_value );
            // Determine error bounds in an interval arithmetic like manner
            double error = std::abs( value - std::atan( m_value + m_error ) );
            double err2 = std::abs( value - std::atan( m_value - m_error ) );
            error = std::max ( error, err2 );
            return( MPNumeric( value, error ) );
        }

        template < typename _SCALAR_TYPE >
        static MPNumeric ArcTan2( const MPNumeric& numer, const MPNumeric& denom )
        {
            // The multiprecision ArcTan2 does not use standard library atan2.  It is too difficult to determine a reasonable error bounds using only that function.
            // Instead, we use the implementation below which is (basically) the same as the internals of some standard library implementations.
            // This allows us to more accurately accumulate the error bounds.
            MPNumeric temp = ( ( numer * numer + denom * denom ).Sqrt() - denom ) / numer;
            MPNumeric result = MPNumeric( _SCALAR_TYPE( 2 ) ) * temp.ArcTan();
            return( result );
        }

        MPNumeric Sqrt() const
        {
            // rsqrte usually return a 12 bits precision in the mantissa.
            // This correspond to a relative error of 1/4096
            // Note: we usually get a better precision however, because of
            // Newton-Raphson iteration add to it.
            double error = 2*GetRelativeError();
            static const double sqrterr = 1./4096.;
            if (error < sqrterr)
                error = sqrterr;

            double value = std::sqrt(m_value);

            error *= std::abs( value );

            if ( !IsTypeValid( value ) )
            {
                return( MPNumeric( value, error ) );
            }

            // This is basically doing the same thing as interval arithmetic in an attempt to find better bounds
            double err1 = std::abs( value - std::sqrt( m_value + m_error ) );
            double err2 = std::abs( value - std::sqrt( std::max( m_value - m_error, 0.0 ) ) );
            double intervalError = std::max( err1, err2 ) * 2.0;
            error = std::max( error, intervalError );

            return( MPNumeric( value, error ) );
        }

        template <class T>
        T GetValue() const
        {
            return static_cast<T>(m_value);
        }

        double GetRelativeError() const
        {
            if ( !IsTypeValid( m_value ) )
            {
                return( m_value );
            }
            // Not sure if we're doing the right thing with zero values here
            return ( m_value != 0.0 ) ? m_error/std::abs(m_value) : m_error;
        }

        double GetAbsoluteError() const
        {
            return m_error;
        }

        template <class T>
        T GetError() const
        {
            FloatingPoint<T> v(std::max<T>(std::abs(static_cast<T>(m_value)), 1.0f));
            return std::max(static_cast<T>(m_error), v.ULPToAbsolute(1));
        }

    protected:
        MPNumeric(double value, double error)
            : m_value(value)
            , m_error(error)
        {}

        double AddRelativeError(const MPNumeric& value) const
        {
            return GetRelativeError()+value.GetRelativeError();
        }

    private:

        double m_value;
        double m_error;
    };
} // namspace G4

// Stream operations are not supported on SPU
namespace G4
{
#if !defined(__SPU__)
    inline std::ostream& operator<<(std::ostream& os, const MPNumeric& value)
    {
        return os << "{ value: " << AsFloatingPoint(value.GetValue<double>())
            << ", error: " << value.GetAbsoluteError() << " }";
    }
#endif

    struct MPNumericPredFormat
    {
        template <class T>
        ::testing::AssertionResult operator()(const char* expExpr,
                                   const char* resExpr,
                                   G4::MPNumeric expected,
                                   T result)
        {
            T diff = std::abs(expected.GetValue<T>()- result);
            T err = expected.GetError<T>();
            if (diff <= err)
                return ::testing::AssertionSuccess();
            return ::testing::AssertionFailure()
                << resExpr << " (" << AsFloatingPoint(result)
                << ") differed from " << expExpr << " (" << AsFloatingPoint(expected.GetValue<T>())
                << ") by " << diff
                << " which is greater than " << AsFloatingPoint(expected.GetError<T>());
        }
    };
} // namspace G4

#if !defined(__SPU__)

# define EXPECT_MP_EQ(TYPE, MP1, V2) ASSERT_PRED_FORMAT2(::G4::MPNumericPredFormat(), MP1, V2)

#else

# define EXPECT_MP_EQ(TYPE, MP1, V2)    __EXPECT_CALL( ::G4::MPNumericPredFormat(), "", "", (MP1), (V2) )

# define ASSERT_MP_EQ(TYPE, MP1, V2)    __ASSERT_CALL( EXPECT_MP_EQ,    (MP1), (V2))

#endif

#endif // #ifndef _GEAR_TESTING__MATH__MULTIPRECISION_H_