#ifndef SHARE_FUNC_H
#define SHARE_FUNC_H

#include <stdlib.h>
#include <stddef.h>
#include <math.h>
#include <stdio.h>
#include <string.h>
/* 主网格：0-based -> 1D */
static inline size_t idx_ex(int i0, int j0, int k0, const int ex[3]) {
    const int ex1 = ex[0], ex2 = ex[1];
    return (size_t)i0 + (size_t)j0 * (size_t)ex1 + (size_t)k0 * (size_t)ex1 * (size_t)ex2;
}

/*
 * fh 对应 Fortran: fh(-1:ex1, -1:ex2, -1:ex3)
 * ord=2 => shift=1
 * iF/jF/kF 为 Fortran 索引（可为 -1,0,1..ex）
 */
static inline size_t idx_fh_F_ord2(int iF, int jF, int kF, const int ex[3]) {
    const int shift = 1;
    const int nx = ex[0] + 2;      // ex1 + ord
    const int ny = ex[1] + 2;

    const int ii = iF + shift;     // 0..ex1+1
    const int jj = jF + shift;     // 0..ex2+1
    const int kk = kF + shift;     // 0..ex3+1

    return (size_t)ii + (size_t)jj * (size_t)nx + (size_t)kk * (size_t)nx * (size_t)ny;
}

/*
 * fh 对应 Fortran: fh(-2:ex1, -2:ex2, -2:ex3)
 * ord=3 => shift=2
 * iF/jF/kF 是 Fortran 索引（可为负）
 */
static inline size_t idx_fh_F(int iF, int jF, int kF, const int ex[3]) {
    const int shift = 2;                 // ord=3 -> -2..ex
    const int nx = ex[0] + 3;            // ex1 + ord
    const int ny = ex[1] + 3;

    const int ii = iF + shift;           // 0..ex1+2
    const int jj = jF + shift;           // 0..ex2+2
    const int kk = kF + shift;           // 0..ex3+2

    return (size_t)ii + (size_t)jj * (size_t)nx + (size_t)kk * (size_t)nx * (size_t)ny;
}

/*
 * func:  (1..extc1, 1..extc2, 1..extc3)   1-based in Fortran
 * funcc: (-ord+1..extc1, -ord+1..extc2, -ord+1..extc3) in Fortran
 *
 * C 里我们把：
 *   func  视为 0-based: i0=0..extc1-1, j0=0..extc2-1, k0=0..extc3-1
 *   funcc 用“平移下标”存为一维数组：
 *     iF in [-ord+1..extc1]  -> ii = iF + (ord-1)  in [0..extc1+ord-1]
 *     总长度 nx = extc1 + ord
 *     同理 ny = extc2 + ord, nz = extc3 + ord
 */

static inline size_t idx_func0(int i0, int j0, int k0, const int extc[3]) {
    const int nx = extc[0], ny = extc[1];
    return (size_t)i0 + (size_t)j0 * (size_t)nx + (size_t)k0 * (size_t)nx * (size_t)ny;
}

static inline size_t idx_funcc_F(int iF, int jF, int kF, int ord, const int extc[3]) {
    const int shift = ord - 1;          // iF = -shift .. extc1
    const int nx = extc[0] + ord;       // [-shift..extc1] 共 extc1+ord 个
    const int ny = extc[1] + ord;

    const int ii = iF + shift;          // 0..extc1+shift
    const int jj = jF + shift;          // 0..extc2+shift
    const int kk = kF + shift;          // 0..extc3+shift

    return (size_t)ii + (size_t)jj * (size_t)nx + (size_t)kk * (size_t)nx * (size_t)ny;
}

/*
 * 等价于 Fortran:
 * funcc(1:extc1,1:extc2,1:extc3)=func
 * do i=0,ord-1
 *   funcc(-i,1:extc2,1:extc3) = funcc(i+1,1:extc2,1:extc3)*SoA(1)
 * enddo
 * do i=0,ord-1
 *   funcc(:,-i,1:extc3) = funcc(:,i+1,1:extc3)*SoA(2)
 * enddo
 * do i=0,ord-1
 *   funcc(:,:,-i) = funcc(:,:,i+1)*SoA(3)
 * enddo
 */
