<commit_before>#include "tinyspline.h"

#include <stdlib.h>
#include <math.h>
#include <string.h>

/********************************************************
*                                                       *
* System dependent configuration                        *
*                                                       *
********************************************************/
#define FLT_MAX_ABS_ERROR 1e-7
#define FLT_MAX_REL_ERROR 1e-5


/********************************************************
*                                                       *
* Internal functions                                    *
*                                                       *
********************************************************/
tsError ts_internal_bspline_find_u(
    const tsBSpline* bspline, const float u,
    size_t* k, size_t* s
)
{
    const size_t deg = bspline->deg;
    const size_t order = bspline->order;
    const size_t n_knots = bspline->n_knots;

    for (*k = 0; *k < n_knots; (*k)++) {
        const float uk = bspline->knots[*k];
        if (ts_fequals(u, uk)) {
            (*s)++;
        } else if (u < uk) {
            break;
        }
    }

    /* keep in mind that currently k is k+1 */
    if (*s > order)
        return TS_MULTIPLICITY;
    if (*k <= deg)                /* u < u_min */
        return TS_U_UNDEFINED;
    if (*k == n_knots && *s == 0) /* u > u_last */
        return TS_U_UNDEFINED;
    if (*k > n_knots-deg + *s-1)  /* u > u_max */
        return TS_U_UNDEFINED;

    (*k)--; /* k+1 - 1 will never underflow */
    return TS_SUCCESS;
}

/**
 * @return ::TS_SUCCESS, ::TS_MULTIPLICITY, ::TS_OVER_UNDERFLOW
 */
tsError ts_internal_bspline_insert_knot(
    const tsBSpline* bspline, const tsDeBoorNet* deBoorNet, const size_t n,
    tsBSpline* result
)
{
    tsError err;

    if (deBoorNet->s+n > bspline->order)
        goto err_multiplicity;

    err = ts_bspline_resize(bspline, (int)n, 1, result);
    if (err < 0)
        return err;
    if (n == 0) /* nothing to insert */
        return TS_SUCCESS;

    const size_t deg = bspline->deg;
    const size_t dim = bspline->dim;
    const size_t k = deBoorNet->k;
    const size_t N = deBoorNet->h+1; /* the number of affected ctrlps
 * n > 0 implies s <= deg implies a regular evaluation implies h+1 is valid */
    const size_t size_ctrlp = dim * sizeof(float);

    /* 1. Copy all necessary control points and knots from
     *    the original B-Spline.
     * 2. Copy all necessary control points and knots from
     *    the de Boor net. */

    /* 1.
     *
     * a) copy left hand side control points from original b-spline
     * b) copy right hand side control points from original b-spline
     * c) copy left hand side knots from original b-spline
     * d) copy right hand side knots form original b-spline */
    const size_t cidx = k-deg+N;
    const size_t kidx = k+1;
    /* copy control points */
    memmove(result->ctrlp, bspline->ctrlp, (k-deg) * size_ctrlp);
    memmove(
        result->ctrlp + (cidx+n)*dim, /* n >= 0 implies cidx+n >= cidx */
        bspline->ctrlp + cidx*dim,
        (result->n_ctrlp-n-cidx) * size_ctrlp
    );
    /* copy knots */
    memmove(result->knots, bspline->knots, (k+1) * sizeof(float));
    memmove(
        result->knots+kidx+n, /* n >= 0 implies kidx+n >= kidx */
        bspline->knots+kidx,
        (result->n_knots-n-kidx) * sizeof(float)
    );

    /* 2.
     *
     * a) copy left hand side control points from de boor net
     * b) copy middle part control points from de boor net
     * c) copy right hand side control points from de boor net
     * d) insert knots with u_k */
    size_t i; /* used in for loops */
    float* from = deBoorNet->points;
    float* to = result->ctrlp + (k-deg)*dim;
    int stride = (int)(N*dim); /* will be negative in c), thus use int */

    /* copy control points */
    for (i = 0; i < n; i++) { /* a) */
        memcpy(to, from, size_ctrlp);
        from += stride;
        to += dim;
        stride -= dim;
    }
    memcpy(to, from, (N-n) * size_ctrlp); /* b) */

    from -= dim;
    to += (N-n)*dim;
    stride = (int)(-(N-n+1)*dim);

    for (i = 0; i < n; i++) { /* c) */
        memcpy(to, from, size_ctrlp);
        from += stride;
        stride -= dim;
        to += dim;
    }
    /* copy knots */
    to = result->knots+k+1;
    for (i = 0; i < n; i++) { /* d) */
        *to = deBoorNet->u;
        to++;
    }

    return TS_SUCCESS;

    /* error handling */
    err_multiplicity:
        err = TS_MULTIPLICITY;
        goto cleanup;
    cleanup:
        if (bspline != result)
            ts_bspline_default(result);
        return err;
}

