  For the Lagrange interpolation over a triangular domain, we propose an
efficient algorithm to rigorously evaluate the interpolation error constant
under the maximum norm by using the finite element method (FEM). In solving the
optimization problem corresponding to the interpolation error constant, the
maximum norm in the constraint condition is the most difficult part to process.
To handle this difficulty, a novel method is proposed by combining the
orthogonality of the interpolation associated to the Fujino--Morley FEM space
and the convex-hull property of the Bernstein representation of functions in
the FEM space. Numerical results for the lower and upper bounds of the
interpolation error constant for triangles of various types are presented to
verify the efficiency of the proposed method.
