view release on metacpan or  search on metacpan

t/lib/TestVec.pm  view on Meta::CPAN

    my $eps = 1.0;
    $eps /= 2 while (1.0 + $eps/2) != 1.0;
    $NV_INFO{machine_epsilon} = $eps;
}

sub nv_info    { return \%NV_INFO }
sub nv_epsilon { return $NV_INFO{machine_epsilon} }
sub nv_digits  { return $NV_INFO{digits} }

sub is_quadmath    { return $NV_INFO{is_quadmath} }
sub is_long_double { return $NV_INFO{is_long_double} }

# ============================================
# Tolerance Calculation
# ============================================

# Get appropriate absolute tolerance based on NV precision
sub get_tolerance {
    my ($strict) = @_;
    if ($NV_INFO{is_quadmath}) {
        return $strict ? 1e-30 : 1e-12;
    } elsif ($NV_INFO{is_long_double}) {
        return $strict ? 1e-18 : 1e-12;
    }
    return $strict ? 1e-14 : 1e-9;
}

# Get relative tolerance (scales with magnitude)
sub relative_tolerance {
    my ($ulps) = @_;
    $ulps //= 4;  # Default: 4 ULPs
    return $NV_INFO{machine_epsilon} * $ulps;
}

# ============================================
# Comparison Methods
# ============================================

# 1. Absolute tolerance (original method)
sub approx_eq {
    my ($got, $expected, $tolerance) = @_;
    $tolerance //= get_tolerance();
    return 1 if !defined $got && !defined $expected;
    return 0 if !defined $got || !defined $expected;
    return abs($got - $expected) < $tolerance;
}

# 2. Relative epsilon comparison (handles different magnitudes)
sub relatively_equal {
    my ($got, $expected, $ulps) = @_;
    $ulps //= 4;

    return 1 if !defined $got && !defined $expected;
    return 0 if !defined $got || !defined $expected;

    # Handle exact zero
    return $got == 0 if $expected == 0;
    return $expected == 0 if $got == 0;

    # Handle infinities
    return $got == $expected if $got != $got || $expected != $expected;  # NaN

    # Relative comparison
    my $max_abs = abs($got) > abs($expected) ? abs($got) : abs($expected);
    my $diff = abs($got - $expected);
    return $diff <= $NV_INFO{machine_epsilon} * $max_abs * $ulps;
}

# 3. ULP (Units in Last Place) comparison - most rigorous
#    Counts how many representable floats apart two values are
sub ulp_distance {
    my ($a, $b) = @_;

    return 0 if !defined $a && !defined $b;
    return -1 if !defined $a || !defined $b;  # Error indicator

    # Handle special cases
    return 0 if $a == $b;  # Exact match (handles infinities)
    return -1 if $a != $a || $b != $b;  # NaN

    # For ULP calculation, we use the difference scaled by epsilon
    # This is a portable approximation that works across NV types
    my $max_abs = abs($a) > abs($b) ? abs($a) : abs($b);
    $max_abs = 1.0 if $max_abs == 0;  # Prevent division by zero

    my $ulps = abs($a - $b) / ($NV_INFO{machine_epsilon} * $max_abs);
    return int($ulps + 0.5);  # Round to nearest integer
}

sub ulp_equal {
    my ($got, $expected, $max_ulps) = @_;
    $max_ulps //= 4;

    my $dist = ulp_distance($got, $expected);
    return 0 if $dist < 0;  # Error case (NaN, undef)
    return $dist <= $max_ulps;
}

# 4. Exact bit comparison - no tolerance
sub bits_equal {
    my ($a, $b) = @_;
    return 0 if !defined $a || !defined $b;
    return pack("F", $a) eq pack("F", $b);
}

sub bits_hex {
    my ($val) = @_;
    return 'undef' unless defined $val;
    return unpack("H*", pack("F", $val));
}

# ============================================
# Test::More Integration
# ============================================

# Primary test function - uses relative comparison by default
sub float_is {
    my ($got, $expected, $name, $opts) = @_;
    $opts //= {};

    my $method = $opts->{method} // 'relative';
    my $ulps = $opts->{ulps} // 4;
    my $tolerance = $opts->{tolerance};

    my $ok;
    my $details = '';

    if ($method eq 'exact') {
        $ok = bits_equal($got, $expected);
        $details = sprintf("got bits: %s, expected bits: %s",
                          bits_hex($got), bits_hex($expected));
    }
    elsif ($method eq 'ulp') {
        my $dist = ulp_distance($got, $expected);
        $ok = $dist >= 0 && $dist <= $ulps;
        $details = sprintf("ULP distance: %d (max allowed: %d)", $dist, $ulps);
    }
    elsif ($method eq 'absolute') {
        $tolerance //= get_tolerance();
        $ok = approx_eq($got, $expected, $tolerance);
        $details = sprintf("diff: %g, tolerance: %g",
                          abs(($got // 0) - ($expected // 0)), $tolerance);
    }
    else {  # relative (default)
        $ok = relatively_equal($got, $expected, $ulps);
        if (defined $got && defined $expected && $expected != 0) {
            my $rel_err = abs($got - $expected) / abs($expected);
            $details = sprintf("relative error: %g, max allowed: %g",
                              $rel_err, $NV_INFO{machine_epsilon} * $ulps);
        }
    }

    local $Test::Builder::Level = $Test::Builder::Level + 1;
    ok($ok, $name);