Great writeup!

I can't tell from your GLSL if these would have forced FMAs for a lot of the intermediate product sums. That would probably be a non-trivial effect, particularly for your large anisotropy cases.

The Heckbert paper also describes the basic theory, but you would want to supplement with some of the offline rendering work that followed it. OpenImageIO (OIIO) is pretty widely used, and has gone through several iterations of bug fixing like https://github.com/AcademySoftwareFoundation/OpenImageIO/pul...

But for your purposes, you probably just need to find all the magic epsilons and sign checks to make it match.