#include "macrodef.h"
#include "share_func.h"
#include <cstddef>

/* Forward declarations — Fortran-mangled names from shell C kernels */
extern "C" {
void fderivs_sh_(const int ex[3], const double *f,
                 double *fx, double *fy, double *fz,
                 const double *X, const double *Y, const double *Z,
                 double SYM1, double SYM2, double SYM3,
                 int Symmetry, int onoff, int sst);

void fdderivs_sh_(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, int sst);

void fdderivs_shc_(int *ex,
                    double *f,
                    double *fxx, double *fxy, double *fxz,
                    double *fyy, double *fyz, double *fzz,
                    double *crho, double *sigma, double *R,
                    double &SYM1, double &SYM2, double &SYM3,
                    int &Symmetry, int &Lev, int &sst,
                    double *drhodx, double *drhody, double *drhodz,
                    double *dsigmadx, double *dsigmady, double *dsigmadz,
                    double *dRdx, double *dRdy, double *dRdz,
                    double *drhodxx, double *drhodxy, double *drhodxz,
                    double *drhodyy, double *drhodyz, double *drhodzz,
                    double *dsigmadxx, double *dsigmadxy, double *dsigmadxz,
                    double *dsigmadyy, double *dsigmadyz, double *dsigmadzz,
                    double *dRdxx, double *dRdxy, double *dRdxz,
                    double *dRdyy, double *dRdyz, double *dRdzz)
{
    const int ex3[3] = { ex[0], ex[1], ex[2] };
    const size_t n = (size_t)ex[0] * (size_t)ex[1] * (size_t)ex[2];

    double *gx  = (double*)malloc(n * sizeof(double));
    double *gy  = (double*)malloc(n * sizeof(double));
    double *gz  = (double*)malloc(n * sizeof(double));
    double *gxx = (double*)malloc(n * sizeof(double));
    double *gxy = (double*)malloc(n * sizeof(double));
    double *gxz = (double*)malloc(n * sizeof(double));
    double *gyy = (double*)malloc(n * sizeof(double));
    double *gyz = (double*)malloc(n * sizeof(double));
    double *gzz = (double*)malloc(n * sizeof(double));

    if (!gx||!gy||!gz||!gxx||!gxy||!gxz||!gyy||!gyz||!gzz) {
        free(gx);free(gy);free(gz);free(gxx);free(gxy);free(gxz);free(gyy);free(gyz);free(gzz);
        return;
    }

    fderivs_sh_(ex3, f, gx, gy, gz, crho, sigma, R, SYM1, SYM2, SYM3, Symmetry, Lev, sst);
    fdderivs_sh_(ex3, f, gxx, gxy, gxz, gyy, gyz, gzz, crho, sigma, R, SYM1, SYM2, SYM3, Symmetry, Lev, sst);

    for (size_t i = 0; i < n; ++i) {
        const double rx=drhodx[i], ry=drhody[i], rz=drhodz[i];
        const double sx=dsigmadx[i], sy=dsigmady[i], sz=dsigmadz[i];
        const double Rx=dRdx[i], Ry=dRdy[i], Rz=dRdz[i];
        const double rxx=drhodxx[i], rxy=drhodxy[i], rxz=drhodxz[i];
        const double ryy=drhodyy[i], ryz=drhodyz[i], rzz=drhodzz[i];
        const double sxx=dsigmadxx[i], sxy=dsigmadxy[i], sxz=dsigmadxz[i];
        const double syy=dsigmadyy[i], syz=dsigmadyz[i], szz=dsigmadzz[i];
        const double Rxx=dRdxx[i], Rxy=dRdxy[i], Rxz=dRdxz[i];
        const double Ryy=dRdyy[i], Ryz=dRdyz[i], Rzz=dRdzz[i];

        const double Gr=gx[i], Gs=gy[i], GR=gz[i];
        const double Grr=gxx[i], Grs=gxy[i], GrR=gxz[i];
        const double Gss=gyy[i], GsR=gyz[i], GRR=gzz[i];

        /* fxx */
        fxx[i] = rx*rx*Grr + sx*sx*Gss + Rx*Rx*GRR
               + 2.0*(rx*sx*Grs + rx*Rx*GrR + sx*Rx*GsR)
               + rxx*Gr + sxx*Gs + Rxx*GR;

        /* fxy */
        fxy[i] = rx*ry*Grr + sx*sy*Gss + Rx*Ry*GRR
               + rx*sy*Grs + ry*sx*Grs + rx*Ry*GrR + ry*Rx*GrR + sx*Ry*GsR + sy*Rx*GsR
               + rxy*Gr + sxy*Gs + Rxy*GR;

        /* fxz */
        fxz[i] = rx*rz*Grr + sx*sz*Gss + Rx*Rz*GRR
               + rx*sz*Grs + rz*sx*Grs + rx*Rz*GrR + rz*Rx*GrR + sx*Rz*GsR + sz*Rx*GsR
               + rxz*Gr + sxz*Gs + Rxz*GR;

        /* fyy */
        fyy[i] = ry*ry*Grr + sy*sy*Gss + Ry*Ry*GRR
               + 2.0*(ry*sy*Grs + ry*Ry*GrR + sy*Ry*GsR)
               + ryy*Gr + syy*Gs + Ryy*GR;

        /* fyz */
        fyz[i] = ry*rz*Grr + sy*sz*Gss + Ry*Rz*GRR
               + ry*sz*Grs + rz*sy*Grs + ry*Rz*GrR + rz*Ry*GrR + sy*Rz*GsR + sz*Ry*GsR
               + ryz*Gr + syz*Gs + Ryz*GR;

        /* fzz */
        fzz[i] = rz*rz*Grr + sz*sz*Gss + Rz*Rz*GRR
               + 2.0*(rz*sz*Grs + rz*Rz*GrR + sz*Rz*GsR)
               + rzz*Gr + szz*Gs + Rzz*GR;
    }

    free(gx);free(gy);free(gz);free(gxx);free(gxy);free(gxz);free(gyy);free(gyz);free(gzz);
}

} // extern "C"