Further optimize the SH code

ScriptableRenderPipeline/HDRenderPipeline/HDRP/Lighting/Volumetrics/Resources/VolumetricLighting.compute

     return probe;
 }
 
-// i = l * (l + 1) + m
-// See PremultiplySH().
-// Ref: "Stupid Spherical Harmonics Tricks", p. 28-29.
-float EvaluatePremultipliedSH(int i, float3 d)
-{
-    switch (i)
-    {
-        case 0: return 1;
-        case 1: return d.y;
-        case 2: return d.z;
-        case 3: return d.x;
-        case 4: return d.x * d.y;
-        case 5: return d.y * d.z;
-        case 6: return (3 * Sq(d.z) - 1);
-        case 7: return d.x * d.z;
-        case 8: return (Sq(d.x) - Sq(d.y));
-        default: return FLT_INF; // Unreachable
-    }
-}
-
+// Warning: since it's a SH, the results can be negative!
-    for (int i = 0; i < 9; i++)
-    {
-        radiance += probe.coeffs[i] * EvaluatePremultipliedSH(i, d);
-    }
-
-    return radiance; // Warning: could be negative
+    // Coefficients are premultiplied and reordered in the MAD form:
+    // HornerForm[c_0 + c_1 y + c_2 z + c_3 x + c_4 x y + c_5 y z + c_6 (3 z^2 - 1) + c_7 x z + c_8 (x^2 - y^2)]
+    // = z (3 c_6 z + c_7 x + c_2) + y (-c_8 y + c_5 z + c_4 x + c_1) + x (c_8 x + c_3) + (c_0 - c_6)
+    return probe.coeffs[0] + d.x * (probe.coeffs[1] + d.x * probe.coeffs[2])
+                           + d.y * (probe.coeffs[3] + d.x * probe.coeffs[4] + d.z * probe.coeffs[5] - d.y * probe.coeffs[2])
+                           + d.z * (probe.coeffs[6] + d.x * probe.coeffs[7] + d.z * probe.coeffs[8]);
 }
 
 // Computes the light integral (in-scattered radiance) within the voxel.

ScriptableRenderPipeline/HDRenderPipeline/HDRP/Lighting/Volumetrics/VolumetricLighting.cs

 
     // Premultiplies the SH with the polynomial coefficients of SH basis functions,
     // which avoids using any constants during SH evaluation.
+    // The resulting evaluation takes the form:
+    // c_0 + c_1 y + c_2 z + c_3 x + c_4 x y + c_5 y z + c_6 (3 z^2 - 1) + c_7 x z + c_8 (x^2 - y^2)
     public static SphericalHarmonicsL2 PremultiplySH(SphericalHarmonicsL2 sh)
     {
         const float k0 = 0.28209479177387814347f; // {0, 0} : 1/2 * sqrt(1/Pi)
         ZonalHarmonicsL2     phaseZH = GetCornetteShanksPhaseFunction(asymmetry);
         SphericalHarmonicsL2 finalSH = PremultiplySH(Convolve(probeSH, phaseZH));
 
+        // Reorder coefficients in the MAD form:
+        // HornerForm[c_0 + c_1 y + c_2 z + c_3 x + c_4 x y + c_5 y z + c_6 (3 z^2 - 1) + c_7 x z + c_8 (x^2 - y^2)]
+        // = z (3 c_6 z + c_7 x + c_2) + y (-c_8 y + c_5 z + c_4 x + c_1) + x (c_8 x + c_3) + (c_0 - c_6)
-        for (int i = 0; i < 9; i++)
-        {
-            coeffs[i].x = finalSH[0, i]; // R
-            coeffs[i].y = finalSH[1, i]; // G
-            coeffs[i].z = finalSH[2, i]; // B
-            coeffs[i].w = 0;             // Unused
-        }
+        const int r = 0, g = 1, b = 2;
+
+        coeffs[0] = new Vector3(finalSH[r, 0], finalSH[g, 0], finalSH[b, 0])
+                  - new Vector3(finalSH[r, 6], finalSH[g, 6], finalSH[b, 6]);
+        coeffs[1] = new Vector3(finalSH[r, 3], finalSH[g, 3], finalSH[b, 3]);
+        coeffs[2] = new Vector3(finalSH[r, 8], finalSH[g, 8], finalSH[b, 8]);
+        coeffs[3] = new Vector3(finalSH[r, 1], finalSH[g, 1], finalSH[b, 1]);
+        coeffs[4] = new Vector3(finalSH[r, 4], finalSH[g, 4], finalSH[b, 4]);
+        coeffs[5] = new Vector3(finalSH[r, 5], finalSH[g, 5], finalSH[b, 5]);
+        // Avoid reduplicating c_8.
+        coeffs[6] = new Vector3(finalSH[r, 2], finalSH[g, 2], finalSH[b, 2]);
+        coeffs[7] = new Vector3(finalSH[r, 7], finalSH[g, 7], finalSH[b, 7]);
+        coeffs[8] = new Vector3(finalSH[r, 6], finalSH[g, 6], finalSH[b, 6]) * 3.0f;
 
         cmd.SetGlobalVectorArray(HDShaderIDs._AmbientProbeCoeffs, coeffs);
     }