From 3b895a3c12cabc60449755fbcd1bdc2813cde74c Mon Sep 17 00:00:00 2001
From: Timo Kreuzer <timo.kreuzer@reactos.org>
Date: Mon, 21 Apr 2025 17:39:19 +0300
Subject: [PATCH] [CRT_APITEST] Add a python tool to generate test tables for
 math functions

Note: The code contains a custom implementation of ldexp in python, because on Windows (below Win 11) ldexp is broken and rounds denormals incorrectly. ldexp is used by mpmath to convert multi-precision to float, which led to incorrect floating point results (specifically denormals).
---
 .../rostests/apitests/crt/gen_math_tests.py   | 154 ++++++++++++++++++
 1 file changed, 154 insertions(+)
 create mode 100644 modules/rostests/apitests/crt/gen_math_tests.py

diff --git a/modules/rostests/apitests/crt/gen_math_tests.py b/modules/rostests/apitests/crt/gen_math_tests.py
new file mode 100644
index 00000000000..fc6528f9f56
--- /dev/null
+++ b/modules/rostests/apitests/crt/gen_math_tests.py
@@ -0,0 +1,154 @@
+#!/usr/bin/env python3
+"""
+ PROJECT:     ReactOS CRT
+ LICENSE:     MIT (https://spdx.org/licenses/MIT)
+ PURPOSE:     Script to generate test data tables for math functions
+ COPYRIGHT:   Copyright 2025 Timo Kreuzer <timo.kreuzer@reactos.org>
+"""
+
+from mpmath import mp
+import struct
+import argparse
+import math
+import array
+
+# Set precision (e.g., 100 decimal places)
+mp.dps = 100
+
+def ldexp_manual(x_float, exp):
+    """
+    Compute x_float * 2**exp for a floating-point number, handling denormals and rounding.
+    This is a workaround for broken ldexp on Windows (before Win 11), which does not handle denormals correctly.
+
+    Args:
+        x_float (float): The floating-point number to scale.
+        exp (int): The integer exponent of 2 by which to scale.
+
+    Returns:
+        float: The result of x_float * 2**exp, respecting IEEE 754 rules.
+    """
+    # Handle special cases: zero, infinity, or NaN
+    if x_float == 0.0 or not math.isfinite(x_float):
+        return x_float
+
+    # Get the 64-bit IEEE 754 representation
+    bits = struct.unpack('Q', struct.pack('d', x_float))[0]
+
+    # Extract components
+    sign = bits >> 63
+    exponent = (bits >> 52) & 0x7FF
+    mantissa = bits & 0xFFFFFFFFFFFFF
+
+    if exponent == 0:
+        # Denormal number
+        if mantissa == 0:
+            return x_float  # Zero
+
+        # Normalize it
+        shift_amount = 52 - mantissa.bit_length()
+        mantissa <<= shift_amount
+        exponent = 1 - shift_amount
+    else:
+        # Normal number, add implicit leading 1
+        mantissa |= (1 << 52)
+
+    # Adjust exponent with exp
+    new_exponent = exponent + exp
+
+    if new_exponent > 2046:
+        # Overflow to infinity
+        return math.copysign(math.inf, x_float)
+    elif new_exponent <= 0:
+        # Denormal or underflow
+        if new_exponent < -52:
+            # Underflow to zero
+            return 0.0
+
+        # Calculate how much we need to shift the mantissa
+        mantissa_shift = 1 - new_exponent
+
+        # Get the highest bit, that would be shifted off
+        round_bit = (mantissa >> (mantissa_shift - 1)) & 1
+
+        # Adjust mantissa for denormal
+        mantissa >>= mantissa_shift
+        mantissa += round_bit
+        new_exponent = 0
+
+    # Reconstruct the float
+    bits = (sign << 63) | (new_exponent << 52) | (mantissa & 0xFFFFFFFFFFFFF)
+    return float(struct.unpack('d', struct.pack('Q', bits))[0])
+
+def mpf_to_float(mpf_value):
+    """
+    Convert an mpmath mpf value to the nearest Python float.
+    We use ldexp_manual, because ldexp is broken on Windows.
+    """
+    sign = mpf_value._mpf_[0]
+    mantissa = mpf_value._mpf_[1]
+    exponent = mpf_value._mpf_[2]
+
+    result = ldexp_manual(mantissa, exponent)
+    if sign == 1:
+        result = -result
+    return result
+
+def generate_range_of_floats(start, end, steps, typecode):
+    if not isinstance(steps, int) or steps < 1:
+        raise ValueError("steps must be an integer >= 1")
+    if steps == 1:
+        return [start]
+    step = (end - start) / (steps - 1)
+    #return [float(start + i * step) for i in range(steps)]
+    return array.array(typecode, [float(start + i * step) for i in range(steps)])
+
+def gen_table_header(func_name):
+    print("static TESTENTRY_DBL_APPROX s_" f"{func_name}" "_approx_tests[] =")
+    print("{")
+    print("//  {    x,                     {    y_rounded,               y_difference           } }")
+
+def gen_table_range(func_name, typecode, func, start, end, steps, max_ulp):
+
+    angles = generate_range_of_floats(start, end, steps, typecode)
+
+    for x in angles:
+        # Compute high-precision cosine
+        high_prec = func(x)
+
+        # Convert to double (rounds to nearest double)
+        rounded_double = mpf_to_float(high_prec)
+
+        # Compute difference
+        difference = high_prec - mp.mpf(rounded_double)
+
+        # Print the table line
+        print("    { " f"{float(x).hex():>24}"
+              ", { " f"{float(rounded_double).hex():>24}"
+              ", " f"{float(difference).hex():>24}"
+              " }, " f"{max_ulp}"
+              " }, // " f"{func_name}" "(" f"{float(x)}" ") == " f"{mp.nstr(high_prec, 20)}")
+
+
+# Dictionary to map math function names to generator functions
+TABLE_FUNCTIONS = {
+}
+
+def main():
+    # Set up argument parser
+    parser = argparse.ArgumentParser(description="Output a specific test table based on the provided function name.")
+    parser.add_argument("function_name", help="Name of the function to output the table for (e.g., sin, cos)")
+    args = parser.parse_args()
+
+    # Get the table name from the command-line argument
+    function_name = args.function_name.lower()
+
+    # Check if the table name is valid
+    if function_name in TABLE_FUNCTIONS:
+        # Call the corresponding function
+        TABLE_FUNCTIONS[function_name]()
+    else:
+        print(f"Error: Unsupported function '{function_name}'. Available tables: {', '.join(TABLE_FUNCTIONS.keys())}")
+        exit(1)
+
+if __name__ == "__main__":
+    main()