static inline void symmetry_bd_impl(int ord,
                 int shift,
                 const int extc[3],
                 const double *__restrict func,
                 double *__restrict funcc,
                 const double SoA[3])
{
    const int extc1 = extc[0], extc2 = extc[1], extc3 = extc[2];
    const int nx = extc1 + ord;
    const int ny = extc2 + ord;

    const size_t snx = (size_t)nx;
    const size_t splane = (size_t)nx * (size_t)ny;
    const size_t interior_i = (size_t)shift + 1u;          /* iF = 1 */
    const size_t interior_j = ((size_t)shift + 1u) * snx;  /* jF = 1 */
    const size_t interior_k = ((size_t)shift + 1u) * splane; /* kF = 1 */
    const size_t interior0 = interior_k + interior_j + interior_i;

    /* 1) funcc(1:extc1,1:extc2,1:extc3) = func */
    for (int k0 = 0; k0 < extc3; ++k0) {
        const double *src_k = func + (size_t)k0 * (size_t)extc2 * (size_t)extc1;
        const size_t dst_k0 = interior0 + (size_t)k0 * splane;
        for (int j0 = 0; j0 < extc2; ++j0) {
            const double *src = src_k + (size_t)j0 * (size_t)extc1;
            double *dst = funcc + dst_k0 + (size_t)j0 * snx;
            memcpy(dst, src, (size_t)extc1 * sizeof(double));
        }
    }

    /* 2) funcc(-i,1:extc2,1:extc3) = funcc(i+1,1:extc2,1:extc3)*SoA(1) */
    const double s1 = SoA[0];
    if (s1 == 1.0) {
        for (int ii = 0; ii < ord; ++ii) {
            const size_t dst_i = (size_t)(shift - ii);
            const size_t src_i = (size_t)(shift + ii + 1);
            for (int k0 = 0; k0 < extc3; ++k0) {
                const size_t kbase = interior_k + (size_t)k0 * splane + interior_j;
                for (int j0 = 0; j0 < extc2; ++j0) {
                    const size_t off = kbase + (size_t)j0 * snx;
                    funcc[off + dst_i] = funcc[off + src_i];
                }
            }
        }
    } else if (s1 == -1.0) {
        for (int ii = 0; ii < ord; ++ii) {
            const size_t dst_i = (size_t)(shift - ii);
            const size_t src_i = (size_t)(shift + ii + 1);
            for (int k0 = 0; k0 < extc3; ++k0) {
                const size_t kbase = interior_k + (size_t)k0 * splane + interior_j;
                for (int j0 = 0; j0 < extc2; ++j0) {
                    const size_t off = kbase + (size_t)j0 * snx;
                    funcc[off + dst_i] = -funcc[off + src_i];
                }
            }
        }
    } else {
        for (int ii = 0; ii < ord; ++ii) {
            const size_t dst_i = (size_t)(shift - ii);
            const size_t src_i = (size_t)(shift + ii + 1);
            for (int k0 = 0; k0 < extc3; ++k0) {
                const size_t kbase = interior_k + (size_t)k0 * splane + interior_j;
                for (int j0 = 0; j0 < extc2; ++j0) {
                    const size_t off = kbase + (size_t)j0 * snx;
                    funcc[off + dst_i] = funcc[off + src_i] * s1;
                }
            }
        }
    }

    /* 3) funcc(:,-j,1:extc3) = funcc(:,j+1,1:extc3)*SoA(2) */
    const double s2 = SoA[1];
    if (s2 == 1.0) {
        for (int jj = 0; jj < ord; ++jj) {
            const size_t dst_j = (size_t)(shift - jj) * snx;
            const size_t src_j = (size_t)(shift + jj + 1) * snx;
            for (int k0 = 0; k0 < extc3; ++k0) {
                const size_t kbase = interior_k + (size_t)k0 * splane;
                double *dst = funcc + kbase + dst_j;
                const double *src = funcc + kbase + src_j;
                for (int i = 0; i < nx; ++i) dst[i] = src[i];
            }
        }
    } else if (s2 == -1.0) {
        for (int jj = 0; jj < ord; ++jj) {
            const size_t dst_j = (size_t)(shift - jj) * snx;
            const size_t src_j = (size_t)(shift + jj + 1) * snx;
            for (int k0 = 0; k0 < extc3; ++k0) {
                const size_t kbase = interior_k + (size_t)k0 * splane;
                double *dst = funcc + kbase + dst_j;
                const double *src = funcc + kbase + src_j;
                for (int i = 0; i < nx; ++i) dst[i] = -src[i];
            }
        }
    } else {
        for (int jj = 0; jj < ord; ++jj) {
            const size_t dst_j = (size_t)(shift - jj) * snx;
            const size_t src_j = (size_t)(shift + jj + 1) * snx;
            for (int k0 = 0; k0 < extc3; ++k0) {
                const size_t kbase = interior_k + (size_t)k0 * splane;
                double *dst = funcc + kbase + dst_j;
                const double *src = funcc + kbase + src_j;
                for (int i = 0; i < nx; ++i) dst[i] = src[i] * s2;
            }
        }
    }

    /* 4) funcc(:,:,-k) = funcc(:,:,k+1)*SoA(3) */
    const double s3 = SoA[2];
    if (s3 == 1.0) {
        for (int kk = 0; kk < ord; ++kk) {
            const size_t dst_k = (size_t)(shift - kk) * splane;
            const size_t src_k = (size_t)(shift + kk + 1) * splane;
            double *dst = funcc + dst_k;
            const double *src = funcc + src_k;
            for (size_t p = 0; p < splane; ++p) dst[p] = src[p];
        }
    } else if (s3 == -1.0) {
        for (int kk = 0; kk < ord; ++kk) {
            const size_t dst_k = (size_t)(shift - kk) * splane;
            const size_t src_k = (size_t)(shift + kk + 1) * splane;
            double *dst = funcc + dst_k;
            const double *src = funcc + src_k;
            for (size_t p = 0; p < splane; ++p) dst[p] = -src[p];
        }
    } else {
        for (int kk = 0; kk < ord; ++kk) {
            const size_t dst_k = (size_t)(shift - kk) * splane;
            const size_t src_k = (size_t)(shift + kk + 1) * splane;
            double *dst = funcc + dst_k;
            const double *src = funcc + src_k;
            for (size_t p = 0; p < splane; ++p) dst[p] = src[p] * s3;
        }
    }
}

static inline void symmetry_bd(int ord,
                 const int extc[3],
                 const double *func,
                 double *funcc,
                 const double SoA[3])
{
    if (ord <= 0) return;

    /* Fast paths used by current C kernels: ord=2 (derivs), ord=3 (lopsided/KO). */
    if (ord == 2) {
        symmetry_bd_impl(2, 1, extc, func, funcc, SoA);
        return;
    }
    if (ord == 3) {
        symmetry_bd_impl(3, 2, extc, func, funcc, SoA);
        return;
    }

    symmetry_bd_impl(ord, ord - 1, extc, func, funcc, SoA);
}
#endif