Bessels.jl Acc test (20260113+849c4cc)

besseli_domains.py

# SPDX-License-Identifier: MIT
# Author: Chengyu Han (inkydragon/awesome-special-functions)
"""
Modified Bessel functions of the first kind of order nu:  besseli(nu, x)

- `nu == 0`:  besseli0(x)
- `nu == 1`:  besseli1(x)
- Neg `nu`: ``I_{-v}(x) = I_v(x) + (2/π)*sin(πv) * K_v(x)``
    - Special case: integer `v=n`: ``I_{-n}(x) = I_n(x)``
- Neg `x`:  ``I_v(-x) = e^{iπv}* I_v(x)``
    - Special case: integer `v=n`: ``I_n(-x) = (-1)^n * I_n(x)``

@assert nu > 0 && x > 0

besseli_positive_args:

if besseli_large_argument_cutoff(nu, x)
    # large argument expansion if x >> nu:   x > (nu^2 / 2 + 19.0)
    return besseli_large_args(nu, x)
elseif besselik_debye_cutoff(nu, x)
    # uniform debye expansion if x or nu is large:  nu > 25.0 || x > 35.0
    return besseli_large_orders(nu, x)
else
    # power series
    return besseli_power_series(nu, x)
end

besseli_large_args:
    Large Argument Expansion
    Hankel’s Expansions
    @assert abs(phase(z)) <= pi/2 - δ
    https://dlmf.nist.gov/10.40.E1

besseli_large_orders:
    Uniform Asymptotic Expansions (Debye)
    https://dlmf.nist.gov/10.41.3

besseli_power_series:
    Standard Solutions
    https://dlmf.nist.gov/10.25.E2

"""

# Plot domains and formula selection for besseli (positive arguments)
# Rules:
#   if x > (nu^2 / 2 + 19) -> Large-Argument (Hankel §10.40)
#   elif (nu > 25 or x > 35) -> Uniform Asymptotic (Debye §10.41)
#   else -> Power Series (§10.25)
#
# Note: main plot shows nu>0, x>0; branch rules for negatives are in text box.

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.colors import ListedColormap
from matplotlib.patches import Patch

# Range and resolution
NU_MAX = 50.0
X_MAX = 1000.0
N_NU = 400
N_X = 800

# Threshold and boundary functions
def large_argument_threshold(nu):
    return (nu**2) / 2.0 + 19.0

DEBYE_NU = 25.0
DEBYE_X = 35.0

# 区域选择编码
REG_POWER = 0
REG_DEBYE = 1
REG_LARGE = 2

def choose_region(nu, x):
    if x > large_argument_threshold(nu):
        return REG_LARGE
    elif (nu > DEBYE_NU) or (x > DEBYE_X):
        return REG_DEBYE
    else:
        return REG_POWER

