[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

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

[xaman27x]

Would like to work upon this! @kgryte

[kgryte]

@xaman27x It looks like you've already been assigned #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 #2018. I suggest moving that along first. If that is blocked, you can circle back here.

On it, Thank you!

[C0ldSmi1e]

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.

[C0ldSmi1e]

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:

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.