Add integrate_1d_gauss_kronord

[avehtari]

integrate_1d uses Boost double exponential quadrature https://live.boost.org/doc/libs/1_74_0/libs/math/doc/html/math_toolkit/double_exponential.html, which is good when there are singularities at the end points, but the problem is that "These routines are very aggressive about approaching the endpoint singularities." which means they are less efficient when there are no singularities at the end points.

A typical smooth function with no end point singularities is normal distribution times log concave likelihood, where the tails are also decreasing fast so that usually there is a negligible contribution near bounds, and double exponential quadrature is wasting computation. For such smooth functions it would be better to use (adaptive) Gauss-Kronord which is also available in Boost https://live.boost.org/doc/libs/1_74_0/libs/math/doc/html/math_toolkit/gauss_kronrod.html. From Boost doc "they should be chosen for the integration of smooth functions with no end point singularities. For difficult functions, or those with end point singularities, please refer to the double-exponential integration schemes"

The function could be named integrate_1d_gauss_kronord

As default this could use adaptive Gauss-Kronord quadrature, and for the best speed it would be nice to allow non-adaptive Gauss-Kronord with the user chosen number of evaluation points provided by one of the arguments.