I analyse the image, finding the max C (of LCh) for each out-of-gamut pixel. I then take the max. of the C_actual / C_at_gamut_boundary, and apply a curve (Tone Curve Math Question - #6 by Thanatomanic) that gradually scales the C values so they all get inside the gamut. For L, I use the same kind of curve, after determining the max brightness; the shoulder is chosen such that the brightest pixel’s L value ends up just under 1.
The code is here, but it’s not written in the style used in image processing. 
I tried to play it (type)-safe, didn’t care about performance. You’ll probably find it waaaay over-engineered (as it is).