#!/usr/bin/env python3.10

"""Compute machine epsilon for single and double precision."""

import numpy as np


def machine_epsilon(float_type):
    """
    Compute machine epsilon for a given floating-point type.

    Parameters
    ----------
    float_type : type
        Floating-point type, such as `numpy.float32` or `numpy.float64`.

    Returns
    -------
    float
        Smallest positive number `epsilon` such that `1 + epsilon > 1`
        in the chosen floating-point arithmetic.
    """
    one_value = float_type(1.0)
    epsilon = float_type(1.0)

    while one_value + float_type(epsilon / 2.0) > one_value:
        epsilon = float_type(epsilon / 2.0)

    return float(epsilon)


def main():
    """Compute machine epsilon in single and double precision."""
    single_precision_epsilon = machine_epsilon(np.float32)
    double_precision_epsilon = machine_epsilon(np.float64)

    print(f"single_precision_epsilon = {single_precision_epsilon:.16e}")
    print(f"double_precision_epsilon = {double_precision_epsilon:.16e}")

    return single_precision_epsilon, double_precision_epsilon


if __name__ == "__main__":
    main()