#!/usr/bin/env python3.10

"""Approximate the standard normal cumulative distribution function."""

import numpy as np


def normal_pdf(sample_points):
    """
    Evaluate the standard normal probability density function.

    Parameters
    ----------
    sample_points : float or numpy.ndarray
        Evaluation point or points.
    """
    return (1.0 / np.sqrt(2.0 * np.pi)) * np.exp(-(sample_points * sample_points) / 2.0)


def composite_simpson_rule(lower_bound, upper_bound, subinterval_count):
    """
    Approximate an integral with the composite Simpson 1/3 rule.

    Parameters
    ----------
    lower_bound, upper_bound : float
        Integration bounds.
    subinterval_count : int
        Number of subintervals. It should be even for Simpson's rule.
    """
    if lower_bound == upper_bound:
        return 0.0

    step_size = (upper_bound - lower_bound) / subinterval_count
    grid_points = np.linspace(lower_bound, upper_bound, subinterval_count + 1)
    function_values = normal_pdf(grid_points)

    integral = function_values[0] + function_values[-1]
    integral += 4.0 * np.sum(function_values[1:-1:2])
    integral += 2.0 * np.sum(function_values[2:-1:2])
    return step_size * integral / 3.0


def normal_cdf(z_value):
    """
    Approximate the standard normal cumulative distribution function.

    Parameters
    ----------
    z_value : float
        Point where the CDF should be evaluated.
    """
    upper_bound = abs(z_value)
    probability = 0.5 + composite_simpson_rule(0.0, upper_bound, 100)
    if z_value < 0.0:
        probability = 1.0 - probability
    return probability


def main():
    """Run the sample computation used in the original project."""
    z_value = 0.3
    probability = normal_cdf(z_value)
    print(probability)
    return probability


if __name__ == "__main__":
    main()