def main():
    # 网格
    nus = np.linspace(0.0, NU_MAX, N_NU)
    xs = np.linspace(0.0, X_MAX, N_X)
    NU, X = np.meshgrid(nus, xs, indexing="xy")  # X: N_X x N_NU

    # 计算区域标签
    region = np.empty_like(NU, dtype=np.int32)
    # 向量化应用选择函数
    region = np.vectorize(choose_region)(NU, X)

    # Colors and legend
    cmap = ListedColormap([
        "#8FD14F",  # Power Series: green
        "#3A7BD5",  # Debye: blue
        "#E74C3C",  # Large-Argument: red
    ])
    labels = {
        REG_POWER: "Power Series (DLMF §10.25)",
        REG_DEBYE: "Uniform Asymptotic (Debye) (DLMF §10.41)",
        REG_LARGE: "Large-Argument (Hankel) (DLMF §10.40)",
    }

    fig, ax = plt.subplots(figsize=(12, 8), dpi=120)

    # Colored domains (semi-transparent for boundary visibility)
    im = ax.pcolormesh(nus, xs, region, cmap=cmap, shading="auto", alpha=0.35)

    # Boundary: x = nu^2/2 + 19
    x_la = large_argument_threshold(nus)
    mask_la = x_la <= X_MAX
    line_hankel, = ax.plot(nus[mask_la], x_la[mask_la], color="#C0392B", lw=2.0, label="x = ν²/2 + 19 (Hankel boundary)")

    # Boundaries: Debye thresholds (ν on x-axis, x on y-axis)
    vline_debye = ax.axvline(DEBYE_NU, color="#2E86C1", lw=2.2, ls="--", alpha=0.95, zorder=3, label="ν = 25 (Debye threshold)")
    hline_debye = ax.axhline(DEBYE_X, color="#884EA0", lw=2.2, ls="-.", alpha=0.95, zorder=3, label="x = 35 (Debye threshold)")

    # Sample points: annotate regions using color+marker
    rng = np.random.default_rng(42)
    def scatter_region(code, marker, color, n_samples=300):
        idx = np.argwhere(region == code)
        if idx.size == 0:
            return
        sel = idx[rng.choice(idx.shape[0], size=min(n_samples, idx.shape[0]), replace=False)]
        xs_sel = xs[sel[:, 0]]
        nus_sel = nus[sel[:, 1]]
        ax.scatter(nus_sel, xs_sel, s=10, marker=marker, color=color, alpha=0.35, linewidths=0)

    scatter_region(REG_POWER, marker="s", color="#2ECC71", n_samples=300)
    scatter_region(REG_DEBYE, marker="o", color="#3498DB", n_samples=300)
    scatter_region(REG_LARGE, marker="^", color="#E74C3C", n_samples=300)

    # Legend (color patches + boundaries)
    legend_patches = [
        Patch(facecolor="#8FD14F", edgecolor="none", alpha=0.35, label=labels[REG_POWER]),
        Patch(facecolor="#3A7BD5", edgecolor="none", alpha=0.35, label=labels[REG_DEBYE]),
        Patch(facecolor="#E74C3C", edgecolor="none", alpha=0.35, label=labels[REG_LARGE]),
    ]
    ax.legend(handles=legend_patches + [line_hankel, vline_debye, hline_debye], loc="upper left", framealpha=0.9)

    # Axis labels and ranges
    ax.set_xlabel("ν (order)", fontsize=12)
    ax.set_ylabel("x (argument)", fontsize=12)
    ax.set_xlim(0, NU_MAX)
    ax.set_ylim(0, X_MAX)

    ax.grid(True, ls="--", alpha=0.35)
    ax.set_title("besseli domains and formula selection (positive arguments)", fontsize=14)
    plt.tight_layout()
    plt.savefig("besseli_domains.png", dpi=120)
    plt.show()

if __name__ == "__main__":
    main()

besseli_relerr.jl

# SPDX-License-Identifier: MIT
# Author: Chengyu Han (inkydragon/awesome-special-functions)
# For BesselI(nu, x)
using Bessels
using ArbNumerics
using Plots
using Statistics
using Random
using Printf

setprecision(ArbFloat, 256)

# include("randf.jl")
# SPDX-License-Identifier: MIT OR Apache-2.0
# Copy from:  inkydragon/PureLibm.jl:test/utils/ranges.jl
"""
Generate random float numbers using `rand(UInt)`.

# Arguments
- `lo::Unsigned`: lower bound
- `hi::Unsigned`: upper bound
- `n::Int`: number of random values

# Returns
- `Vector{T}`: Generator for random float values in the range `[lo, hi]`
"""
function rand_float(xu_lo::T, xu_hi::T, n::Int) where {T<:Unsigned}
    @assert xu_lo < xu_hi
    FloatType = Base.floattype(T)
    xu_range = xu_lo:xu_hi
    xu_rand = rand(xu_range, n)
    f_rand = Iterators.map(xu->reinterpret(FloatType, xu), xu_rand)
    f_rand
end

"""
Generate random float numbers using `rand(UInt)`.

# Arguments
- `lo::AbstractFloat`: lower bound
- `hi::AbstractFloat`: upper bound
- `n::Int`: number of random values

# Returns
- `Vector{T}`: Generator for random float values in the range `[lo, hi]`
"""
function rand_float(lo::T, hi::T, n::Int) where {T<:AbstractFloat}
    @assert abs(lo) < abs(hi)
    UIntBaseType = Base.uinttype(T)
    xu_lo = reinterpret(UIntBaseType, lo)
    xu_hi = reinterpret(UIntBaseType, hi)
    rand_float(xu_lo, xu_hi, n)
end

