/*
 * Polynomial interpolation template
 * Copyright 2019 Anton Khirnov <anton@khirnov.net>
 *
 * This program is free software: you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation, either version 3 of the License, or
 * (at your option) any later version.
 *
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program.  If not, see <http://www.gnu.org/licenses/>.
 */

#if CONTIGUOUS
# define CONT cont
#else
# define CONT generic
#endif

#define JOIN3(a, b, c) a ## _ ## b ## _ ## c
#define FUNC2(name, cont, stencil) JOIN3(name, cont, stencil)
#define FUNC(name) FUNC2(name, CONT, STENCIL)

static void
FUNC(interp2d_transfer_line)(double *dst, ptrdiff_t dst_len,
                             const double *src, ptrdiff_t src_stride,
                             const ptrdiff_t *idx_x, const double *fact_x, const double *fact_y
#if !CONTIGUOUS
                           , ptrdiff_t dst_stride0, ptrdiff_t src_stride0
# define SSTRIDE1 src_stride0
# define DSTRIDE1 dst_stride0
#else
# define SSTRIDE1 1
# define DSTRIDE1 1
#endif
                           )
{
    for (size_t x = 0; x < dst_len; x++) {
        const double *src_data = src + SSTRIDE1 * idx_x[x];

        double val = 0.0;

        for (ptrdiff_t idx0 = 0; idx0 < STENCIL; idx0++) {
            double tmp = 0.0;
            for (ptrdiff_t idx1 = 0; idx1 < STENCIL; idx1++)
                tmp += src_data[idx0 * src_stride + idx1 * SSTRIDE1] * fact_x[STENCIL * x + idx1];
            val += tmp * fact_y[idx0];
        }

        dst[x * DSTRIDE1] = val;
    }
}

static void
FUNC(interp1d_transfer_line)(double *dst, ptrdiff_t dst_len,
                             const double *src, const ptrdiff_t *idx_x, const double *fact_x
#if !CONTIGUOUS
                           , ptrdiff_t dst_stride0, ptrdiff_t src_stride0
#endif
                           )
{
    for (size_t x = 0; x < dst_len; x++) {
        const double *src_data = src + SSTRIDE1 * idx_x[x];

        double val = 0.0;

        for (ptrdiff_t idx = 0; idx < STENCIL; idx++)
            val += src_data[idx * SSTRIDE1] * fact_x[STENCIL * x + idx];

        dst[x * DSTRIDE1] = val;
    }
}

#undef SSTRIDE1
#undef DSTRIDE1
#undef CONT
#undef JOIN3
#undef FUNC2
#undef FUNC