[RFC]: add `math/base/tools/normhermitepolyf`

[kgryte]

Description

This RFC proposes adding the package @stdlib/math/base/tools/normhermitepolyf for evaluating a normalized Hermite polynomial using single-precision floating-point arithmetic.

Package: @stdlib/math/base/tools/normhermitepolyf Alias: normhermitepolyf

Related Issues

None.

Questions

No.

Other

This RFC has prior art. Namely,

• @stdlib/math/base/tools/normhermitepoly: uses double-precision floating-point arithmetic
• @stdlib/math/base/tools/evalpolyf: emulates single-precision floating-point arithmetic using float64ToFloat32.

This RFC proposes effectively copying the normhermitepoly package and modifying the implementation, tests, and examples (as in evalpolyf) to use float64ToFloat32 for emulating single-precision arithmetic.

Checklist

• [X] I have read and understood the Code of Conduct.
• [X] Searched for existing issues and pull requests.
• [X] The issue name begins with RFC:.

[xaman27x]

Would like to work upon this! @kgryte

[kgryte]

@xaman27x It looks like you've already been assigned https://github.com/stdlib-js/stdlib/issues/2018. I suggest moving that along first. If that is blocked, you can circle back here.

[xaman27x]

@xaman27x It looks like you've already been assigned https://github.com/stdlib-js/stdlib/issues/2018. I suggest moving that along first. If that is blocked, you can circle back here.

On it, Thank you!

[Daniel777y]

Hi @kgryte , can I work on this? I will try it. Thank you!

[kgryte]

@Daniel777y As @xaman27x is working on #2018, feel free to work on this issue. I'll go ahead and assign you.

[Daniel777y]

Hi @kgryte , I've implemented all the files but didn't pass the test (about 5% not passed). I guess I get stuck somewhere and need help.

I am not quite sure if it's because of the precision of test cases or of my implementation. Do I need to convert those cases to np.float32 in python?

Here's my implementation:

function normhermitepolyf( n, x ) {
    var y1;
    var y2;
    var y3;
    var i;

    x = float64ToFloat32( x );
    if ( isnan( n ) || isnanf( x ) || n < 0 || !isint( n ) ) {
        return NaN;
    }
    if ( n === 0 ) {
        // `x` is completely canceled from the expression:
        return 1.0;
    }
    if ( n === 1 ) {
        return x;
    }
    y2 = 1.0;
    y3 = 0.0;
    for ( i = n; i > 1; i-- ) {
        y1 = float64ToFloat32( float64ToFloat32( x * y2 ) - float64ToFloat32( i * y3 ) );
        y3 = y2;
        y2 = y1;
    }
    return float64ToFloat32( float64ToFloat32( x * y2 ) - y3 );
}

Here're the fixtures generator:

"""Generate fixtures."""

import os
import json
import numpy as np
from scipy.special import eval_hermitenorm

# Get the file path:
FILE = os.path.realpath(__file__)

# Extract the directory in which this file resides:
DIR = os.path.dirname(FILE)

def gen(n, x, name):
    """Generate fixture data and write to file.

    # Arguments

    * `n`: degree(s)
    * `x`: domain
    * `name::str`: output filename

    # Examples

    ``` python
    python> n = 1
    python> x = linspace(-1000, 1000, 2001)
    python> gen(n, x, './data.json')

"""
y = eval_hermitenorm(n, x)
if isinstance(n, np.ndarray):
    data = {
        "n": n.tolist(),
        "x": x.tolist(),
        "expected": y.tolist()
    }
else:
    data = {
        "n": n,
        "x": x.tolist(),
        "expected": y.tolist()
    }

# Based on the script directory, create an output filepath:
filepath = os.path.join(DIR, name)

# Write the data to the output filepath as JSON:
with open(filepath, "w", encoding="utf-8") as outfile:
    json.dump(data, outfile)

def main(): """Generate fixture data."""

Random values across n and x:

n = np.random.randint(1, 100, 1000)
x = np.float32(np.random.random(1000)*100.0)
gen(n, x, "random2.json")

# Medium negative:
x = np.linspace(-709.78, -1.0, 1000, dtype=np.float32)
gen(1, x, "medium_negative_1.json")
gen(2, x, "medium_negative_2.json")
gen(5, x, "medium_negative_5.json")

# Medium positive:
x = np.linspace(1.0, 709.78, 1000, dtype=np.float32)
gen(1, x, "medium_positive_1.json")
gen(2, x, "medium_positive_2.json")
gen(5, x, "medium_positive_5.json")

# Small positive:
x = np.linspace(2.0**-51, 1.0, 1000, dtype=np.float32)
gen(1, x, "small_positive_1.json")
gen(2, x, "small_positive_2.json")
gen(5, x, "small_positive_5.json")

# Small negative:
x = np.linspace(-1.0, -2.0**-51, 1000, dtype=np.float32)
gen(1, x, "small_negative_1.json")
gen(2, x, "small_negative_2.json")
gen(5, x, "small_negative_5.json")

# Tiny values:
x = np.linspace(-2.0**-51, 2.0**-51, 1000, dtype=np.float32)
gen(1, x, "tiny_1.json")
gen(2, x, "tiny_2.json")
gen(5, x, "tiny_5.json")

if name == "main": main()

Here's a unit of tests:
```javascript
tape( 'the function accurately computes `He1(x)` for small positive numbers', function test( t ) {
    var expected;
    var delta;
    var tol;
    var x;
    var v;
    var i;
    var n;
    var e;

    n = smallPositive1.n;
    x = smallPositive1.x;
    expected = smallPositive1.expected;

    for ( i = 0; i < x.length; i++ ) {
        v = normhermitepoly( n, x[ i ] );
        // I also convert the expected value to float32 here
        e = float64ToFloat32( expected[ i ] );
        delta = abs( v - e );
        tol = EPS * abs( e );
        t.ok( delta <= tol, 'within tolerance. n: ' + n + '. x: ' + x[ i ] + '. Value: ' + v + '. Expected: ' + expected[ i ] + '. Delta: ' + delta + '. Tolerance: ' + tol + '.' );
    }
    t.end();
});

If I need to include more information or there're other examples for reference, please let me know.

Thank you!

[kgryte]

@Daniel777y You can tweak the tolerance level a bit. E.g., 2.0 * EPS. That is fine.

[kgryte]

@Daniel777y Your implementation looks sound. It's possible that SciPy upcasts inputs to float64, resulting in higher intermediate precision. As long as the accuracy is not significantly far off, we can live with some minor divergence.

[Daniel777y]

Thank you for response @kgryte ! After tweaking the tolerance level, it's better.

But because of float32, some test cases would be large (or small) enough to be Infinity (or -Infinity). The calculation contains subtraction, so the return value can be NaN. Do I need to make sure that it will not return NaN?

Here're some potential ideas to pass tests:

1. Don't convert middle variables, only convert the input and output to float32. Then if the return value and expected value are both Infinity, accept it.
2. Adjust the test functions, try to figure out the whether the NaN is due to Infinity subtraction, and then accept those with expected value of Infinity.
3. Simply not include large and small cases temporarily.

But for now I am not sure if these methods can accurately calculate the results.

I also walked through the packages evalpolyf and evalrationalf. But the range of the test cases are kinda small, so they might not include the Infinity and NaN cases.

Or maybe there're other good ideas to deal with this problem, so I am looking forward to your suggestions.