Ziggurat法の定数計算用関数

main.jl

# The Ziggurat Method for generating random variables - Marsaglia and Tsang
# Paper and reference code: http://www.jstatsoft.org/v05/i08/

module ZigguratTable

    function compute_core( table_size::Int, pdf0::T, r::T, pdf::Core.Function, integral_pdf::Core.Function, inverse_pdf::Core.Function ) where { T <: Core.AbstractFloat }

        v     =    ( r * pdf( r ) ) + integral_pdf( r )
        x_ref =    r
        x_new =    Base.zero( T )
        z_err =  - Base.one( T )

        for _ = ( table_size - 2 ) : -1 : 1

            try
                x_new = inverse_pdf( pdf( x_ref ) + ( v / x_ref ) )
            catch
                return v, z_err
            finally

                if x_new < Base.zero( T )
                    return v, z_err
                else
                    x_ref = x_new
                end
            end

        end

        return v, ( ( x_new * ( pdf0 - pdf( x_new ) ) ) - v )

    end



    function compute( n_cycle::Int, table_size::Int, r_min_ini::T, r_max_ini::T, pdf::Core.Function, integral_pdf::Core.Function, inverse_pdf::Core.Function, flag_display::Core.Bool ) where { T <: Core.AbstractFloat }

        r_min =   r_min_ini
        r_mid = ( r_min_ini + r_max_ini ) / 2
        r_max =               r_max_ini

        pdf0  = pdf( Base.zero( T ) )

        v_min , z_min = compute_core( table_size, pdf0, r_min, pdf, integral_pdf, inverse_pdf )
        v_mid , z_mid = compute_core( table_size, pdf0, r_mid, pdf, integral_pdf, inverse_pdf )
        v_max , z_max = compute_core( table_size, pdf0, r_max, pdf, integral_pdf, inverse_pdf )

        for _ in Base.OneTo( n_cycle )

            if ( z_min * z_mid ) > Base.zero( T )

                r_min = r_mid
                z_min = z_mid

            elseif ( z_max * z_mid ) > Base.zero( T )

                r_max = r_mid
                z_max = z_mid

            end

            r_mid        = ( r_min + r_max ) / 2
            v_mid, z_mid = compute_core( table_size, pdf0, r_mid, pdf, integral_pdf, inverse_pdf )

            if flag_display
                println( r_min, " ", r_mid, " ", r_max, " ", z_min, " ", z_mid, " ", z_max )
            end

            if isequal( r_min, r_mid ) || isequal( r_mid, r_max )
                return r_mid, v_mid, z_mid, true
            end

        end

        return r_mid, v_mid, z_mid, false

    end



    function test_r_vz( filename::Core.AbstractString, r_step::Core.AbstractFloat, r_oneto::Int, pdf::Core.Function, integral_pdf::Core.Function, inverse_pdf::Core.Function )

        open( filename, "w" ) do file

            pdf0 = pdf( Base.zero( r_step ) )

            for i in Base.oneto( r_oneto )

                r = r_step * i

                v, z =
                    ZigguratTable.compute_core(
                        128,
                        pdf0,
                        r,
                        ZigguratTable.Normal.pdf,
                        ZigguratTable.Normal.integral_pdf,
                        ZigguratTable.Normal.inverse_pdf
                    )
        
                Base.write( file, r, v, z )

            end

        end

    end



    module Exponential

        using ..ZigguratTable



        function pdf( x::Core.AbstractFloat )
            return Base.exp( -x )
        end



        function integral_pdf( x::Core.AbstractFloat )
            return pdf(x)
        end



        function inverse_pdf( x::Core.AbstractFloat )
            return ( -Base.log( x ) )
        end



        function specialized_compute( table_size::Int, T::Core.Type, flag_display::Core.Bool )

            r_min_ini = Base.one( T )
            r_max_ini = Base.one( T ) * 10

            return (
                    ZigguratTable.compute(
                        512,
                        table_size,
                        r_min_ini,
                        r_max_ini,
                        pdf,
                        integral_pdf,
                        inverse_pdf,
                        flag_display
                    )
                )

        end

    end



    module Normal

        using ..ZigguratTable
        using SpecialFunctions



        function pdf( x::Core.AbstractFloat )
            return Base.exp( - ( x * x ) / 2 )
        end



        function integral_pdf( x::Core.AbstractFloat )
            return Base.sqrt( Base.MathConstants.pi / 2 ) * SpecialFunctions.erfc( x / Base.sqrt( 2 * Base.one(x) ) )
        end



        function inverse_pdf( x::Core.AbstractFloat )
            return Base.sqrt( -2 * Base.log( x ) )
        end



        function specialized_compute( table_size::Int, T::Core.Type, flag_display::Core.Bool )

            r_min_ini = Base.one( T )
            r_max_ini = r_min_ini * 10

            return (
                    ZigguratTable.compute(
                        512,
                        table_size,
                        r_min_ini,
                        r_max_ini,
                        pdf,
                        integral_pdf,
                        inverse_pdf,
                        flag_display
                    )
                )

        end
    end

end



using Random

# ZigguratTable.test_r_vz(
#     "log.dat",
#     ( Base.one( Core.Float64 ) / 100 ),
#     1000,
#     ZigguratTable.Normal.pdf,
#     ZigguratTable.Normal.integral_pdf,
#     ZigguratTable.Normal.inverse_pdf
# )

# Base.println( ZigguratTable.Exponential.specialized_compute( 128, Base.MPFR.BigFloat, false ) )
# Base.println( ZigguratTable.Exponential.specialized_compute( 256, Base.MPFR.BigFloat, false ) )
# Base.println( ZigguratTable.Exponential.specialized_compute( 512, Base.MPFR.BigFloat, false ) )

# Base.println( ZigguratTable.Normal.specialized_compute( 128, Base.MPFR.BigFloat, false ) )
# Base.println( ZigguratTable.Normal.specialized_compute( 256, Base.MPFR.BigFloat, false ) )
# Base.println( ZigguratTable.Normal.specialized_compute( 512, Base.MPFR.BigFloat, false ) )

@timev r, v, z, flag = ZigguratTable.Normal.specialized_compute( 256, Core.Float64, false )
Base.println( r - Random.ziggurat_nor_r )

@timev r, v, z, flag = ZigguratTable.Exponential.specialized_compute( 256, Core.Float64, false )
Base.println( r - Random.ziggurat_exp_r )