#include "integrate.h"

static inline int bad_args(integrand_f f, float *a, float *b, unsigned *n,
                           float *out)
{
    return (f == NULL || out == NULL || *n == 0u || !(*b > *a));
}

integ_status_t midpoint(integrand_f f, void *ctx, float a, float b, unsigned n,
                        float *out)
{
    /*check if all necessary parameters have been given */
    if (bad_args(f, &a, &b, &n, out))
        return INTEG_ERR_BAD_ARGS;

    float sum = 0.0f;
    float fx = 0.0f;

    /*intervall width h: */
    const float h = (b - a) / (float)n;
    float x = a + 0.5f * h;

    for (unsigned i = 0; i < n; ++i)
    {
        fx = f(x, ctx);
        sum = sum + fx;
        x = x + h;
    }

    *out = sum * h;
    return INTEG_OK;
}

integ_status_t trapezoid(integrand_f f, void *ctx, float a, float b, unsigned n,
                         float *out)
{
    /*check if all necessary parameters have been given */
    if (bad_args(f, &a, &b, &n, out))
    {
        return INTEG_ERR_BAD_ARGS;
    }

    /*intervall width h: */
    const float h = (b - a) / (float)n;
    float sum = 0.5f * (f(a, ctx) + f(b, ctx));

    float x = a + h;
    for (unsigned i = 1; i < n; ++i)
    {
        sum = sum + f(x, ctx);
        x = x + h;
    }

    *out = sum * h;
    return INTEG_OK;
}

integ_status_t simpson(integrand_f f, void *ctx, float a, float b, unsigned n,
                       float *out)
{
    /*check if all necessary parameters have been given */
    if (bad_args(f, &a, &b, &n, out))
    {
        return INTEG_ERR_BAD_ARGS;
    }
    if (n & 1u) /*n must be even*/
    {
        return INTEG_ERR_N_EVEN_REQUIRED;
    }

    float sum_odd = 0.0f;
    float sum_even = 0.0f;
    /*intervall width h: */
    const float h = (b - a) / (float)n;

    float x = a + h;
    for (unsigned i = 1; i < n; ++i)
    {
        const float fx = f(x, ctx);
        if (i & 1u)
            sum_odd = sum_odd + fx;
        else
            sum_even = sum_even + fx;
        x = x + h;
    }

    *out =
        (h / 3.0f) * (f(a, ctx) + 4.0f * sum_odd + 2.0f * sum_even + f(b, ctx));
    return INTEG_OK;
}