"""
Generate random float numbers around `x` with integer radius `r`.

# Arguments
- `x::AbstractFloat`: center value
- `r::Int`: radius in integer representation (ULPs)
- `n::Int`: number of random values

# Returns
- `Vector{T}`: Generator for random float values
"""
function rand_float_at(x::T; r::Int, n::Int) where {T<:AbstractFloat}
    UIntType = Base.uinttype(T)
    xu = reinterpret(UIntType, x)
    r_u = convert(UIntType, r)
    
    # Check for underflow
    xu_lo = (xu < r_u) ? zero(UIntType) : xu - r_u
    
    # Check for overflow
    xu_hi = (xu > typemax(UIntType) - r_u) ? typemax(UIntType) : xu + r_u
    
    rand_float(xu_lo, xu_hi, n)
end
#= END randf.jl =#

const NUM_POINTS = 500_000
const X_START = eps(eps(1.0))
const X_END = 100.0

"""
SKIP:  For large `nu`, we need to start at a larger value to avoid underflow.
"""
function x_start(nu::Int)
    # if log10(nu+1) > 1.5
    #     return 1.0
    # end

    return X_START
end
function x_start(fname::String)
    @assert startswith(fname, "besseli")
    nu = parse(Int, fname[8:end])
    return x_start(nu)
end

function generate_test_points(n_points::Int, func_name::String)
    Xoshiro(42)
    x_vals = Vector{Float64}(undef, 0)
    sizehint!(x_vals, n_points)

    # [X_START, 1.0]
    n_small = 0
    # if func_name != "besseli90"
        n_small = n_points ÷ 10
        append!(x_vals, rand_float(X_START, 1.0, n_small))
    # end
    # [1.0, X_END]
    n_large = n_points - n_small
    append!(x_vals, rand_float(1.0, X_END, n_large))
    
    return sort!(x_vals)
end

function calculate_errors(x_vals::Vector{Float64}, test_fn::Function, mp_ref_fn::Function)
    n_points = length(x_vals)

    y_approx_all = test_fn.(x_vals)
    
    diffs = Vector{Float64}(undef, n_points)
    y_exacts = Vector{Float64}(undef, n_points)
    for i in 1:n_points
        x_mp = ArbFloat(x_vals[i])

        y_exact_mp = mp_ref_fn(x_mp)
        y_exact_f64 = Float64(y_exact_mp)
        diff_mp = ArbFloat(y_approx_all[i]) - y_exact_mp
        
        diffs[i] = Float64(diff_mp)
        y_exacts[i] = y_exact_f64
    end
    
    # rel_errors
    rel_errors = Vector{Float64}(undef, n_points)
    for i in 1:n_points
        if y_exacts[i] == 0.0
            rel_errors[i] = 0.0
        else
            rel_errors[i] = diffs[i] / y_exacts[i]
        end
    end

    return x_vals, rel_errors
end

function rel_stats(rel::Vector{Float64})
    abs_rel = abs.(rel)
    valid_mask = isfinite.(abs_rel)
    
    if !any(valid_mask)
        @warn "No valid points for relative error calculation"
        return (max=NaN, median=NaN, mean=NaN, p99=NaN)
    end
    
    vals = abs_rel[valid_mask]
    return (
        max = maximum(vals),
        median = median(vals),
        mean = mean(vals),
        p99 = quantile(vals, 0.99)
    )
end