tsError ts_internal_bspline_thomas_algorithm(
    const float* points, const size_t n, const size_t dim,
    float* output
)
{
    if (points == output)
        return TS_INPUT_EQ_OUTPUT;
    if (dim == 0)
        return TS_DIM_ZERO;
    if (n == 0)
        return TS_DEG_GE_NCTRLP;

    const size_t size_flt = sizeof(float);
    const size_t ndsf = n*dim*size_flt;

    if (n <= 2) {
        memcpy(output, points, ndsf);
        return TS_SUCCESS;
    }

    /* in the following n >= 3 applies */
    size_t i, d; /* loop counter */
    size_t j, k, l; /* temporary indices */

    /* m_0 = 1/4, m_{k+1} = 1/(4-m_k), for k = 0,...,n-2 */
    const size_t len_m = n-2; /* n >= 3 implies n-2 >= 1 */
    float* m = malloc(len_m*size_flt);
    if (m == NULL)
        return TS_MALLOC;
    m[0] = 0.25f;
    for (i = 1; i < len_m+1; i++)
        m[i] = 1.f/(4 - m[i-1]);

    const size_t lst = (n-1)*dim; /* n >= 3 implies n-1 >= 2 */
    memset(output, 0, ndsf);
    memcpy(output+lst, points+lst, dim*size_flt);

    /* forward sweep */
    for (d = 0; d < dim; d++) {
        k = dim+d;
        output[k] = 6*points[k];
        output[k] -= points[d];
    }
    for (i = 2; i <= n-2; i++) {
        for (d = 0; d < dim; d++) {
            j = (i-1)*dim+d;
            k = i*dim+d;
            l = (i+1)*dim+d;
            output[k] = 6*points[k];
            output[k] -= output[l];
            output[k] -= m[i-2]*output[j];
        }
    }

    /* back substitution */
    memset(output+lst, 0, dim*size_flt);
    for (i = n-2; i >= 1; i--) {
        for (d = 0; d < dim; d++) {<commit_msg>fixed buffer overflow in interpolation<commit_after>#include "tinyspline.h"

#include <stdlib.h>
#include <math.h>
#include <string.h>

/********************************************************
*                                                       *
* System dependent configuration                        *
*                                                       *
********************************************************/
#define FLT_MAX_ABS_ERROR 1e-7
#define FLT_MAX_REL_ERROR 1e-5


/********************************************************
*                                                       *
* Internal functions                                    *
*                                                       *
********************************************************/
tsError ts_internal_bspline_find_u(
    const tsBSpline* bspline, const float u,
    size_t* k, size_t* s
)
{
    const size_t deg = bspline->deg;
    const size_t order = bspline->order;
    const size_t n_knots = bspline->n_knots;

    for (*k = 0; *k < n_knots; (*k)++) {
        const float uk = bspline->knots[*k];
        if (ts_fequals(u, uk)) {
            (*s)++;
        } else if (u < uk) {
            break;
        }
    }

    /* keep in mind that currently k is k+1 */
    if (*s > order)
        return TS_MULTIPLICITY;
    if (*k <= deg)                /* u < u_min */
        return TS_U_UNDEFINED;
    if (*k == n_knots && *s == 0) /* u > u_last */
        return TS_U_UNDEFINED;
    if (*k > n_knots-deg + *s-1)  /* u > u_max */
        return TS_U_UNDEFINED;

    (*k)--; /* k+1 - 1 will never underflow */
    return TS_SUCCESS;
}

/**
 * @return ::TS_SUCCESS, ::TS_MULTIPLICITY, ::TS_OVER_UNDERFLOW
 */
tsError ts_internal_bspline_insert_knot(
    const tsBSpline* bspline, const tsDeBoorNet* deBoorNet, const size_t n,
    tsBSpline* result
)
{
    tsError err;

    if (deBoorNet->s+n > bspline->order)
        goto err_multiplicity;

    err = ts_bspline_resize(bspline, (int)n, 1, result);
    if (err < 0)
        return err;
    if (n == 0) /* nothing to insert */
        return TS_SUCCESS;

    const size_t deg = bspline->deg;
    const size_t dim = bspline->dim;
    const size_t k = deBoorNet->k;
    const size_t N = deBoorNet->h+1; /* the number of affected ctrlps
 * n > 0 implies s <= deg implies a regular evaluation implies h+1 is valid */
    const size_t size_ctrlp = dim * sizeof(float);

    /* 1. Copy all necessary control points and knots from
     *    the original B-Spline.
     * 2. Copy all necessary control points and knots from
     *    the de Boor net. */

    /* 1.
     *
     * a) copy left hand side control points from original b-spline
     * b) copy right hand side control points from original b-spline
     * c) copy left hand side knots from original b-spline
     * d) copy right hand side knots form original b-spline */
    const size_t cidx = k-deg+N;
    const size_t kidx = k+1;
    /* copy control points */
    memmove(result->ctrlp, bspline->ctrlp, (k-deg) * size_ctrlp);
    memmove(
        result->ctrlp + (cidx+n)*dim, /* n >= 0 implies cidx+n >= cidx */
        bspline->ctrlp + cidx*dim,
        (result->n_ctrlp-n-cidx) * size_ctrlp
    );
    /* copy knots */
    memmove(result->knots, bspline->knots, (k+1) * sizeof(float));
    memmove(
        result->knots+kidx+n, /* n >= 0 implies kidx+n >= kidx */
        bspline->knots+kidx,
        (result->n_knots-n-kidx) * sizeof(float)
    );

    /* 2.
     *
     * a) copy left hand side control points from de boor net
     * b) copy middle part control points from de boor net
     * c) copy right hand side control points from de boor net
     * d) insert knots with u_k */
    size_t i; /* used in for loops */
    float* from = deBoorNet->points;
    float* to = result->ctrlp + (k-deg)*dim;
    int stride = (int)(N*dim); /* will be negative in c), thus use int */

    /* copy control points */
    for (i = 0; i < n; i++) { /* a) */
        memcpy(to, from, size_ctrlp);
        from += stride;
        to += dim;
        stride -= dim;
    }
    memcpy(to, from, (N-n) * size_ctrlp); /* b) */

    from -= dim;
    to += (N-n)*dim;
    stride = (int)(-(N-n+1)*dim);

    for (i = 0; i < n; i++) { /* c) */
        memcpy(to, from, size_ctrlp);
        from += stride;
        stride -= dim;
        to += dim;
    }
    /* copy knots */
    to = result->knots+k+1;
    for (i = 0; i < n; i++) { /* d) */
        *to = deBoorNet->u;
        to++;
    }

    return TS_SUCCESS;

    /* error handling */
    err_multiplicity:
        err = TS_MULTIPLICITY;
        goto cleanup;
    cleanup:
        if (bspline != result)
            ts_bspline_default(result);
        return err;
}

tsError ts_internal_bspline_thomas_algorithm(
    const float* points, const size_t n, const size_t dim,
    float* output
)
{
    if (points == output)
        return TS_INPUT_EQ_OUTPUT;
    if (dim == 0)
        return TS_DIM_ZERO;
    if (n == 0)
        return TS_DEG_GE_NCTRLP;

    const size_t size_flt = sizeof(float);
    const size_t ndsf = n*dim*size_flt;

    if (n <= 2) {
        memcpy(output, points, ndsf);
        return TS_SUCCESS;
    }

    /* in the following n >= 3 applies */
    size_t i, d; /* loop counter */
    size_t j, k, l; /* temporary indices */

    /* m_0 = 1/4, m_{k+1} = 1/(4-m_k), for k = 0,...,n-2 */
    const size_t len_m = n-2; /* n >= 3 implies n-2 >= 1 */
    float* m = malloc(len_m*size_flt);
    if (m == NULL)
        return TS_MALLOC;
    m[0] = 0.25f;
    for (i = 1; i < len_m; i++)
        m[i] = 1.f/(4 - m[i-1]);

    const size_t lst = (n-1)*dim; /* n >= 3 implies n-1 >= 2 */
    memset(output, 0, ndsf);
    memcpy(output+lst, points+lst, dim*size_flt);

    /* forward sweep */
    for (d = 0; d < dim; d++) {
        k = dim+d;
        output[k] = 6*points[k];
        output[k] -= points[d];
    }
    for (i = 2; i <= n-2; i++) {
        for (d = 0; d < dim; d++) {
            j = (i-1)*dim+d;
            k = i*dim+d;
            l = (i+1)*dim+d;
            output[k] = 6*points[k];
            output[k] -= output[l];
            output[k] -= m[i-2]*output[j];
        }
    }

    /* back substitution */
    memset(output+lst, 0, dim*size_flt);
    for (i = n-2; i >= 1; i--) {
        for (d = 0; d < dim; d++) {
{'fraction_non_alphanumeric': 0.1748028987815561, 'fraction_numerical': 0.011388070398980649, 'mean_word_length': 1.8599407879754042, 'pattern_counts': {'":': 0, '<': 39, '<?xml version=': 0, '>': 89, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 18, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '878287', 'n_tokens_mistral': 4540, 'n_tokens_neox': 4340, 'n_words': 1202}