The linked source mention "change its fusion mode to < Multiply >", so the operation to do is not a simple average of input hillshades (for this, see also How to average gdal_hillshades?). It's something else. Yet, let's create the 3 different-sunlight-directions hillshades :
gdaldem hillshade input.tif hillshades_A.tmp.tif -s 111120 -z 5 -az 315 -alt 60 -compute_edges
gdaldem hillshade input.tif hillshades_B.tmp.tif -s 111120 -z 5 -az 355 -alt 60 -compute_edges
gdaldem hillshade input.tif hillshades_C.tmp.tif -s 111120 -z 5 -az 275 -alt 60 -compute_edges
Keep lowest value of A,B,C
The first algorithm I though about is to filter and keep the darkest pixels, aka the pixels with lower values among input A,B,C.A boolean can do that :
gdal_calc.py -A hillshades_A.tmp.tif -B hillshades_B.tmp.tif -C hillshades_C.tmp.tif --outfile=./hillshades_xl.tmp.tif \
--calc="(A*(A<=B)*(A<=C)+ B*(B<A)*(B<=C)+ C*(C<A)*(C<B))"
The area dominated by shadows now make up more than the opposite side of one central lignt, it been increase by 40⁰ on each side. Not as the link provided, this current algorithm seems to lost too much the enlighten area.
Angle of 315±30⁰ (smaller angle variation) rather that current 315±40⁰ would do nicer.
The diagram below is a the basis of the equation. It shows light sources A, B, C, and the Boolean comparisons for pixels values A, B, C in each part. Equality lines need a special attention to be included into the Boolean. The median lines have value 221 for the perpendicular light source. Think about areas of influence, the closest source of light is the main influence, and the farest the weakest influence.
Keep extremes values for A,B,C
An other Boolean algorithm could be to keep the most extreme values, both the darkest and the whitest pixels. The following diagram helps to thing about the Boolean formula. For each sixth of the circle, it identifies the value to keep from A,B,C, and a Boolean to select the triangular area plus the clockwise equality line, and only that. It gives (from top and clockwise) :
--calc="A*(A>B)*(A>=C) + C*(C>A)*(A>=221) + B*(C>A)*(C>=221) + A*(A<B)*(B<=221) + C*(B>A)*(A<=221) + B*(B>A)*(C<=221)"
If the angles variations is not too important, it could give good results.
Other Booleans
You may create more complex Boolean to cover the whole circles using any combination the border segments. It stays important that only one value out of A,B,C is kept for one segment.
Multiply
I made several failled attempts to multiply pixels values without a proven formula nor final success.
@Radouxju pointed out that (a*b*c)^(1/3) (GEOMETRIC mean) instead of the ARITHMETIC mean (a*b*c)/(255*255) may work. Geometric mean is lower or equal to arithmetic mean, which accentuates the darkness of shaded areas. I haven't tested it yet.
