Optimize the SH code

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

 }
 
 // i = l * (l + 1) + m
+// See PremultiplySH().
-float EvaluateSphericalHarmonics(int i, float3 d)
+float EvaluatePremultipliedSH(int i, float3 d)
-        case 0: return rsqrt(4 * PI);
-        case 1: return -sqrt(0.75 / PI) * d.y;
-        case 2: return  sqrt(0.75 / PI) * d.z;
-        case 3: return -sqrt(0.75 / PI) * d.x;
-        case 4: return  sqrt(3.75 / PI) * d.x * d.y;
-        case 5: return -sqrt(3.75 / PI) * d.y * d.z;
-        case 6: return  sqrt(0.3125 / PI) * (3 * Sq(d.z) - 1);
-        case 7: return -sqrt(3.75 / PI) * d.x * d.z;
-        case 8: return  sqrt(0.9375 / PI) * (Sq(d.x) - Sq(d.y));
+        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
     }
 }
 
     for (int i = 0; i < 9; i++)
     {
-        radiance += probe.coeffs[i] * EvaluateSphericalHarmonics(i, d);
+        radiance += probe.coeffs[i] * EvaluatePremultipliedSH(i, d);
-    return max(0, radiance);
+    return radiance; // Warning: could be negative
 }
 
 // Computes the light integral (in-scattered radiance) within the voxel.

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

         return zh;
     }
 
+    // Ref: "Stupid Spherical Harmonics Tricks", p. 6.
+    public static unsafe SphericalHarmonicsL2 Convolve(SphericalHarmonicsL2 sh, ZonalHarmonicsL2 zh)
+    {
+        for (int l = 0; l <= 2; l++)
+        {
+            float n = Mathf.Sqrt((4.0f * Mathf.PI) / (2 * l + 1));
+            float k = zh.coeffs[l];
+            float p = n * k;
+
+            for (int m = -l; m <= l; m++)
+            {
+                int i = l * (l + 1) + m;
+
+                for (int c = 0; c < 3; c++)
+                {
+                    sh[c, i] *= p;
+                }
+            }
+        }
+
+        return sh;
+    }
+
-    public static unsafe SphericalHarmonicsL2 DenormalizeSH(SphericalHarmonicsL2 sh)
+    public static SphericalHarmonicsL2 DenormalizeSH(SphericalHarmonicsL2 sh)
-        float c0 = 1.0f / (2.0f * sqrtPi);
-        float c1 = Mathf.Sqrt( 3.0f) / ( 3.0f * sqrtPi);
-        float c2 = Mathf.Sqrt(15.0f) / ( 8.0f * sqrtPi);
-        float c3 = Mathf.Sqrt( 5.0f) / (16.0f * sqrtPi);
-        float c4 = 0.5f * c2;
+        const float c0 = 0.28209479177387814347f; // 1/2  * sqrt(1/Pi)
+        const float c1 = 0.32573500793527994772f; // 1/3  * sqrt(3/Pi)
+        const float c2 = 0.27313710764801976764f; // 1/8  * sqrt(15/Pi)
+        const float c3 = 0.07884789131313000151f; // 1/16 * sqrt(5/Pi)
+        const float c4 = 0.13656855382400988382f; // 1/16 * sqrt(15/Pi)
-        float[] invNormConsts = { 1.0f / c0, -1.0f / c1, 1.0f / c1, -1.0f / c1, 1.0f / c2, -1.0f / c2, 1.0f / c3, -1.0f / c2, 1.0f / c4 };
+        float[] invNormConsts ={ 1.0f / c0, -1.0f / c1, 1.0f / c1, -1.0f / c1, 1.0f / c2, -1.0f / c2, 1.0f / c3, -1.0f / c2, 1.0f / c4 };
-        for (int i = 0; i < 9; i++)
+        for (int c = 0; c < 3; c++)
-            for (int c = 0; c < 3; c++)
+            for (int i = 0; i < 9; i++)
             {
                 sh[c, i] *= invNormConsts[i];
             }
     }
 
-    // Ref: "Stupid Spherical Harmonics Tricks", p. 6.
-    public static unsafe SphericalHarmonicsL2 Convolve(SphericalHarmonicsL2 sh, ZonalHarmonicsL2 zh)
+    // Premultiplies the SH with the polynomial coefficients of SH basis functions,
+    // which avoids using any constants during SH evaluation.
+    public static SphericalHarmonicsL2 PremultiplySH(SphericalHarmonicsL2 sh)
-        for (int l = 0; l <= 2; l++)
+        const float k0 = 0.28209479177387814347f; // {0, 0} : 1/2 * sqrt(1/Pi)
+        const float k1 = 0.48860251190291992159f; // {1, 0} : 1/2 * sqrt(3/Pi)
+        const float k2 = 1.09254843059207907054f; // {2,-2} : 1/2 * sqrt(15/Pi)
+        const float k3 = 0.31539156525252000603f; // {2, 0} : 1/4 * sqrt(5/Pi)
+        const float k4 = 0.54627421529603953527f; // {2, 2} : 1/4 * sqrt(15/Pi)
+
+        float[] ks ={ k0, -k1, k1, -k1, k2, -k2, k3, -k2, k4 };
+
+        for (int c = 0; c < 3; c++)
-            float n = Mathf.Sqrt((4.0f * Mathf.PI) / (2 * l + 1));
-            float k = zh.coeffs[l];
-            float p = n * k;
-
-            for (int m = -l; m <= l; m++)
+            for (int i = 0; i < 9; i++)
-                int i = l * (l + 1) + m;
-
-                sh[0, i] *= p;
-                sh[1, i] *= p;
-                sh[2, i] *= p;
+                sh[c, i] *= ks[i];
             }
         }
 
     {
         SphericalHarmonicsL2 probeSH = DenormalizeSH(RenderSettings.ambientProbe);
         ZonalHarmonicsL2     phaseZH = GetCornetteShanksPhaseFunction(asymmetry);
-
-        SphericalHarmonicsL2 convolvedSH = Convolve(probeSH, phaseZH);
+        SphericalHarmonicsL2 finalSH = PremultiplySH(Convolve(probeSH, phaseZH));
-            coeffs[i].x = convolvedSH[0, i]; // R
-            coeffs[i].y = convolvedSH[1, i]; // G
-            coeffs[i].z = convolvedSH[2, i]; // B
-            coeffs[i].w = 0;                 // Unused
+            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
         }
 
         cmd.SetGlobalVectorArray(HDShaderIDs._AmbientProbeCoeffs, coeffs);