#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;
}

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

/*
 * fh 对应 Fortran: fh(-3:ex1, -3:ex2, -3:ex3)
 * ord=4 => shift=3
 */
static inline size_t idx_fh_F_ord4(int iF, int jF, int kF, const int ex[3]) {
    const int shift = 3;
    const int nx = ex[0] + 4;            // ex1 + ord
    const int ny = ex[1] + 4;
    const int ii = iF + shift;           // 0..ex1+3
    const int jj = jF + shift;           // 0..ex2+3
    const int kk = kF + shift;           // 0..ex3+3
    return (size_t)ii + (size_t)jj * (size_t)nx + (size_t)kk * (size_t)nx * (size_t)ny;
}

/*
 * fh 对应 Fortran: fh(-4:ex1, -4:ex2, -4:ex3)
 * ord=5 => shift=4
 */
static inline size_t idx_fh_F_ord5(int iF, int jF, int kF, const int ex[3]) {
    const int shift = 4;
    const int nx = ex[0] + 5;            // ex1 + ord
    const int ny = ex[1] + 5;
    const int ii = iF + shift;           // 0..ex1+4
    const int jj = jF + shift;           // 0..ex2+4
    const int kk = kF + shift;           // 0..ex3+4
    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;

    if (ord == 1) {
        symmetry_bd_impl(1, 0, extc, func, funcc, SoA);
        return;
    }
    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;
    }
    if (ord == 4) {
        symmetry_bd_impl(4, 3, extc, func, funcc, SoA);
        return;
    }

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

/*
 * symmetry_stbd — shell-patch (staggered boundary) ghost fill.
 *
 * Fortran: funcc(-ord+1:extc1+ord, -ord+1:extc2+ord, extc3)
 * Only 2 SoA values (x/y). No z symmetry fill.
 * Ghost on BOTH positive and negative sides of x and y.
 * Reflection uses i+2 (skips boundary) instead of i+1.
 * nx = extc1 + 2*ord, ny = extc2 + 2*ord
 */
static inline void symmetry_stbd(int ord,
                                 const int extc[3],
                                 const double *func,
                                 double *funcc,
                                 const double SoA[2])
{
    const int extc1 = extc[0], extc2 = extc[1], extc3 = extc[2];
    const int nx = extc1 + 2 * ord;
    const int ny = extc2 + 2 * ord;
    const int sh = ord - 1;
    const size_t snx = (size_t)nx;
    const size_t splane = snx * (size_t)ny;

    /* 1) Copy interior: funcc(1:extc1, 1:extc2, 1:extc3) = func */
    for (int k0 = 0; k0 < extc3; ++k0) {
        const double *src = func + (size_t)k0 * (size_t)extc2 * (size_t)extc1;
        const size_t kbase = (size_t)k0 * splane;
        for (int j0 = 0; j0 < extc2; ++j0) {
            double *dst = funcc + kbase + (size_t)(sh + j0 + 1) * snx + (size_t)(sh + 1);
            const double *s = src + (size_t)j0 * (size_t)extc1;
            for (int i0 = 0; i0 < extc1; ++i0) dst[i0] = s[i0];
        }
    }

    /* 2) x-direction ghost fill */
    const double s1 = SoA[0];
    for (int k0 = 0; k0 < extc3; ++k0) {
        const size_t kbase = (size_t)k0 * splane;
        for (int j0 = 0; j0 < extc2; ++j0) {
            const size_t off = kbase + (size_t)(sh + j0 + 1) * snx;
            /* left side: funcc(-i) = funcc(i+2) * s1 */
            for (int i = 0; i < ord; ++i) {
                funcc[off + (size_t)(sh - i)]       = funcc[off + (size_t)(sh + i + 2)] * s1;
                /* right side: funcc(extc1+1+i) = funcc(extc1-1-i) * s1 */
                funcc[off + (size_t)(sh + extc1 + 1 + i)] = funcc[off + (size_t)(sh + extc1 - 1 - i)] * s1;
            }
        }
    }

    /* 3) y-direction ghost fill */
    const double s2 = SoA[1];
    for (int i = 0; i < nx; ++i) {
        for (int k0 = 0; k0 < extc3; ++k0) {
            const size_t kbase = (size_t)k0 * splane;
            /* bottom: funcc(:,-i,:) = funcc(:,i+2,:) * s2 */
            for (int jj = 0; jj < ord; ++jj) {
                funcc[kbase + (size_t)(sh - jj) * snx + (size_t)i] =
                    funcc[kbase + (size_t)(sh + jj + 2) * snx + (size_t)i] * s2;
                /* top: funcc(:,extc2+1+jj,:) = funcc(:,extc2-1-jj,:) * s2 */
                funcc[kbase + (size_t)(sh + extc2 + 1 + jj) * snx + (size_t)i] =
                    funcc[kbase + (size_t)(sh + extc2 - 1 - jj) * snx + (size_t)i] * s2;
            }
        }
    }
}

/*
 * Indexing for shell fh buffer: Fortran fh(-ord+1:extc1+ord, -ord+1:extc2+ord, extc3)
 * C 0-based: ii = iF + ord - 1
 * nx = extc1 + 2*ord, ny = extc2 + 2*ord
 */
static inline size_t idx_fh_stbd(int iF, int jF, int kF, int ord, const int extc[3]) {
    const int sh = ord - 1;
    const int nx = extc[0] + 2 * ord;
    const int ny = extc[1] + 2 * ord;
    const int ii = iF + sh;
    const int jj = jF + sh;
    const int kk = kF - 1;  // Fortran 1-based kF → C 0-based
    return (size_t)ii + (size_t)jj * (size_t)nx + (size_t)kk * (size_t)nx * (size_t)ny;
}
#endif