#include "macrodef.h"
#include "tool.h"

/*
 * C 版 fdderivs — second derivatives d2f/dx2, d2f/dxdy, d2f/dxdz, d2f/dy2, d2f/dydz, d2f/dz2.
 *
 * Finite difference order selected at compile time via ghost_width macro.
 * Multi-pass skip strategy: lowest-order computes widest region while skipping
 * the union of higher-order regions, then each higher pass overwrites its interior.
 */
void fdderivs(const int ex[3],
              const double *f,
              double *fxx, double *fxy, double *fxz,
              double *fyy, double *fyz, double *fzz,
              const double *X, const double *Y, const double *Z,
              double SYM1, double SYM2, double SYM3,
              int Symmetry, int onoff)
{
    (void)onoff;

    const int NO_SYMM = 0, EQ_SYMM = 1;
    const double ZEO = 0.0, ONE = 1.0, TWO = 2.0;
    const double F1o4   = 2.5e-1;
    const double F8     = 8.0;
    const double F16    = 16.0;
    const double F30    = 30.0;
    const double F1o12  = ONE / 12.0;
    const double F1o144 = ONE / 144.0;
    const double F9  = 9.0,  F45 = 45.0,  F60 = 60.0;
    const double F27 = 27.0, F270 = 270.0, F490 = 490.0;
    const double F1o180  = ONE / 180.0;
    const double F1o3600 = ONE / 3600.0;
    const double F32 = 32.0, F128 = 128.0, F168 = 168.0, F672 = 672.0;
    const double F840  = 840.0, F1008 = 1008.0, F8064 = 8064.0, F14350 = 14350.0;
    const double F1o5040   = ONE / 5040.0;
    const double F1o705600 = ONE / 705600.0;

    const int ex1 = ex[0], ex2 = ex[1], ex3 = ex[2];
    const double dX = X[1] - X[0];
    const double dY = Y[1] - Y[0];
    const double dZ = Z[1] - Z[0];

    const int imaxF = ex1, jmaxF = ex2, kmaxF = ex3;

#if (ghost_width == 2)
    /* ---- 2nd-order ------------------------------------------------------ */
    {
        const int ord = 1;

        int iminF = 1, jminF = 1, kminF = 1;
        if (Symmetry > NO_SYMM && fabs(Z[0]) < dZ) kminF = 0;
        if (Symmetry > EQ_SYMM && fabs(X[0]) < dX) iminF = 0;
        if (Symmetry > EQ_SYMM && fabs(Y[0]) < dY) jminF = 0;

        const double SoA[3] = { SYM1, SYM2, SYM3 };

        const size_t nx = (size_t)ex1 + ord;
        const size_t ny = (size_t)ex2 + ord;
        const size_t nz = (size_t)ex3 + ord;
        const size_t fh_size = nx * ny * nz;

        static double *fh_buf = NULL;
        static size_t cap = 0;
        if (fh_size > cap) {
            free(fh_buf);
            fh_buf = (double*)aligned_alloc(64, fh_size * sizeof(double));
            cap = fh_size;
        }
        double *fh = fh_buf;
        if (!fh) return;

        symmetry_bd(ord, ex, f, fh, SoA);

        const double Sdxdx = ONE / (dX * dX);
        const double Sdydy = ONE / (dY * dY);
        const double Sdzdz = ONE / (dZ * dZ);
        const double Sdxdy = F1o4 / (dX * dY);
        const double Sdxdz = F1o4 / (dX * dZ);
        const double Sdydz = F1o4 / (dY * dZ);

        const size_t all = (size_t)ex1 * (size_t)ex2 * (size_t)ex3;
        for (size_t p = 0; p < all; ++p) {
            fxx[p]=ZEO; fyy[p]=ZEO; fzz[p]=ZEO;
            fxy[p]=ZEO; fxz[p]=ZEO; fyz[p]=ZEO;
        }

        const int i2_lo = (iminF > 0) ? iminF : 0;
        const int j2_lo = (jminF > 0) ? jminF : 0;
        const int k2_lo = (kminF > 0) ? kminF : 0;
        const int i2_hi = ex1 - 2;
        const int j2_hi = ex2 - 2;
        const int k2_hi = ex3 - 2;

        if (i2_lo <= i2_hi && j2_lo <= j2_hi && k2_lo <= k2_hi) {
            for (int k0 = k2_lo; k0 <= k2_hi; ++k0) {
                const int kF = k0;
                for (int j0 = j2_lo; j0 <= j2_hi; ++j0) {
                    const int jF = j0;
                    for (int i0 = i2_lo; i0 <= i2_hi; ++i0) {
                        const int iF = i0;
                        const size_t p = idx_ex(i0, j0, k0, ex);

                        fxx[p] = Sdxdx * (
                            fh[idx_fh_F_ord1(iF - 1, jF,     kF,     ex)] -
                            TWO * fh[idx_fh_F_ord1(iF,     jF,     kF,     ex)] +
                            fh[idx_fh_F_ord1(iF + 1, jF,     kF,     ex)]);

                        fyy[p] = Sdydy * (
                            fh[idx_fh_F_ord1(iF,     jF - 1, kF,     ex)] -
                            TWO * fh[idx_fh_F_ord1(iF,     jF,     kF,     ex)] +
                            fh[idx_fh_F_ord1(iF,     jF + 1, kF,     ex)]);

                        fzz[p] = Sdzdz * (
                            fh[idx_fh_F_ord1(iF,     jF,     kF - 1, ex)] -
                            TWO * fh[idx_fh_F_ord1(iF,     jF,     kF,     ex)] +
                            fh[idx_fh_F_ord1(iF,     jF,     kF + 1, ex)]);

                        fxy[p] = Sdxdy * (
                            fh[idx_fh_F_ord1(iF - 1, jF - 1, kF, ex)] -
                            fh[idx_fh_F_ord1(iF + 1, jF - 1, kF, ex)] -
                            fh[idx_fh_F_ord1(iF - 1, jF + 1, kF, ex)] +
                            fh[idx_fh_F_ord1(iF + 1, jF + 1, kF, ex)]);

                        fxz[p] = Sdxdz * (
                            fh[idx_fh_F_ord1(iF - 1, jF, kF - 1, ex)] -
                            fh[idx_fh_F_ord1(iF + 1, jF, kF - 1, ex)] -
                            fh[idx_fh_F_ord1(iF - 1, jF, kF + 1, ex)] +
                            fh[idx_fh_F_ord1(iF + 1, jF, kF + 1, ex)]);

                        fyz[p] = Sdydz * (
                            fh[idx_fh_F_ord1(iF, jF - 1, kF - 1, ex)] -
                            fh[idx_fh_F_ord1(iF, jF + 1, kF - 1, ex)] -
                            fh[idx_fh_F_ord1(iF, jF - 1, kF + 1, ex)] +
                            fh[idx_fh_F_ord1(iF, jF + 1, kF + 1, ex)]);
                    }
                }
            }
        }
        return;
    }
#elif (ghost_width == 3)
    /* ---- 4th-order (original code) ------------------------------------ */
    {
        const int ord = 2;

        int iminF = 1, jminF = 1, kminF = 1;
        if (Symmetry > NO_SYMM && fabs(Z[0]) < dZ) kminF = -1;
        if (Symmetry > EQ_SYMM && fabs(X[0]) < dX) iminF = -1;
        if (Symmetry > EQ_SYMM && fabs(Y[0]) < dY) jminF = -1;

        const double SoA[3] = { SYM1, SYM2, SYM3 };

        const size_t nx = (size_t)ex1 + ord;
        const size_t ny = (size_t)ex2 + ord;
        const size_t nz = (size_t)ex3 + ord;
        const size_t fh_size = nx * ny * nz;

        static double *fh_buf = NULL;
        static size_t cap = 0;
        if (fh_size > cap) {
            free(fh_buf);
            fh_buf = (double*)aligned_alloc(64, fh_size * sizeof(double));
            cap = fh_size;
        }
        double *fh = fh_buf;
        if (!fh) return;

        symmetry_bd(ord, ex, f, fh, SoA);

        const double Sdxdx = ONE / (dX * dX);
        const double Sdydy = ONE / (dY * dY);
        const double Sdzdz = ONE / (dZ * dZ);
        const double Fdxdx = F1o12 / (dX * dX);
        const double Fdydy = F1o12 / (dY * dY);
        const double Fdzdz = F1o12 / (dZ * dZ);
        const double Sdxdy = F1o4 / (dX * dY);
        const double Sdxdz = F1o4 / (dX * dZ);
        const double Sdydz = F1o4 / (dY * dZ);
        const double Fdxdy = F1o144 / (dX * dY);
        const double Fdxdz = F1o144 / (dX * dZ);
        const double Fdydz = F1o144 / (dY * dZ);

        /* zero high-boundary faces (points the loops below won't cover) */
        for (int j0 = 0; j0 < ex2; ++j0)
            for (int i0 = 0; i0 < ex1; ++i0) {
                const size_t p = idx_ex(i0, j0, ex3 - 1, ex);
                fxx[p]=ZEO; fyy[p]=ZEO; fzz[p]=ZEO;
                fxy[p]=ZEO; fxz[p]=ZEO; fyz[p]=ZEO;
            }
        for (int k0 = 0; k0 < ex3 - 1; ++k0)
            for (int i0 = 0; i0 < ex1; ++i0) {
                const size_t p = idx_ex(i0, ex2 - 1, k0, ex);
                fxx[p]=ZEO; fyy[p]=ZEO; fzz[p]=ZEO;
                fxy[p]=ZEO; fxz[p]=ZEO; fyz[p]=ZEO;
            }
        for (int k0 = 0; k0 < ex3 - 1; ++k0)
            for (int j0 = 0; j0 < ex2 - 1; ++j0) {
                const size_t p = idx_ex(ex1 - 1, j0, k0, ex);
                fxx[p]=ZEO; fyy[p]=ZEO; fzz[p]=ZEO;
                fxy[p]=ZEO; fxz[p]=ZEO; fyz[p]=ZEO;
            }

        if (kminF == 1) {
            for (int j0 = 0; j0 < ex2; ++j0)
                for (int i0 = 0; i0 < ex1; ++i0) {
                    const size_t p = idx_ex(i0, j0, 0, ex);
                    fxx[p]=ZEO; fyy[p]=ZEO; fzz[p]=ZEO;
                    fxy[p]=ZEO; fxz[p]=ZEO; fyz[p]=ZEO;
                }
        }
        if (jminF == 1) {
            for (int k0 = 0; k0 < ex3; ++k0)
                for (int i0 = 0; i0 < ex1; ++i0) {
                    const size_t p = idx_ex(i0, 0, k0, ex);
                    fxx[p]=ZEO; fyy[p]=ZEO; fzz[p]=ZEO;
                    fxy[p]=ZEO; fxz[p]=ZEO; fyz[p]=ZEO;
                }
        }
        if (iminF == 1) {
            for (int k0 = 0; k0 < ex3; ++k0)
                for (int j0 = 0; j0 < ex2; ++j0) {
                    const size_t p = idx_ex(0, j0, k0, ex);
                    fxx[p]=ZEO; fyy[p]=ZEO; fzz[p]=ZEO;
                    fxy[p]=ZEO; fxz[p]=ZEO; fyz[p]=ZEO;
                }
        }

        const int i2_lo = (iminF > 0) ? iminF : 0;
        const int j2_lo = (jminF > 0) ? jminF : 0;
        const int k2_lo = (kminF > 0) ? kminF : 0;
        const int i2_hi = ex1 - 2;
        const int j2_hi = ex2 - 2;
        const int k2_hi = ex3 - 2;

        const int i4_lo = (iminF + 1 > 0) ? (iminF + 1) : 0;
        const int j4_lo = (jminF + 1 > 0) ? (jminF + 1) : 0;
        const int k4_lo = (kminF + 1 > 0) ? (kminF + 1) : 0;
        const int i4_hi = ex1 - 3;
        const int j4_hi = ex2 - 3;
        const int k4_hi = ex3 - 3;

        const int has4 = (i4_lo <= i4_hi && j4_lo <= j4_hi && k4_lo <= k4_hi);

        if (i2_lo <= i2_hi && j2_lo <= j2_hi && k2_lo <= k2_hi) {
            for (int k0 = k2_lo; k0 <= k2_hi; ++k0) {
                const int kF = k0 + 1;
                for (int j0 = j2_lo; j0 <= j2_hi; ++j0) {
                    const int jF = j0 + 1;
                    for (int i0 = i2_lo; i0 <= i2_hi; ++i0) {
                        if (has4 &&
                            i0 >= i4_lo && i0 <= i4_hi &&
                            j0 >= j4_lo && j0 <= j4_hi &&
                            k0 >= k4_lo && k0 <= k4_hi) {
                            continue;
                        }
                        const int iF = i0 + 1;
                        const size_t p = idx_ex(i0, j0, k0, ex);

                        fxx[p] = Sdxdx * (
                            fh[idx_fh_F_ord2(iF - 1, jF,     kF,     ex)] -
                            TWO * fh[idx_fh_F_ord2(iF,     jF,     kF,     ex)] +
                            fh[idx_fh_F_ord2(iF + 1, jF,     kF,     ex)]);

                        fyy[p] = Sdydy * (
                            fh[idx_fh_F_ord2(iF,     jF - 1, kF,     ex)] -
                            TWO * fh[idx_fh_F_ord2(iF,     jF,     kF,     ex)] +
                            fh[idx_fh_F_ord2(iF,     jF + 1, kF,     ex)]);

                        fzz[p] = Sdzdz * (
                            fh[idx_fh_F_ord2(iF,     jF,     kF - 1, ex)] -
                            TWO * fh[idx_fh_F_ord2(iF,     jF,     kF,     ex)] +
                            fh[idx_fh_F_ord2(iF,     jF,     kF + 1, ex)]);

                        fxy[p] = Sdxdy * (
                            fh[idx_fh_F_ord2(iF - 1, jF - 1, kF, ex)] -
                            fh[idx_fh_F_ord2(iF + 1, jF - 1, kF, ex)] -
                            fh[idx_fh_F_ord2(iF - 1, jF + 1, kF, ex)] +
                            fh[idx_fh_F_ord2(iF + 1, jF + 1, kF, ex)]);

                        fxz[p] = Sdxdz * (
                            fh[idx_fh_F_ord2(iF - 1, jF, kF - 1, ex)] -
                            fh[idx_fh_F_ord2(iF + 1, jF, kF - 1, ex)] -
                            fh[idx_fh_F_ord2(iF - 1, jF, kF + 1, ex)] +
                            fh[idx_fh_F_ord2(iF + 1, jF, kF + 1, ex)]);

                        fyz[p] = Sdydz * (
                            fh[idx_fh_F_ord2(iF, jF - 1, kF - 1, ex)] -
                            fh[idx_fh_F_ord2(iF, jF + 1, kF - 1, ex)] -
                            fh[idx_fh_F_ord2(iF, jF - 1, kF + 1, ex)] +
                            fh[idx_fh_F_ord2(iF, jF + 1, kF + 1, ex)]);
                    }
                }
            }
        }

        if (has4) {
            for (int k0 = k4_lo; k0 <= k4_hi; ++k0) {
                const int kF = k0 + 1;
                for (int j0 = j4_lo; j0 <= j4_hi; ++j0) {
                    const int jF = j0 + 1;
                    for (int i0 = i4_lo; i0 <= i4_hi; ++i0) {
                        const int iF = i0 + 1;
                        const size_t p = idx_ex(i0, j0, k0, ex);

                        fxx[p] = Fdxdx * (
                            -fh[idx_fh_F_ord2(iF - 2, jF,     kF,     ex)] +
                            F16 * fh[idx_fh_F_ord2(iF - 1, jF,     kF,     ex)] -
                            F30 * fh[idx_fh_F_ord2(iF,     jF,     kF,     ex)] -
                            fh[idx_fh_F_ord2(iF + 2, jF,     kF,     ex)] +
                            F16 * fh[idx_fh_F_ord2(iF + 1, jF,     kF,     ex)]);

                        fyy[p] = Fdydy * (
                            -fh[idx_fh_F_ord2(iF,     jF - 2, kF,     ex)] +
                            F16 * fh[idx_fh_F_ord2(iF,     jF - 1, kF,     ex)] -
                            F30 * fh[idx_fh_F_ord2(iF,     jF,     kF,     ex)] -
                            fh[idx_fh_F_ord2(iF,     jF + 2, kF,     ex)] +
                            F16 * fh[idx_fh_F_ord2(iF,     jF + 1, kF,     ex)]);

                        fzz[p] = Fdzdz * (
                            -fh[idx_fh_F_ord2(iF,     jF,     kF - 2, ex)] +
                            F16 * fh[idx_fh_F_ord2(iF,     jF,     kF - 1, ex)] -
                            F30 * fh[idx_fh_F_ord2(iF,     jF,     kF,     ex)] -
                            fh[idx_fh_F_ord2(iF,     jF,     kF + 2, ex)] +
                            F16 * fh[idx_fh_F_ord2(iF,     jF,     kF + 1, ex)]);

                        /* fxy: 5x5 outer product */
                        {
                            const double t_jm2 = (
                                fh[idx_fh_F_ord2(iF - 2, jF - 2, kF, ex)]
                                -F8*fh[idx_fh_F_ord2(iF - 1, jF - 2, kF, ex)]
                                +F8*fh[idx_fh_F_ord2(iF + 1, jF - 2, kF, ex)]
                                -    fh[idx_fh_F_ord2(iF + 2, jF - 2, kF, ex)]);

                            const double t_jm1 = (
                                fh[idx_fh_F_ord2(iF - 2, jF - 1, kF, ex)]
                                -F8*fh[idx_fh_F_ord2(iF - 1, jF - 1, kF, ex)]
                                +F8*fh[idx_fh_F_ord2(iF + 1, jF - 1, kF, ex)]
                                -    fh[idx_fh_F_ord2(iF + 2, jF - 1, kF, ex)]);

                            const double t_jp1 = (
                                fh[idx_fh_F_ord2(iF - 2, jF + 1, kF, ex)]
                                -F8*fh[idx_fh_F_ord2(iF - 1, jF + 1, kF, ex)]
                                +F8*fh[idx_fh_F_ord2(iF + 1, jF + 1, kF, ex)]
                                -    fh[idx_fh_F_ord2(iF + 2, jF + 1, kF, ex)]);

                            const double t_jp2 = (
                                fh[idx_fh_F_ord2(iF - 2, jF + 2, kF, ex)]
                                -F8*fh[idx_fh_F_ord2(iF - 1, jF + 2, kF, ex)]
                                +F8*fh[idx_fh_F_ord2(iF + 1, jF + 2, kF, ex)]
                                -    fh[idx_fh_F_ord2(iF + 2, jF + 2, kF, ex)]);

                            fxy[p] = Fdxdy * ( t_jm2 - F8 * t_jm1 + F8 * t_jp1 - t_jp2 );
                        }

                        /* fxz: 5x5 outer product */
                        {
                            const double t_km2 = (
                                fh[idx_fh_F_ord2(iF - 2, jF, kF - 2, ex)]
                                -F8*fh[idx_fh_F_ord2(iF - 1, jF, kF - 2, ex)]
                                +F8*fh[idx_fh_F_ord2(iF + 1, jF, kF - 2, ex)]
                                -    fh[idx_fh_F_ord2(iF + 2, jF, kF - 2, ex)]);

                            const double t_km1 = (
                                fh[idx_fh_F_ord2(iF - 2, jF, kF - 1, ex)]
                                -F8*fh[idx_fh_F_ord2(iF - 1, jF, kF - 1, ex)]
                                +F8*fh[idx_fh_F_ord2(iF + 1, jF, kF - 1, ex)]
                                -    fh[idx_fh_F_ord2(iF + 2, jF, kF - 1, ex)]);

                            const double t_kp1 = (
                                fh[idx_fh_F_ord2(iF - 2, jF, kF + 1, ex)]
                                -F8*fh[idx_fh_F_ord2(iF - 1, jF, kF + 1, ex)]
                                +F8*fh[idx_fh_F_ord2(iF + 1, jF, kF + 1, ex)]
                                -    fh[idx_fh_F_ord2(iF + 2, jF, kF + 1, ex)]);

                            const double t_kp2 = (
                                fh[idx_fh_F_ord2(iF - 2, jF, kF + 2, ex)]
                                -F8*fh[idx_fh_F_ord2(iF - 1, jF, kF + 2, ex)]
                                +F8*fh[idx_fh_F_ord2(iF + 1, jF, kF + 2, ex)]
                                -    fh[idx_fh_F_ord2(iF + 2, jF, kF + 2, ex)]);

                            fxz[p] = Fdxdz * ( t_km2 - F8 * t_km1 + F8 * t_kp1 - t_kp2 );
                        }

                        /* fyz: 5x5 outer product */
                        {
                            const double t_km2 = (
                                fh[idx_fh_F_ord2(iF, jF - 2, kF - 2, ex)]
                                -F8*fh[idx_fh_F_ord2(iF, jF - 1, kF - 2, ex)]
                                +F8*fh[idx_fh_F_ord2(iF, jF + 1, kF - 2, ex)]
                                -    fh[idx_fh_F_ord2(iF, jF + 2, kF - 2, ex)]);

                            const double t_km1 = (
                                fh[idx_fh_F_ord2(iF, jF - 2, kF - 1, ex)]
                                -F8*fh[idx_fh_F_ord2(iF, jF - 1, kF - 1, ex)]
                                +F8*fh[idx_fh_F_ord2(iF, jF + 1, kF - 1, ex)]
                                -    fh[idx_fh_F_ord2(iF, jF + 2, kF - 1, ex)]);

                            const double t_kp1 = (
                                fh[idx_fh_F_ord2(iF, jF - 2, kF + 1, ex)]
                                -F8*fh[idx_fh_F_ord2(iF, jF - 1, kF + 1, ex)]
                                +F8*fh[idx_fh_F_ord2(iF, jF + 1, kF + 1, ex)]
                                -    fh[idx_fh_F_ord2(iF, jF + 2, kF + 1, ex)]);

                            const double t_kp2 = (
                                fh[idx_fh_F_ord2(iF, jF - 2, kF + 2, ex)]
                                -F8*fh[idx_fh_F_ord2(iF, jF - 1, kF + 2, ex)]
                                +F8*fh[idx_fh_F_ord2(iF, jF + 1, kF + 2, ex)]
                                -    fh[idx_fh_F_ord2(iF, jF + 2, kF + 2, ex)]);

                            fyz[p] = Fdydz * ( t_km2 - F8 * t_km1 + F8 * t_kp1 - t_kp2 );
                        }
                    }
                }
            }
        }
        return;
    }
#elif (ghost_width == 4)
    /* ---- 6th-order ----------------------------------------------------- */
    {
        const int ord = 3;

        int iminF = 1, jminF = 1, kminF = 1;
        if (Symmetry > NO_SYMM && fabs(Z[0]) < dZ) kminF = -2;
        if (Symmetry > EQ_SYMM && fabs(X[0]) < dX) iminF = -2;
        if (Symmetry > EQ_SYMM && fabs(Y[0]) < dY) jminF = -2;

        const double SoA[3] = { SYM1, SYM2, SYM3 };

        const size_t nx = (size_t)ex1 + ord;
        const size_t ny = (size_t)ex2 + ord;
        const size_t nz = (size_t)ex3 + ord;
        const size_t fh_size = nx * ny * nz;
        static double *fh_buf = NULL;
        static size_t cap = 0;
        if (fh_size > cap) {
            free(fh_buf);
            fh_buf = (double*)aligned_alloc(64, fh_size * sizeof(double));
            cap = fh_size;
        }
        double *fh = fh_buf;
        if (!fh) return;

        symmetry_bd(ord, ex, f, fh, SoA);

        /* Denominators */
        const double Sdxdx = ONE / (dX * dX);           // 2nd
        const double Sdydy = ONE / (dY * dY);
        const double Sdzdz = ONE / (dZ * dZ);
        const double Fdxdx = F1o12 / (dX * dX);          // 4th
        const double Fdydy = F1o12 / (dY * dY);
        const double Fdzdz = F1o12 / (dZ * dZ);
        const double Xdxdx = F1o180 / (dX * dX);         // 6th
        const double Xdydy = F1o180 / (dY * dY);
        const double Xdzdz = F1o180 / (dZ * dZ);
        const double Sdxdy = F1o4    / (dX * dY);
        const double Sdxdz = F1o4    / (dX * dZ);
        const double Sdydz = F1o4    / (dY * dZ);
        const double Fdxdy = F1o144  / (dX * dY);
        const double Fdxdz = F1o144  / (dX * dZ);
        const double Fdydz = F1o144  / (dY * dZ);
        const double Xdxdy = F1o3600 / (dX * dY);
        const double Xdxdz = F1o3600 / (dX * dZ);
        const double Xdydz = F1o3600 / (dY * dZ);

        /* zero everything first */
        const size_t all = (size_t)ex1 * (size_t)ex2 * (size_t)ex3;
        for (size_t p = 0; p < all; ++p) {
            fxx[p]=ZEO; fyy[p]=ZEO; fzz[p]=ZEO;
            fxy[p]=ZEO; fxz[p]=ZEO; fyz[p]=ZEO;
        }

        /* loop bounds for each pass (from widest to narrowest) */
        const int i2_lo = (iminF + 1 > 0) ? (iminF + 1) : 0;
        const int j2_lo = (jminF + 1 > 0) ? (jminF + 1) : 0;
        const int k2_lo = (kminF + 1 > 0) ? (kminF + 1) : 0;
        const int i2_hi = ex1 - 2, j2_hi = ex2 - 2, k2_hi = ex3 - 2;

        const int i4_lo = (iminF + 2 > 0) ? (iminF + 2) : 0;
        const int j4_lo = (jminF + 2 > 0) ? (jminF + 2) : 0;
        const int k4_lo = (kminF + 2 > 0) ? (kminF + 2) : 0;
        const int i4_hi = ex1 - 3, j4_hi = ex2 - 3, k4_hi = ex3 - 3;

        const int i6_lo = (iminF + 3 > 0) ? (iminF + 3) : 0;
        const int j6_lo = (jminF + 3 > 0) ? (jminF + 3) : 0;
        const int k6_lo = (kminF + 3 > 0) ? (kminF + 3) : 0;
        const int i6_hi = ex1 - 4, j6_hi = ex2 - 4, k6_hi = ex3 - 4;

        const int has4 = (i4_lo <= i4_hi && j4_lo <= j4_hi && k4_lo <= k4_hi);
        const int has6 = (i6_lo <= i6_hi && j6_lo <= j6_hi && k6_lo <= k6_hi);

        /* 2nd-order: skip 4th+6th overlap */
        if (i2_lo <= i2_hi && j2_lo <= j2_hi && k2_lo <= k2_hi) {
            for (int k0 = k2_lo; k0 <= k2_hi; ++k0) {
                const int kF = k0 + 2;
                for (int j0 = j2_lo; j0 <= j2_hi; ++j0) {
                    const int jF = j0 + 2;
                    for (int i0 = i2_lo; i0 <= i2_hi; ++i0) {
                        _Bool in4 = has4 && i0>=i4_lo && i0<=i4_hi && j0>=j4_lo && j0<=j4_hi && k0>=k4_lo && k0<=k4_hi;
                        if (in4) continue;
                        const int iF = i0 + 2;
                        const size_t p = idx_ex(i0, j0, k0, ex);

                        fxx[p] = Sdxdx * (fh[idx_fh_F(iF - 1, jF,     kF,     ex)] - TWO*fh[idx_fh_F(iF,jF,kF,ex)] + fh[idx_fh_F(iF + 1, jF,     kF,     ex)]);
                        fyy[p] = Sdydy * (fh[idx_fh_F(iF,     jF - 1, kF,     ex)] - TWO*fh[idx_fh_F(iF,jF,kF,ex)] + fh[idx_fh_F(iF,     jF + 1, kF,     ex)]);
                        fzz[p] = Sdzdz * (fh[idx_fh_F(iF,     jF,     kF - 1, ex)] - TWO*fh[idx_fh_F(iF,jF,kF,ex)] + fh[idx_fh_F(iF,     jF,     kF + 1, ex)]);

                        fxy[p] = Sdxdy * (fh[idx_fh_F(iF - 1, jF - 1, kF, ex)] - fh[idx_fh_F(iF + 1, jF - 1, kF, ex)] - fh[idx_fh_F(iF - 1, jF + 1, kF, ex)] + fh[idx_fh_F(iF + 1, jF + 1, kF, ex)]);
                        fxz[p] = Sdxdz * (fh[idx_fh_F(iF - 1, jF, kF - 1, ex)] - fh[idx_fh_F(iF + 1, jF, kF - 1, ex)] - fh[idx_fh_F(iF - 1, jF, kF + 1, ex)] + fh[idx_fh_F(iF + 1, jF, kF + 1, ex)]);
                        fyz[p] = Sdydz * (fh[idx_fh_F(iF, jF - 1, kF - 1, ex)] - fh[idx_fh_F(iF, jF + 1, kF - 1, ex)] - fh[idx_fh_F(iF, jF - 1, kF + 1, ex)] + fh[idx_fh_F(iF, jF + 1, kF + 1, ex)]);
                    }
                }
            }
        }

        /* 4th-order: skip 6th overlap */
        if (has4) {
            for (int k0 = k4_lo; k0 <= k4_hi; ++k0) {
                const int kF = k0 + 2;
                for (int j0 = j4_lo; j0 <= j4_hi; ++j0) {
                    const int jF = j0 + 2;
                    for (int i0 = i4_lo; i0 <= i4_hi; ++i0) {
                        if (has6 && i0>=i6_lo && i0<=i6_hi && j0>=j6_lo && j0<=j6_hi && k0>=k6_lo && k0<=k6_hi) continue;
                        const int iF = i0 + 2;
                        const size_t p = idx_ex(i0, j0, k0, ex);

                        fxx[p] = Fdxdx * (-fh[idx_fh_F(iF - 2, jF,     kF,     ex)] + F16*fh[idx_fh_F(iF-1,jF,kF,ex)] - F30*fh[idx_fh_F(iF,jF,kF,ex)] - fh[idx_fh_F(iF+2,jF,kF,ex)] + F16*fh[idx_fh_F(iF+1,jF,kF,ex)]);
                        fyy[p] = Fdydy * (-fh[idx_fh_F(iF,     jF - 2, kF,     ex)] + F16*fh[idx_fh_F(iF,jF-1,kF,ex)] - F30*fh[idx_fh_F(iF,jF,kF,ex)] - fh[idx_fh_F(iF,jF+2,kF,ex)] + F16*fh[idx_fh_F(iF,jF+1,kF,ex)]);
                        fzz[p] = Fdzdz * (-fh[idx_fh_F(iF,     jF,     kF - 2, ex)] + F16*fh[idx_fh_F(iF,jF,kF-1,ex)] - F30*fh[idx_fh_F(iF,jF,kF,ex)] - fh[idx_fh_F(iF,jF,kF+2,ex)] + F16*fh[idx_fh_F(iF,jF,kF+1,ex)]);

                        {
                            const double t_jm2 = (fh[idx_fh_F(iF-2,jF-2,kF,ex)]-F8*fh[idx_fh_F(iF-1,jF-2,kF,ex)]+F8*fh[idx_fh_F(iF+1,jF-2,kF,ex)]-fh[idx_fh_F(iF+2,jF-2,kF,ex)]);
                            const double t_jm1 = (fh[idx_fh_F(iF-2,jF-1,kF,ex)]-F8*fh[idx_fh_F(iF-1,jF-1,kF,ex)]+F8*fh[idx_fh_F(iF+1,jF-1,kF,ex)]-fh[idx_fh_F(iF+2,jF-1,kF,ex)]);
                            const double t_jp1 = (fh[idx_fh_F(iF-2,jF+1,kF,ex)]-F8*fh[idx_fh_F(iF-1,jF+1,kF,ex)]+F8*fh[idx_fh_F(iF+1,jF+1,kF,ex)]-fh[idx_fh_F(iF+2,jF+1,kF,ex)]);
                            const double t_jp2 = (fh[idx_fh_F(iF-2,jF+2,kF,ex)]-F8*fh[idx_fh_F(iF-1,jF+2,kF,ex)]+F8*fh[idx_fh_F(iF+1,jF+2,kF,ex)]-fh[idx_fh_F(iF+2,jF+2,kF,ex)]);
                            fxy[p] = Fdxdy * (t_jm2 - F8*t_jm1 + F8*t_jp1 - t_jp2);
                        }
                        {
                            const double t_km2 = (fh[idx_fh_F(iF-2,jF,kF-2,ex)]-F8*fh[idx_fh_F(iF-1,jF,kF-2,ex)]+F8*fh[idx_fh_F(iF+1,jF,kF-2,ex)]-fh[idx_fh_F(iF+2,jF,kF-2,ex)]);
                            const double t_km1 = (fh[idx_fh_F(iF-2,jF,kF-1,ex)]-F8*fh[idx_fh_F(iF-1,jF,kF-1,ex)]+F8*fh[idx_fh_F(iF+1,jF,kF-1,ex)]-fh[idx_fh_F(iF+2,jF,kF-1,ex)]);
                            const double t_kp1 = (fh[idx_fh_F(iF-2,jF,kF+1,ex)]-F8*fh[idx_fh_F(iF-1,jF,kF+1,ex)]+F8*fh[idx_fh_F(iF+1,jF,kF+1,ex)]-fh[idx_fh_F(iF+2,jF,kF+1,ex)]);
                            const double t_kp2 = (fh[idx_fh_F(iF-2,jF,kF+2,ex)]-F8*fh[idx_fh_F(iF-1,jF,kF+2,ex)]+F8*fh[idx_fh_F(iF+1,jF,kF+2,ex)]-fh[idx_fh_F(iF+2,jF,kF+2,ex)]);
                            fxz[p] = Fdxdz * (t_km2 - F8*t_km1 + F8*t_kp1 - t_kp2);
                        }
                        {
                            const double t_km2 = (fh[idx_fh_F(iF,jF-2,kF-2,ex)]-F8*fh[idx_fh_F(iF,jF-1,kF-2,ex)]+F8*fh[idx_fh_F(iF,jF+1,kF-2,ex)]-fh[idx_fh_F(iF,jF+2,kF-2,ex)]);
                            const double t_km1 = (fh[idx_fh_F(iF,jF-2,kF-1,ex)]-F8*fh[idx_fh_F(iF,jF-1,kF-1,ex)]+F8*fh[idx_fh_F(iF,jF+1,kF-1,ex)]-fh[idx_fh_F(iF,jF+2,kF-1,ex)]);
                            const double t_kp1 = (fh[idx_fh_F(iF,jF-2,kF+1,ex)]-F8*fh[idx_fh_F(iF,jF-1,kF+1,ex)]+F8*fh[idx_fh_F(iF,jF+1,kF+1,ex)]-fh[idx_fh_F(iF,jF+2,kF+1,ex)]);
                            const double t_kp2 = (fh[idx_fh_F(iF,jF-2,kF+2,ex)]-F8*fh[idx_fh_F(iF,jF-1,kF+2,ex)]+F8*fh[idx_fh_F(iF,jF+1,kF+2,ex)]-fh[idx_fh_F(iF,jF+2,kF+2,ex)]);
                            fyz[p] = Fdydz * (t_km2 - F8*t_km1 + F8*t_kp1 - t_kp2);
                        }
                    }
                }
            }
        }

        /* 6th-order: interior only */
        if (has6) {
            for (int k0 = k6_lo; k0 <= k6_hi; ++k0) {
                const int kF = k0 + 2;
                for (int j0 = j6_lo; j0 <= j6_hi; ++j0) {
                    const int jF = j0 + 2;
                    for (int i0 = i6_lo; i0 <= i6_hi; ++i0) {
                        const int iF = i0 + 2;
                        const size_t p = idx_ex(i0, j0, k0, ex);

                        /* Diagonal: [+2,-27,+270,-490,+270,-27,+2] / (180*dx^2) */
                        fxx[p] = Xdxdx * (
                            TWO  * fh[idx_fh_F(iF - 3, jF,     kF,     ex)] -
                            F27  * fh[idx_fh_F(iF - 2, jF,     kF,     ex)] +
                            F270 * fh[idx_fh_F(iF - 1, jF,     kF,     ex)] -
                            F490 * fh[idx_fh_F(iF,     jF,     kF,     ex)] +
                            F270 * fh[idx_fh_F(iF + 1, jF,     kF,     ex)] -
                            F27  * fh[idx_fh_F(iF + 2, jF,     kF,     ex)] +
                            TWO  * fh[idx_fh_F(iF + 3, jF,     kF,     ex)]);

                        fyy[p] = Xdydy * (
                            TWO  * fh[idx_fh_F(iF,     jF - 3, kF,     ex)] -
                            F27  * fh[idx_fh_F(iF,     jF - 2, kF,     ex)] +
                            F270 * fh[idx_fh_F(iF,     jF - 1, kF,     ex)] -
                            F490 * fh[idx_fh_F(iF,     jF,     kF,     ex)] +
                            F270 * fh[idx_fh_F(iF,     jF + 1, kF,     ex)] -
                            F27  * fh[idx_fh_F(iF,     jF + 2, kF,     ex)] +
                            TWO  * fh[idx_fh_F(iF,     jF + 3, kF,     ex)]);

                        fzz[p] = Xdzdz * (
                            TWO  * fh[idx_fh_F(iF,     jF,     kF - 3, ex)] -
                            F27  * fh[idx_fh_F(iF,     jF,     kF - 2, ex)] +
                            F270 * fh[idx_fh_F(iF,     jF,     kF - 1, ex)] -
                            F490 * fh[idx_fh_F(iF,     jF,     kF,     ex)] +
                            F270 * fh[idx_fh_F(iF,     jF,     kF + 1, ex)] -
                            F27  * fh[idx_fh_F(iF,     jF,     kF + 2, ex)] +
                            TWO  * fh[idx_fh_F(iF,     jF,     kF + 3, ex)]);

                        /* Mixed: 7x7 outer product. Compute 1D x-stencil at each y/z offset,
                           then combine using 1D y/z weights [-1,+9,-45,0,+45,-9,+1] / (3600*dx*dy) */
                        {
                            // x-stencil: -f(i-3)+9f(i-2)-45f(i-1)+45f(i+1)-9f(i+2)+f(i+3)
                            // Helper macro would help but explicit is safer
                            #define XSTEN6(JF, KF_DUMMY) \
                                (-fh[idx_fh_F(iF-3,JF,KF_DUMMY,ex)] + F9*fh[idx_fh_F(iF-2,JF,KF_DUMMY,ex)] - F45*fh[idx_fh_F(iF-1,JF,KF_DUMMY,ex)] + F45*fh[idx_fh_F(iF+1,JF,KF_DUMMY,ex)] - F9*fh[idx_fh_F(iF+2,JF,KF_DUMMY,ex)] + fh[idx_fh_F(iF+3,JF,KF_DUMMY,ex)])

                            fxy[p] = Xdxdy * (
                                -XSTEN6(jF-3, kF) + F9*XSTEN6(jF-2, kF) - F45*XSTEN6(jF-1, kF) + F45*XSTEN6(jF+1, kF) - F9*XSTEN6(jF+2, kF) + XSTEN6(jF+3, kF));

                            fxz[p] = Xdxdz * (
                                -XSTEN6(jF, kF-3) + F9*XSTEN6(jF, kF-2) - F45*XSTEN6(jF, kF-1) + F45*XSTEN6(jF, kF+1) - F9*XSTEN6(jF, kF+2) + XSTEN6(jF, kF+3));

                            #undef XSTEN6
                        }

                        /* fyz: apply 1D y-stencil at each z offset */
                        {
                            #define YSTEN6(JF, KF_DUMMY) \
                                (-fh[idx_fh_F(iF,JF-3,KF_DUMMY,ex)] + F9*fh[idx_fh_F(iF,JF-2,KF_DUMMY,ex)] - F45*fh[idx_fh_F(iF,JF-1,KF_DUMMY,ex)] + F45*fh[idx_fh_F(iF,JF+1,KF_DUMMY,ex)] - F9*fh[idx_fh_F(iF,JF+2,KF_DUMMY,ex)] + fh[idx_fh_F(iF,JF+3,KF_DUMMY,ex)])

                            fyz[p] = Xdydz * (
                                -YSTEN6(jF, kF-3) + F9*YSTEN6(jF, kF-2) - F45*YSTEN6(jF, kF-1) + F45*YSTEN6(jF, kF+1) - F9*YSTEN6(jF, kF+2) + YSTEN6(jF, kF+3));

                            #undef YSTEN6
                        }
                    }
                }
            }
        }
        return;
    }
#elif (ghost_width == 5)
    /* ---- 8th-order ----------------------------------------------------- */
    {
        const int ord = 4;

        int iminF = 1, jminF = 1, kminF = 1;
        if (Symmetry > NO_SYMM && fabs(Z[0]) < dZ) kminF = -3;
        if (Symmetry > EQ_SYMM && fabs(X[0]) < dX) iminF = -3;
        if (Symmetry > EQ_SYMM && fabs(Y[0]) < dY) jminF = -3;

        const double SoA[3] = { SYM1, SYM2, SYM3 };

        const size_t nx = (size_t)ex1 + ord;
        const size_t ny = (size_t)ex2 + ord;
        const size_t nz = (size_t)ex3 + ord;
        const size_t fh_size = nx * ny * nz;
        static double *fh_buf = NULL;
        static size_t cap = 0;
        if (fh_size > cap) {
            free(fh_buf);
            fh_buf = (double*)aligned_alloc(64, fh_size * sizeof(double));
            cap = fh_size;
        }
        double *fh = fh_buf;
        if (!fh) return;

        symmetry_bd(ord, ex, f, fh, SoA);

        const double Sdxdx = ONE / (dX * dX);
        const double Sdydy = ONE / (dY * dY);
        const double Sdzdz = ONE / (dZ * dZ);
        const double Fdxdx = F1o12  / (dX * dX);
        const double Fdydy = F1o12  / (dY * dY);
        const double Fdzdz = F1o12  / (dZ * dZ);
        const double Xdxdx = F1o180 / (dX * dX);
        const double Xdydy = F1o180 / (dY * dY);
        const double Xdzdz = F1o180 / (dZ * dZ);
        const double Edxdx = F1o5040 / (dX * dX);
        const double Edydy = F1o5040 / (dY * dY);
        const double Edzdz = F1o5040 / (dZ * dZ);
        const double Sdxdy = F1o4    / (dX * dY);
        const double Sdxdz = F1o4    / (dX * dZ);
        const double Sdydz = F1o4    / (dY * dZ);
        const double Fdxdy = F1o144  / (dX * dY);
        const double Fdxdz = F1o144  / (dX * dZ);
        const double Fdydz = F1o144  / (dY * dZ);
        const double Xdxdy = F1o3600 / (dX * dY);
        const double Xdxdz = F1o3600 / (dX * dZ);
        const double Xdydz = F1o3600 / (dY * dZ);
        const double Edxdy = F1o705600 / (dX * dY);
        const double Edxdz = F1o705600 / (dX * dZ);
        const double Edydz = F1o705600 / (dY * dZ);

        const size_t all = (size_t)ex1 * (size_t)ex2 * (size_t)ex3;
        for (size_t p = 0; p < all; ++p) {
            fxx[p]=ZEO; fyy[p]=ZEO; fzz[p]=ZEO;
            fxy[p]=ZEO; fxz[p]=ZEO; fyz[p]=ZEO;
        }

        /* Loop bounds for each pass */
        const int i2_lo = (iminF + 2 > 0) ? (iminF + 2) : 0;
        const int j2_lo = (jminF + 2 > 0) ? (jminF + 2) : 0;
        const int k2_lo = (kminF + 2 > 0) ? (kminF + 2) : 0;
        const int i2_hi = ex1 - 2, j2_hi = ex2 - 2, k2_hi = ex3 - 2;

        const int i4_lo = (iminF + 3 > 0) ? (iminF + 3) : 0;
        const int j4_lo = (jminF + 3 > 0) ? (jminF + 3) : 0;
        const int k4_lo = (kminF + 3 > 0) ? (kminF + 3) : 0;
        const int i4_hi = ex1 - 3, j4_hi = ex2 - 3, k4_hi = ex3 - 3;

        const int i6_lo = (iminF + 4 > 0) ? (iminF + 4) : 0;
        const int j6_lo = (jminF + 4 > 0) ? (jminF + 4) : 0;
        const int k6_lo = (kminF + 4 > 0) ? (kminF + 4) : 0;
        const int i6_hi = ex1 - 4, j6_hi = ex2 - 4, k6_hi = ex3 - 4;

        const int i8_lo = (iminF + 5 > 0) ? (iminF + 5) : 0;
        const int j8_lo = (jminF + 5 > 0) ? (jminF + 5) : 0;
        const int k8_lo = (kminF + 5 > 0) ? (kminF + 5) : 0;
        const int i8_hi = ex1 - 5, j8_hi = ex2 - 5, k8_hi = ex3 - 5;

        const int has4 = (i4_lo <= i4_hi && j4_lo <= j4_hi && k4_lo <= k4_hi);
        const int has6 = (i6_lo <= i6_hi && j6_lo <= j6_hi && k6_lo <= k6_hi);
        const int has8 = (i8_lo <= i8_hi && j8_lo <= j8_hi && k8_lo <= k8_hi);

        /* 2nd-order: skip 4th+6th+8th overlap */
        if (i2_lo <= i2_hi && j2_lo <= j2_hi && k2_lo <= k2_hi) {
            for (int k0 = k2_lo; k0 <= k2_hi; ++k0) {
                const int kF = k0 + 3;
                for (int j0 = j2_lo; j0 <= j2_hi; ++j0) {
                    const int jF = j0 + 3;
                    for (int i0 = i2_lo; i0 <= i2_hi; ++i0) {
                        _Bool in4 = has4 && i0>=i4_lo && i0<=i4_hi && j0>=j4_lo && j0<=j4_hi && k0>=k4_lo && k0<=k4_hi;
                        if (in4) continue;
                        const int iF = i0 + 3;
                        const size_t p = idx_ex(i0, j0, k0, ex);

                        fxx[p] = Sdxdx * (fh[idx_fh_F_ord4(iF-1,jF,kF,ex)] - TWO*fh[idx_fh_F_ord4(iF,jF,kF,ex)] + fh[idx_fh_F_ord4(iF+1,jF,kF,ex)]);
                        fyy[p] = Sdydy * (fh[idx_fh_F_ord4(iF,jF-1,kF,ex)] - TWO*fh[idx_fh_F_ord4(iF,jF,kF,ex)] + fh[idx_fh_F_ord4(iF,jF+1,kF,ex)]);
                        fzz[p] = Sdzdz * (fh[idx_fh_F_ord4(iF,jF,kF-1,ex)] - TWO*fh[idx_fh_F_ord4(iF,jF,kF,ex)] + fh[idx_fh_F_ord4(iF,jF,kF+1,ex)]);

                        fxy[p] = Sdxdy * (fh[idx_fh_F_ord4(iF-1,jF-1,kF,ex)] - fh[idx_fh_F_ord4(iF+1,jF-1,kF,ex)] - fh[idx_fh_F_ord4(iF-1,jF+1,kF,ex)] + fh[idx_fh_F_ord4(iF+1,jF+1,kF,ex)]);
                        fxz[p] = Sdxdz * (fh[idx_fh_F_ord4(iF-1,jF,kF-1,ex)] - fh[idx_fh_F_ord4(iF+1,jF,kF-1,ex)] - fh[idx_fh_F_ord4(iF-1,jF,kF+1,ex)] + fh[idx_fh_F_ord4(iF+1,jF,kF+1,ex)]);
                        fyz[p] = Sdydz * (fh[idx_fh_F_ord4(iF,jF-1,kF-1,ex)] - fh[idx_fh_F_ord4(iF,jF+1,kF-1,ex)] - fh[idx_fh_F_ord4(iF,jF-1,kF+1,ex)] + fh[idx_fh_F_ord4(iF,jF+1,kF+1,ex)]);
                    }
                }
            }
        }

        /* 4th-order: skip 6th+8th overlap */
        if (has4) {
            for (int k0 = k4_lo; k0 <= k4_hi; ++k0) {
                const int kF = k0 + 3;
                for (int j0 = j4_lo; j0 <= j4_hi; ++j0) {
                    const int jF = j0 + 3;
                    for (int i0 = i4_lo; i0 <= i4_hi; ++i0) {
                        _Bool in6 = has6 && i0>=i6_lo && i0<=i6_hi && j0>=j6_lo && j0<=j6_hi && k0>=k6_lo && k0<=k6_hi;
                        if (in6) continue;
                        const int iF = i0 + 3;
                        const size_t p = idx_ex(i0, j0, k0, ex);

                        fxx[p] = Fdxdx * (-fh[idx_fh_F_ord4(iF-2,jF,kF,ex)] + F16*fh[idx_fh_F_ord4(iF-1,jF,kF,ex)] - F30*fh[idx_fh_F_ord4(iF,jF,kF,ex)] - fh[idx_fh_F_ord4(iF+2,jF,kF,ex)] + F16*fh[idx_fh_F_ord4(iF+1,jF,kF,ex)]);
                        fyy[p] = Fdydy * (-fh[idx_fh_F_ord4(iF,jF-2,kF,ex)] + F16*fh[idx_fh_F_ord4(iF,jF-1,kF,ex)] - F30*fh[idx_fh_F_ord4(iF,jF,kF,ex)] - fh[idx_fh_F_ord4(iF,jF+2,kF,ex)] + F16*fh[idx_fh_F_ord4(iF,jF+1,kF,ex)]);
                        fzz[p] = Fdzdz * (-fh[idx_fh_F_ord4(iF,jF,kF-2,ex)] + F16*fh[idx_fh_F_ord4(iF,jF,kF-1,ex)] - F30*fh[idx_fh_F_ord4(iF,jF,kF,ex)] - fh[idx_fh_F_ord4(iF,jF,kF+2,ex)] + F16*fh[idx_fh_F_ord4(iF,jF,kF+1,ex)]);

                        {
                            const double t_jm2 = (fh[idx_fh_F_ord4(iF-2,jF-2,kF,ex)]-F8*fh[idx_fh_F_ord4(iF-1,jF-2,kF,ex)]+F8*fh[idx_fh_F_ord4(iF+1,jF-2,kF,ex)]-fh[idx_fh_F_ord4(iF+2,jF-2,kF,ex)]);
                            const double t_jm1 = (fh[idx_fh_F_ord4(iF-2,jF-1,kF,ex)]-F8*fh[idx_fh_F_ord4(iF-1,jF-1,kF,ex)]+F8*fh[idx_fh_F_ord4(iF+1,jF-1,kF,ex)]-fh[idx_fh_F_ord4(iF+2,jF-1,kF,ex)]);
                            const double t_jp1 = (fh[idx_fh_F_ord4(iF-2,jF+1,kF,ex)]-F8*fh[idx_fh_F_ord4(iF-1,jF+1,kF,ex)]+F8*fh[idx_fh_F_ord4(iF+1,jF+1,kF,ex)]-fh[idx_fh_F_ord4(iF+2,jF+1,kF,ex)]);
                            const double t_jp2 = (fh[idx_fh_F_ord4(iF-2,jF+2,kF,ex)]-F8*fh[idx_fh_F_ord4(iF-1,jF+2,kF,ex)]+F8*fh[idx_fh_F_ord4(iF+1,jF+2,kF,ex)]-fh[idx_fh_F_ord4(iF+2,jF+2,kF,ex)]);
                            fxy[p] = Fdxdy * (t_jm2 - F8*t_jm1 + F8*t_jp1 - t_jp2);
                        }
                        {
                            const double t_km2 = (fh[idx_fh_F_ord4(iF-2,jF,kF-2,ex)]-F8*fh[idx_fh_F_ord4(iF-1,jF,kF-2,ex)]+F8*fh[idx_fh_F_ord4(iF+1,jF,kF-2,ex)]-fh[idx_fh_F_ord4(iF+2,jF,kF-2,ex)]);
                            const double t_km1 = (fh[idx_fh_F_ord4(iF-2,jF,kF-1,ex)]-F8*fh[idx_fh_F_ord4(iF-1,jF,kF-1,ex)]+F8*fh[idx_fh_F_ord4(iF+1,jF,kF-1,ex)]-fh[idx_fh_F_ord4(iF+2,jF,kF-1,ex)]);
                            const double t_kp1 = (fh[idx_fh_F_ord4(iF-2,jF,kF+1,ex)]-F8*fh[idx_fh_F_ord4(iF-1,jF,kF+1,ex)]+F8*fh[idx_fh_F_ord4(iF+1,jF,kF+1,ex)]-fh[idx_fh_F_ord4(iF+2,jF,kF+1,ex)]);
                            const double t_kp2 = (fh[idx_fh_F_ord4(iF-2,jF,kF+2,ex)]-F8*fh[idx_fh_F_ord4(iF-1,jF,kF+2,ex)]+F8*fh[idx_fh_F_ord4(iF+1,jF,kF+2,ex)]-fh[idx_fh_F_ord4(iF+2,jF,kF+2,ex)]);
                            fxz[p] = Fdxdz * (t_km2 - F8*t_km1 + F8*t_kp1 - t_kp2);
                        }
                        {
                            const double t_km2 = (fh[idx_fh_F_ord4(iF,jF-2,kF-2,ex)]-F8*fh[idx_fh_F_ord4(iF,jF-1,kF-2,ex)]+F8*fh[idx_fh_F_ord4(iF,jF+1,kF-2,ex)]-fh[idx_fh_F_ord4(iF,jF+2,kF-2,ex)]);
                            const double t_km1 = (fh[idx_fh_F_ord4(iF,jF-2,kF-1,ex)]-F8*fh[idx_fh_F_ord4(iF,jF-1,kF-1,ex)]+F8*fh[idx_fh_F_ord4(iF,jF+1,kF-1,ex)]-fh[idx_fh_F_ord4(iF,jF+2,kF-1,ex)]);
                            const double t_kp1 = (fh[idx_fh_F_ord4(iF,jF-2,kF+1,ex)]-F8*fh[idx_fh_F_ord4(iF,jF-1,kF+1,ex)]+F8*fh[idx_fh_F_ord4(iF,jF+1,kF+1,ex)]-fh[idx_fh_F_ord4(iF,jF+2,kF+1,ex)]);
                            const double t_kp2 = (fh[idx_fh_F_ord4(iF,jF-2,kF+2,ex)]-F8*fh[idx_fh_F_ord4(iF,jF-1,kF+2,ex)]+F8*fh[idx_fh_F_ord4(iF,jF+1,kF+2,ex)]-fh[idx_fh_F_ord4(iF,jF+2,kF+2,ex)]);
                            fyz[p] = Fdydz * (t_km2 - F8*t_km1 + F8*t_kp1 - t_kp2);
                        }
                    }
                }
            }
        }

        /* 6th-order: skip 8th overlap */
        if (has6) {
            for (int k0 = k6_lo; k0 <= k6_hi; ++k0) {
                const int kF = k0 + 3;
                for (int j0 = j6_lo; j0 <= j6_hi; ++j0) {
                    const int jF = j0 + 3;
                    for (int i0 = i6_lo; i0 <= i6_hi; ++i0) {
                        if (has8 && i0>=i8_lo && i0<=i8_hi && j0>=j8_lo && j0<=j8_hi && k0>=k8_lo && k0<=k8_hi) continue;
                        const int iF = i0 + 3;
                        const size_t p = idx_ex(i0, j0, k0, ex);

                        fxx[p] = Xdxdx * (
                            TWO  * fh[idx_fh_F_ord4(iF-3,jF,kF,ex)] - F27*fh[idx_fh_F_ord4(iF-2,jF,kF,ex)] + F270*fh[idx_fh_F_ord4(iF-1,jF,kF,ex)] - F490*fh[idx_fh_F_ord4(iF,jF,kF,ex)] + F270*fh[idx_fh_F_ord4(iF+1,jF,kF,ex)] - F27*fh[idx_fh_F_ord4(iF+2,jF,kF,ex)] + TWO*fh[idx_fh_F_ord4(iF+3,jF,kF,ex)]);
                        fyy[p] = Xdydy * (
                            TWO  * fh[idx_fh_F_ord4(iF,jF-3,kF,ex)] - F27*fh[idx_fh_F_ord4(iF,jF-2,kF,ex)] + F270*fh[idx_fh_F_ord4(iF,jF-1,kF,ex)] - F490*fh[idx_fh_F_ord4(iF,jF,kF,ex)] + F270*fh[idx_fh_F_ord4(iF,jF+1,kF,ex)] - F27*fh[idx_fh_F_ord4(iF,jF+2,kF,ex)] + TWO*fh[idx_fh_F_ord4(iF,jF+3,kF,ex)]);
                        fzz[p] = Xdzdz * (
                            TWO  * fh[idx_fh_F_ord4(iF,jF,kF-3,ex)] - F27*fh[idx_fh_F_ord4(iF,jF,kF-2,ex)] + F270*fh[idx_fh_F_ord4(iF,jF,kF-1,ex)] - F490*fh[idx_fh_F_ord4(iF,jF,kF,ex)] + F270*fh[idx_fh_F_ord4(iF,jF,kF+1,ex)] - F27*fh[idx_fh_F_ord4(iF,jF,kF+2,ex)] + TWO*fh[idx_fh_F_ord4(iF,jF,kF+3,ex)]);

                        {
                            #define XSTEN6_8(JF, KF_DUMMY) \
                                (-fh[idx_fh_F_ord4(iF-3,JF,KF_DUMMY,ex)] + F9*fh[idx_fh_F_ord4(iF-2,JF,KF_DUMMY,ex)] - F45*fh[idx_fh_F_ord4(iF-1,JF,KF_DUMMY,ex)] + F45*fh[idx_fh_F_ord4(iF+1,JF,KF_DUMMY,ex)] - F9*fh[idx_fh_F_ord4(iF+2,JF,KF_DUMMY,ex)] + fh[idx_fh_F_ord4(iF+3,JF,KF_DUMMY,ex)])
                            fxy[p] = Xdxdy * (
                                -XSTEN6_8(jF-3,kF) + F9*XSTEN6_8(jF-2,kF) - F45*XSTEN6_8(jF-1,kF) + F45*XSTEN6_8(jF+1,kF) - F9*XSTEN6_8(jF+2,kF) + XSTEN6_8(jF+3,kF));
                            fxz[p] = Xdxdz * (
                                -XSTEN6_8(jF,kF-3) + F9*XSTEN6_8(jF,kF-2) - F45*XSTEN6_8(jF,kF-1) + F45*XSTEN6_8(jF,kF+1) - F9*XSTEN6_8(jF,kF+2) + XSTEN6_8(jF,kF+3));
                            #undef XSTEN6_8
                        }
                        {
                            #define YSTEN6_8(JF, KF_DUMMY) \
                                (-fh[idx_fh_F_ord4(iF,JF-3,KF_DUMMY,ex)] + F9*fh[idx_fh_F_ord4(iF,JF-2,KF_DUMMY,ex)] - F45*fh[idx_fh_F_ord4(iF,JF-1,KF_DUMMY,ex)] + F45*fh[idx_fh_F_ord4(iF,JF+1,KF_DUMMY,ex)] - F9*fh[idx_fh_F_ord4(iF,JF+2,KF_DUMMY,ex)] + fh[idx_fh_F_ord4(iF,JF+3,KF_DUMMY,ex)])
                            fyz[p] = Xdydz * (
                                -YSTEN6_8(jF,kF-3) + F9*YSTEN6_8(jF,kF-2) - F45*YSTEN6_8(jF,kF-1) + F45*YSTEN6_8(jF,kF+1) - F9*YSTEN6_8(jF,kF+2) + YSTEN6_8(jF,kF+3));
                            #undef YSTEN6_8
                        }
                    }
                }
            }
        }

        /* 8th-order: interior only */
        if (has8) {
            for (int k0 = k8_lo; k0 <= k8_hi; ++k0) {
                const int kF = k0 + 3;
                for (int j0 = j8_lo; j0 <= j8_hi; ++j0) {
                    const int jF = j0 + 3;
                    for (int i0 = i8_lo; i0 <= i8_hi; ++i0) {
                        const int iF = i0 + 3;
                        const size_t p = idx_ex(i0, j0, k0, ex);

                        /* Diagonal: [-9,+128,-1008,+8064,-14350,+8064,-1008,+128,-9] / (5040*dx^2) */
                        fxx[p] = Edxdx * (
                            -(double)9 * fh[idx_fh_F_ord4(iF - 4, jF,     kF,     ex)] +
                            F128  * fh[idx_fh_F_ord4(iF - 3, jF,     kF,     ex)] -
                            F1008 * fh[idx_fh_F_ord4(iF - 2, jF,     kF,     ex)] +
                            F8064 * fh[idx_fh_F_ord4(iF - 1, jF,     kF,     ex)] -
                            F14350* fh[idx_fh_F_ord4(iF,     jF,     kF,     ex)] +
                            F8064 * fh[idx_fh_F_ord4(iF + 1, jF,     kF,     ex)] -
                            F1008 * fh[idx_fh_F_ord4(iF + 2, jF,     kF,     ex)] +
                            F128  * fh[idx_fh_F_ord4(iF + 3, jF,     kF,     ex)] -
                            (double)9 * fh[idx_fh_F_ord4(iF + 4, jF,     kF,     ex)]);

                        fyy[p] = Edydy * (
                            -(double)9 * fh[idx_fh_F_ord4(iF,     jF - 4, kF,     ex)] +
                            F128  * fh[idx_fh_F_ord4(iF,     jF - 3, kF,     ex)] -
                            F1008 * fh[idx_fh_F_ord4(iF,     jF - 2, kF,     ex)] +
                            F8064 * fh[idx_fh_F_ord4(iF,     jF - 1, kF,     ex)] -
                            F14350* fh[idx_fh_F_ord4(iF,     jF,     kF,     ex)] +
                            F8064 * fh[idx_fh_F_ord4(iF,     jF + 1, kF,     ex)] -
                            F1008 * fh[idx_fh_F_ord4(iF,     jF + 2, kF,     ex)] +
                            F128  * fh[idx_fh_F_ord4(iF,     jF + 3, kF,     ex)] -
                            (double)9 * fh[idx_fh_F_ord4(iF,     jF + 4, kF,     ex)]);

                        fzz[p] = Edzdz * (
                            -(double)9 * fh[idx_fh_F_ord4(iF,     jF,     kF - 4, ex)] +
                            F128  * fh[idx_fh_F_ord4(iF,     jF,     kF - 3, ex)] -
                            F1008 * fh[idx_fh_F_ord4(iF,     jF,     kF - 2, ex)] +
                            F8064 * fh[idx_fh_F_ord4(iF,     jF,     kF - 1, ex)] -
                            F14350* fh[idx_fh_F_ord4(iF,     jF,     kF,     ex)] +
                            F8064 * fh[idx_fh_F_ord4(iF,     jF,     kF + 1, ex)] -
                            F1008 * fh[idx_fh_F_ord4(iF,     jF,     kF + 2, ex)] +
                            F128  * fh[idx_fh_F_ord4(iF,     jF,     kF + 3, ex)] -
                            (double)9 * fh[idx_fh_F_ord4(iF,     jF,     kF + 4, ex)]);

                        /* Mixed: 9x9 outer product.
                           x-stencil: +3*f(i-4)-32*f(i-3)+168*f(i-2)-672*f(i-1)+672*f(i+1)-168*f(i+2)+32*f(i+3)-3*f(i+4)
                           y/z weights: same [+3,-32,+168,-672,+672,-168,+32,-3] / 705600 */
                        {
                            #define XSTEN8(JF, KF_DUMMY) \
                                (+(double)3*fh[idx_fh_F_ord4(iF-4,JF,KF_DUMMY,ex)] - F32*fh[idx_fh_F_ord4(iF-3,JF,KF_DUMMY,ex)] + F168*fh[idx_fh_F_ord4(iF-2,JF,KF_DUMMY,ex)] - F672*fh[idx_fh_F_ord4(iF-1,JF,KF_DUMMY,ex)] + F672*fh[idx_fh_F_ord4(iF+1,JF,KF_DUMMY,ex)] - F168*fh[idx_fh_F_ord4(iF+2,JF,KF_DUMMY,ex)] + F32*fh[idx_fh_F_ord4(iF+3,JF,KF_DUMMY,ex)] - (double)3*fh[idx_fh_F_ord4(iF+4,JF,KF_DUMMY,ex)])

                            fxy[p] = Edxdy * (
                                +(double)3*XSTEN8(jF-4,kF) - F32*XSTEN8(jF-3,kF) + F168*XSTEN8(jF-2,kF) - F672*XSTEN8(jF-1,kF) + F672*XSTEN8(jF+1,kF) - F168*XSTEN8(jF+2,kF) + F32*XSTEN8(jF+3,kF) - (double)3*XSTEN8(jF+4,kF));

                            fxz[p] = Edxdz * (
                                +(double)3*XSTEN8(jF,kF-4) - F32*XSTEN8(jF,kF-3) + F168*XSTEN8(jF,kF-2) - F672*XSTEN8(jF,kF-1) + F672*XSTEN8(jF,kF+1) - F168*XSTEN8(jF,kF+2) + F32*XSTEN8(jF,kF+3) - (double)3*XSTEN8(jF,kF+4));

                            #undef XSTEN8
                        }
                        {
                            #define YSTEN8(JF, KF_DUMMY) \
                                (+(double)3*fh[idx_fh_F_ord4(iF,JF-4,KF_DUMMY,ex)] - F32*fh[idx_fh_F_ord4(iF,JF-3,KF_DUMMY,ex)] + F168*fh[idx_fh_F_ord4(iF,JF-2,KF_DUMMY,ex)] - F672*fh[idx_fh_F_ord4(iF,JF-1,KF_DUMMY,ex)] + F672*fh[idx_fh_F_ord4(iF,JF+1,KF_DUMMY,ex)] - F168*fh[idx_fh_F_ord4(iF,JF+2,KF_DUMMY,ex)] + F32*fh[idx_fh_F_ord4(iF,JF+3,KF_DUMMY,ex)] - (double)3*fh[idx_fh_F_ord4(iF,JF+4,KF_DUMMY,ex)])

                            fyz[p] = Edydz * (
                                +(double)3*YSTEN8(jF,kF-4) - F32*YSTEN8(jF,kF-3) + F168*YSTEN8(jF,kF-2) - F672*YSTEN8(jF,kF-1) + F672*YSTEN8(jF,kF+1) - F168*YSTEN8(jF,kF+2) + F32*YSTEN8(jF,kF+3) - (double)3*YSTEN8(jF,kF+4));

                            #undef YSTEN8
                        }
                    }
                }
            }
        }
        return;
    }
#else
#error "fdderivs_c.C: unsupported ghost_width (must be 2, 3, 4, or 5)"
#endif
}
