pow.h

#pragma once

#include <funcy/concepts.h>
#include <funcy/util/chainer.h>
#include <funcy/util/exceptions.h>
#include <funcy/util/macros.h>

#include <cmath>

namespace funcy
{
    template < int dividend, int divisor = 1 >
    struct Pow : Chainer< Pow< dividend, divisor > >
    {
        explicit Pow( double x = 1 ) FUNCY_NOEXCEPT
        {
            update( x );
        }

        void update( double x ) FUNCY_NOEXCEPT
        {
#ifdef FUNCY_ENABLE_EXCEPTIONS
            if ( k < 3 && x == 0 )
                throw OutOfDomainException( "Pow<" + std::to_string( dividend ) + "," +
                                                std::to_string( divisor ) + ">",
                                            "]-inf,inf[ \\ {0}", x, __FILE__, __LINE__ );
#endif
            xk = x * ( xk1 = x * ( xk2 = x * ( xk3 = ::pow( x, k - 3 ) ) ) );
        }

        [[nodiscard]] double d0() const noexcept
        {
            return xk;
        }

        [[nodiscard]] double d1( double dx = 1. ) const noexcept
        {
            return k * xk1 * dx;
        }

        [[nodiscard]] double d2( double dx = 1., double dy = 1. ) const noexcept
        {
            return k * ( k - 1 ) * xk2 * dx * dy;
        }

        [[nodiscard]] double d3( double dx = 1., double dy = 1., double dz = 1. ) const noexcept
        {
            return k * ( k - 1 ) * ( k - 2 ) * xk3 * dx * dy * dz;
        }

    private:
        const double k = static_cast< double >( dividend ) / divisor;
        double xk = 0;
        double xk1 = 0;
        double xk2 = 0;
        double xk3 = 0;
    };

    template <>
    struct Pow< 2, 1 > : Chainer< Pow< 2, 1 > >
    {
        explicit Pow( double x_ = 0 ) FUNCY_NOEXCEPT
        {
            update( x_ );
        }

        void update( const double& x_ ) FUNCY_NOEXCEPT
        {
            x = 2 * x_;
            x2 = x_ * x_;
        }

        [[nodiscard]] double d0() const noexcept
        {
            return x2;
        }

        [[nodiscard]] double d1( double dx = 1. ) const noexcept
        {
            return x * dx;
        }

        [[nodiscard]] double d2( double dx = 1, double dy = 1 ) const noexcept
        {
            return 2 * dx * dy;
        }

    private:
        double x = 0., x2 = 0.;
    };

    template <>
    struct Pow< 3, 1 > : Chainer< Pow< 3, 1 > >
    {
        explicit Pow( double x_ = 0 ) FUNCY_NOEXCEPT
        {
            update( x_ );
        }

        void update( double x_ ) FUNCY_NOEXCEPT
        {
            x = x_;
            x2 = x * x;
            x3 = x2 * x;
        }

        [[nodiscard]] double d0() const noexcept
        {
            return x3;
        }

        [[nodiscard]] double d1( double dx = 1 ) const noexcept
        {
            return 3 * x2 * dx;
        }

        [[nodiscard]] double d2( double dx = 1, double dy = 1 ) const noexcept
        {
            return 6 * x * dx * dy;
        }

        [[nodiscard]] double d3( double dx = 1, double dy = 1, double dz = 1 ) const noexcept
        {
            return 6 * dx * dy * dz;
        }

    private:
        double x = 0;
        double x2 = 0;
        double x3 = 0.;
    };

    template <>
    struct Pow< -1, 1 > : Chainer< Pow< -1, 1 > >
    {
        explicit Pow( double x = 1. ) FUNCY_NOEXCEPT
        {
            update( x );
        }

        void update( double x ) FUNCY_NOEXCEPT
        {
#ifdef FUNCY_ENABLE_EXCEPTIONS
            if ( x == 0 )
                throw OutOfDomainException( "Pow<-1,1>", "]-inf,inf[ \\ {0}", x, __FILE__,
                                            __LINE__ );
#endif
            x_inv = 1. / x;
            x_inv2 = x_inv * x_inv;
        }

        [[nodiscard]] double d0() const noexcept
        {
            return x_inv;
        }

        [[nodiscard]] double d1( double dx = 1. ) const noexcept
        {
            return -1 * x_inv2 * dx;
        }

        [[nodiscard]] double d2( double dx = 1., double dy = 1. ) const noexcept
        {
            return 2 * x_inv2 * x_inv * dx * dy;
        }

        [[nodiscard]] double d3( double dx = 1., double dy = 1., double dz = 1. ) const noexcept
        {
            return -6 * x_inv2 * x_inv2 * dx * dy * dz;
        }

    private:
        double x_inv = 1;
        double x_inv2 = 1;
    };