function plot_relative_error(x, rel, func_name::String)
    # Max abs rel error
    stat = rel_stats(rel)
    if isnan(stat.max)
        return
    end

    # ---- Calculate xlims
    x_limit_min = floor(x_start(func_name)) - 0.5
    x_limit_max = ceil(X_END) + 0.5

    # ---- Calculate ylims
    # 10^-13.5 => 10^-13.0
    top_yy = floor(-log10(stat.max))
    # We don't care about the smallest errors
    # 10^-20.1 => 10^-21.0
    bottom_yy = ceil(-log10(stat.median)) + 1
    top_exp = top_yy
    bottom_exp = min(bottom_yy, 20)
    y_limits = (10.0^(-bottom_exp), 10.0^(-top_exp))

    title = @sprintf(
        """%s(x) (N=%g k points)
        (max=%.3g, p99=%.3g, mean=%.3g, median=%.3g)
        """,
        func_name, NUM_POINTS/1000,
        stat.max, stat.p99, stat.mean, stat.median)
    p = plot(
        size=(800,600),
        dpi=300,
        title=title,
        xlabel="x",
        ylabel="log10(|Relative Error|)",
        yscale=:log10,
        legend=:bottomright,
        grid=:true,
        gridalpha=0.5,
        minorgrid=:true,
        xlims=(x_limit_min, x_limit_max),
        # xticks=floor(x_limit_min):1.0:ceil(x_limit_max),
        ylims=y_limits
    )

    pos_mask = (rel .> 0) .& isfinite.(rel)
    neg_mask = (rel .< 0) .& isfinite.(rel)

    # Plot positive errors
    scatter!(p, x[pos_mask], rel[pos_mask],
        label="Pos Error",
        markersize=1.2,
        markerstrokewidth=0,
        color=:dodgerblue,
        alpha=0.6
    )

    # Plot negative errors as absolute values
    scatter!(p, x[neg_mask], abs.(rel[neg_mask]),
        label="Neg Error (abs)",
        markersize=1.2,
        markerstrokewidth=0,
        color=:darkorange,
        alpha=0.6
    )
    
    # Add horizontal line for max_rel
    hline!(p, [stat.max], linestyle=:dash, color=:red, label=@sprintf("max(|rel|) = %.2e", stat.max))
    # Add horizontal line for p99
    hline!(p, [stat.p99], linestyle=:solid, color=:forestgreen, linewidth=2.0, alpha=0.8,
        label=@sprintf("p99(|rel|) = %.2e", stat.p99))
    # Add horizontal line for mean
    hline!(p, [stat.mean], linestyle=:solid, color=:red, label=@sprintf("mean(|rel|) = %.2e", stat.mean))

    x_range = @sprintf("x=%g~%g", floor(x_start(func_name)), (X_END))
    npoints_k = NUM_POINTS ÷ 1000
    fname = "out/jl-$(func_name)-$(x_range)-N=$(npoints_k)k.png"
    savefig(p, fname)
    @info "Plot saved to $fname"

    # Show plot
    display(p)
end

function run_benchmark(test_fn, mp_ref_fn, func_name::String=repr(test_fn))
    if isempty(func_name)
        func_name = repr(test_fn)
    end

    @info @sprintf("Start calculating %s(x) errors for %g points...\n", func_name, NUM_POINTS)
    x_vals = generate_test_points(NUM_POINTS, func_name)
    t_start = time()
    x, rel = calculate_errors(x_vals, test_fn, mp_ref_fn)
    t_elapsed = time() - t_start
    @info @sprintf("Calc finished in %.2f seconds.\n", t_elapsed)
    
    # Print stats
    stats = rel_stats(rel)
    @info @sprintf("Rel max   = %.3e\n", stats.max)
    @info @sprintf("Rel p99   = %.3e\n", stats.p99)
    @info @sprintf("Rel mean  = %.3e\n", stats.mean)
    @info @sprintf("Rel median= %.3e\n", stats.median)
    
    # Plot relative error
    plot_relative_error(x, rel, func_name)
end


# --- Main Execution ---

function main()
    bench_list = [
        # Iv(z): modified Bessel function of the first kind
        (test = Bessels.besseli0, ref = x -> ArbNumerics.besseli(ArbReal("0"), ArbReal(x)), name="besseli0"),
        (test = Bessels.besseli1, ref = x -> ArbNumerics.besseli(ArbReal("1"), ArbReal(x)), name="besseli1"),
        (test = x -> Bessels.besseli(9, x), ref = x -> ArbNumerics.besseli(ArbReal("9"), ArbReal(x)), name="besseli9"),
        (test = x -> Bessels.besseli(24, x), ref = x -> ArbNumerics.besseli(ArbReal("24"), ArbReal(x)), name="besseli24"),
        (test = x -> Bessels.besseli(30, x), ref = x -> ArbNumerics.besseli(ArbReal("30"), ArbReal(x)), name="besseli30"),
        (test = x -> Bessels.besseli(90, x), ref = x -> ArbNumerics.besseli(ArbReal("90"), ArbReal(x)), name="besseli90"),
    ]

    for (test_fn, mp_ref_fn, name) in bench_list
        run_benchmark(test_fn, mp_ref_fn, name)
        # break  # DEBUG: gen one plot
    end
end

main()