Going from hillshade to percent shaded using gdaldem?

I assumed hillshades output percent shaded (scaled to 255) at the given azimuth and altitude but this is incorrect. A flat raster at 80 degrees altitude has a maximum value of 251 instead of 255. How can I go from a hillshade to percent shaded and can I simply divide by the cosine of the 10 degrees in this case to get the correct value or is it more complicated?

Here is some R code with gdaldem that may help you see what I'm seeing:

library(raster)
r <- raster(ncol=100,nrow=100)
values(r) <- rep(1,ncell(r))
writeRaster(r,"flatelevation.tif")
system("gdaldem hillshade -az 180 -alt 80 flatelevation.tif hillshade_flat.tif")
summary(values(hs))
summary(values(hs)/cos((90.0-80.0)*pi/180.0))

I will try to parse the source code but it's been ten years since I've used C++.

Can someone with general knowledge about hillshades point me the right direction?

I thought I could get a "percent of the tile that is shaded" from running a hillshade. Upon further reading a hillshade outputs the illuminance of each tile. So, the cosine is only part of the equation, I would need to pull in the slope and aspect as well for additional trigonometric readjustments, and maybe more.

Turns out there are tools that are functionally equivalent for that I need called "shadow masks" that only output 1 or 0, indicating whether the tile is shaded.

This nail biting article referenced in gdal's source code helped enlighten me, in addition to this question on cloud shading, this question on solar maps, and this question on r.sun where a contributor indicates that r.sunmask is slow. A couple available tools: