Как ускорить математический код C#

У меня есть код трехмерной интерполяции, который занимает 90% времени выполнения моих проектов и не может быть предварительно вычислен.

Какие методы можно использовать, чтобы ускорить это? Алгоритмическая или микрооптимизация?

Вот код для интересующихся.

Он в основном берет данные, которые были помещены в массивы 2 3d, и интерполирует остальные данные.

РЕДАКТИРОВАТЬ: Кроме того, я уже разбиваю это на потоки на более высоком уровне для повышения производительности, но это не помогает на Windows Phone, поскольку все они одноядерные...

Я, вероятно, сделаю что-то вроде ( Single[] DensityMap = new Single[128 * 128 * 128];), чтобы удалить попадание в массив multi-D. Я обращаюсь к массиву в сотнях мест и надеялся, что мне не придется этого делать (обтекание функцией не помогает, так как телефон Windows не будет встраивать вызов функции, и тогда это не поможет...)

float[, ,] DensityMap = new float[128, 128, 128];
float[, ,] PressureMap = new float[128, 128, 128];

unchecked
{
    for (int x = 0; x < g_CraftWorldConstants.RegionSizeX; x++)
    {
        int offsetX = (x / SAMPLE_RATE_3D_HOR) * SAMPLE_RATE_3D_HOR;
        int plusOffsetX = SAMPLE_RATE_3D_HOR + offsetX;
        int poxox = plusOffsetX - offsetX;
        double poxxpoxox = ((plusOffsetX - x) / (double)poxox);
        double xoxpoxox = ((x - offsetX) / (double)poxox);

        for (int y = 0; y < g_CraftWorldSettings.GET.RegionSizeY; y++)
        {
            int offsetY = (y / SAMPLE_RATE_3D_VERT) * SAMPLE_RATE_3D_VERT;
            int plusOffsetY = SAMPLE_RATE_3D_VERT + offsetY;
            int poyoy = plusOffsetY - offsetY;
            double poyypoyoy = ((plusOffsetY - y) / (double)poyoy);
            double yoypoyoy = ((y - offsetY) / (double)poyoy);

            for (int z = 0; z < g_CraftWorldConstants.RegionSizeZ; z++)
            {
                if (!(x % SAMPLE_RATE_3D_HOR == 0 && y % SAMPLE_RATE_3D_VERT == 0 && z % SAMPLE_RATE_3D_HOR == 0))
                {
                    int offsetZ = (z / SAMPLE_RATE_3D_HOR) * SAMPLE_RATE_3D_HOR;
                    int plusOffsetZ = SAMPLE_RATE_3D_HOR + offsetZ;
                    int pozoz = plusOffsetZ - offsetZ;
                    double pozzpozoz = ((plusOffsetZ - z) / (double)pozoz);
                    double zozpozoz = ((z - offsetZ) / (double)pozoz);

                    double x00 = poxxpoxox * in_DensityMap[offsetX, offsetY, offsetZ] + xoxpoxox * in_DensityMap[plusOffsetX, offsetY, offsetZ];
                    double x10 = poxxpoxox * in_DensityMap[offsetX, offsetY, plusOffsetZ] + xoxpoxox * in_DensityMap[plusOffsetX, offsetY, plusOffsetZ];
                    double x01 = poxxpoxox * in_DensityMap[offsetX, plusOffsetY, offsetZ] + xoxpoxox * in_DensityMap[plusOffsetX, plusOffsetY, offsetZ];
                    double x11 = poxxpoxox * in_DensityMap[offsetX, plusOffsetY, plusOffsetZ] + xoxpoxox * in_DensityMap[plusOffsetX, plusOffsetY, plusOffsetZ];

                    double r0 = poyypoyoy * x00 + yoypoyoy * x01;
                    double r1 = poyypoyoy * x10 + yoypoyoy * x11;
                    in_DensityMap[x, y, z] = (float)(pozzpozoz * r0 + zozpozoz * r1);

                    double x02 = poxxpoxox * in_CaveDensity[offsetX, offsetY, offsetZ] + xoxpoxox * in_CaveDensity[plusOffsetX, offsetY, offsetZ];
                    double x12 = poxxpoxox * in_CaveDensity[offsetX, offsetY, plusOffsetZ] + xoxpoxox * in_CaveDensity[plusOffsetX, offsetY, plusOffsetZ];
                    double x03 = poxxpoxox * in_CaveDensity[offsetX, plusOffsetY, offsetZ] + xoxpoxox * in_CaveDensity[plusOffsetX, plusOffsetY, offsetZ];
                    double x13 = poxxpoxox * in_CaveDensity[offsetX, plusOffsetY, plusOffsetZ] + xoxpoxox * in_CaveDensity[plusOffsetX, plusOffsetY, plusOffsetZ];

                    double r2 = poyypoyoy * x02 + yoypoyoy * x03;
                    double r3 = poyypoyoy * x12 + yoypoyoy * x13;
                    in_CaveDensity[x, y, z] = (float)(pozzpozoz * r2 + zozpozoz * r3);
                }
            }
        }
    }
}