Optimize and add Mathematica notes

Assets/ScriptableRenderPipeline/ShaderLibrary/AreaLighting.hlsl

 // N.b.: this function accounts for horizon clipping.
 float DiffuseSphereLightIrradiance(float sinSqSigma, float cosOmega)
 {
-    float irradiance;
-
 #if 0 // Ref: Area Light Sources for Real-Time Graphics, page 4 (1996).
     float sinSqOmega = saturate(1 - cosOmega * cosOmega);
     float cosSqSigma = saturate(1 - sinSqSigma);
     if (omega < HALF_PI - sigma)
     {
         // No horizon occlusion (case #1).
-        irradiance = e;
+        return e;
     }
     else
     {
         if (omega < HALF_PI)
         {
             // Partial horizon occlusion (case #2).
-            irradiance = e + INV_PI * (g - h);
+            return saturate(e + INV_PI * (g - h));
-            irradiance = INV_PI * (g + h);
+            return saturate(INV_PI * (g + h));
-#else // Ref: Moving Frostbite to Physically Based Rendering, page 47 (2015).
-    float cosSqOmega = cosOmega * cosOmega;
+#else // Ref: Moving Frostbite to Physically Based Rendering, page 47 (2015, optimized).
+    float cosSqOmega = cosOmega * cosOmega;                     // y^2
-    if (cosSqOmega > sinSqSigma)
+    if (cosSqOmega > sinSqSigma)                                // (y^2)>x
-        irradiance = sinSqSigma * saturate(cosOmega);
+        return sinSqSigma * saturate(cosOmega);                 // x*Clip[y,0,1]
-        float cotSqSigma = rcp(sinSqSigma) - 1;
-        float tanSqSigma = rcp(cotSqSigma);
-        float sinSqOmega = 1 - cosSqOmega;
+        float cotSqSigma = rcp(sinSqSigma) - 1;                 // 1/x-1
+        float tanSqSigma = rcp(cotSqSigma);                     // (1-x)/x
+        float sinSqOmega = 1 - cosSqOmega;                      // 1-y^2
-        float x = sinSqOmega * tanSqSigma;
-        float y = -cosOmega * rsqrt(x);
+        float w = sinSqOmega * tanSqSigma;                      // (1-y^2)*((1-x)/x)
+        float x = -cosOmega * rsqrt(w);                         // -y*Sqrt[(1/x-1)/(1-y^2)]
+        float y = sqrt(sinSqOmega * tanSqSigma - cosSqOmega);   // Sqrt[(1-y^2)*((1-x)/x)-y^2]
+        float z = y * cotSqSigma;                               // Sqrt[(1-y^2)*((1-x)/x)-y^2]*(1/x-1)
-        float a = sqrt(sinSqOmega * tanSqSigma - cosSqOmega);
-        float b = a * cotSqSigma;
+        float a = cosOmega * acos(x) - z;                       // y*ArcCos[-y*Sqrt[(1/x-1)/(1-y^2)]]-Sqrt[(1-y^2)*((1-x)/x)-y^2]*(1/x-1)
+        float b = atan(y);                                      // ArcTan[Sqrt[(1-y^2)*((1-x)/x)-y^2]]
-        irradiance = INV_PI * ((cosOmega * acos(y) - b) * sinSqSigma + atan(a));
+        return saturate(INV_PI * (a * sinSqSigma + b));         // (a/Pi)*x+(b/Pi)
-    return max(irradiance, 0);
 }
 
 // Expects non-normalized vertex positions.