cnd.h

#pragma once

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

#include <cmath>

namespace funcy
{
    struct CumulativeNormalDistribution : Chainer< CumulativeNormalDistribution >
    {
        explicit CumulativeNormalDistribution( double x = 0. ) noexcept
        {
            update( x );
        }

        void update( double x ) noexcept
        {
            x_ = x;
            value_ = 0.5 * ( 1 + std::erf( x_ * div_sqrt2 ) );
            d1_ = std::exp( -0.5 * x_ * x_ );
        }

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

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

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

        [[nodiscard]] double d3( double dx = 1., double dy = 1., double dz = 1. ) const noexcept
        {
            return div_sqrt2pi * ( -d1_ + x_ * x_ * d1_ ) * dx * dy * dz;
        }

    private:
        static const double div_sqrt2;
        static const double div_sqrt2pi;
        double value_;
        double d1_;
        double x_;
    };

    inline const double CumulativeNormalDistribution::div_sqrt2 = 1 / std::sqrt( 2 );
    inline const double CumulativeNormalDistribution::div_sqrt2pi = 1 / std::sqrt( 2 * M_PI );

    [[nodiscard]] inline auto cnd( double x )
    {
        return CumulativeNormalDistribution( x )();
    }

    template < Function F >
    [[nodiscard]] auto cnd( const F& f )
    {
        return CumulativeNormalDistribution()( f );
    }
} // namespace funcy