    template <>
    struct Pow< 1, 2 > : Chainer< Pow< 1, 2 > >
    {
        explicit Pow( double x = 0 ) FUNCY_NOEXCEPT
        {
            update( x );
        }

        void update( double x )
        {
#ifdef FUNCY_ENABLE_EXCEPTIONS
            if ( x < 0 )
                throw OutOfDomainException( "Pow<1,2>", "[0,inf[", x, __FILE__, __LINE__ );
#endif
            x_ = x;
            sqrt_x = ::sqrt( x );
        }

        [[nodiscard]] double d0() const noexcept
        {
            return sqrt_x;
        }

        [[nodiscard]] double d1( double dx = 1. ) const noexcept
        {
            return 0.5 / sqrt_x * dx;
        }

        [[nodiscard]] double d2( double dx = 1., double dy = 1. ) const noexcept
        {
            return -0.25 / ( x_ * sqrt_x ) * dx * dy;
        }

        [[nodiscard]] double d3( double dx = 1., double dy = 1., double dz = 1. ) const noexcept
        {
            return 0.375 / ( x_ * x_ * sqrt_x ) * dx * dy * dz;
        }

    private:
        double x_ = 0;
        double sqrt_x = 1;
    };

    template <>
    struct Pow< -1, 3 > : Chainer< Pow< -1, 3 > >
    {
        explicit Pow( double t = 1 ) FUNCY_NOEXCEPT
        {
            update( t );
        }

        void update( double x ) FUNCY_NOEXCEPT
        {
#ifdef FUNCY_ENABLE_EXCEPTIONS
            if ( x < 0 )
                throw OutOfDomainException( "Pow<1,3>", "[0,inf[", x, __FILE__, __LINE__ );
#endif
            auto p = cbrt( x );
            d0val = 1 / p;
            p *= x;
            d1val = -1 / ( 3 * p );
            p *= x;
            d2val = 4 / ( 9 * p );
            p *= x;
            d3val = -28 / ( 27 * p );
        }

        [[nodiscard]] double d0() const noexcept
        {
            return d0val;
        }

        [[nodiscard]] double d1( double dt = 1 ) const noexcept
        {
            return d1val * dt;
        }

        [[nodiscard]] double d2( double dt0 = 1, double dt1 = 1 ) const noexcept
        {
            return d2val * dt0 * dt1;
        }

        [[nodiscard]] double d3( double dt0 = 1, double dt1 = 1, double dt2 = 1 ) const noexcept
        {
            return d3val * dt0 * dt1 * dt2;
        }

    private:
        double d0val = 0;
        double d1val = 0;
        double d2val = 0;
        double d3val = 0;
    };

    template <>
    struct Pow< -2, 3 > : Chainer< Pow< -2, 3 > >
    {
        explicit Pow( double t = 1. ) FUNCY_NOEXCEPT
        {
            update( t );
        }

        void update( double x ) FUNCY_NOEXCEPT
        {
#ifdef FUNCY_ENABLE_EXCEPTIONS
            if ( x < 0 )
                throw OutOfDomainException( "Pow<2,3>", "[0,inf[", x, __FILE__, __LINE__ );
#endif
            auto p0 = cbrt( x );
            auto p = p0 * p0;
            d0val = 1 / p;
            p *= x;
            d1val = -2 / ( 3 * p );
            p *= x;
            d2val = 10 / ( 9 * p );
            p *= x;
            d3val = -80 / ( 27 * p );
        }

        [[nodiscard]] double d0() const noexcept
        {
            return d0val;
        }

        [[nodiscard]] double d1( double dt = 1 ) const noexcept
        {
            return d1val * dt;
        }

        [[nodiscard]] double d2( double dt0 = 1, double dt1 = 1 ) const noexcept
        {
            return d2val * dt0 * dt1;
        }

        [[nodiscard]] double d3( double dt0 = 1, double dt1 = 1, double dt2 = 1 ) const noexcept
        {
            return d3val * dt0 * dt1 * dt2;
        }

    private:
        double d0val = 0;
        double d1val = 0;
        double d2val = 0;
        double d3val = 0;
    };

    using Sqrt = Pow< 1, 2 >;

    using Cbrt = Pow< 1, 3 >;

    using Cbrt2 = Pow< 2, 3 >;

    template < Function F >
    [[nodiscard]] auto sqrt( const F& f )
    {
        return Sqrt()( f );
    }

    template < Function F >
    [[nodiscard]] auto cbrt( const F& f )
    {
        return Cbrt()( f );
    }

    template < Function F >
    [[nodiscard]] auto cbrt2( const F& f )
    {
        return Cbrt2()( f );
    }

    template < int k, int l, Function F >
    [[nodiscard]] auto pow( const F& f )
    {
        return Pow< k, l >()( f );
    }

    template < int k, Function F >
    [[nodiscard]] auto pow( const F& f )
    {
        return Pow< k >()( f );
    }
} // namespace funcy