 //! The developer has chosen not to document this variable
 __device__ __constant__ float GPU_rho_coeff[CONSTANT_SIZE];
 
+//! Coefficients of a power expansion of sin(x)/x
+const Scalar sinc_coeff[] = {Scalar(1.0), Scalar(-1.0/6.0), Scalar(1.0/120.0),
+                        Scalar(-1.0/5040.0),Scalar(1.0/362880.0),
+                        Scalar(-1.0/39916800.0)};
+
+//! Coefficients of a power expansion of sin(x)/x
+__device__ __constant__ float gpu_sinc_coeff[6];
+
 /*! \file PPPMForceGPU.cu
   \brief Defines GPU kernel code for calculating the Fourier space forces for the Coulomb interaction. Used by PPPMForceComputeGPU.
 */

                qx = (j.x+(float)(Nx*ix));
                kn1 = b1 * qx;
                argx = Scalar(0.5)*qx*kH.x;
                if (argx != 0.0f)
                    wx = powf(sinf(argx)/argx,order);
                else
                    wx = Scalar(1.0);

               for (iy = -nby; iy <= nby; iy++) {
                    qy = (j.y+(float)(Ny*iy));
                    kn2 = b2 * qy;
                    argy = Scalar(0.5)*qy*kH.y;
                    if (argy != 0.0f)
                        wy = powf(sinf(argy)/argy,order);
                    else
                        wy = Scalar(1.0);

                    for (iz = -nbz; iz <= nbz; iz++) {
                        qz = (j.z+(float)(Nz*iz));

                        kn = kn1+kn2+kn3;

                        argz = Scalar(0.5)*qz*kH.z;
                        wz = powf(gpu_sinc(argz),order);

                        dot1 = dot(kn,k);
                        dot2 = dot(kn,kn);
