# How is the smoothing spline penalty computed in practice?

I'm digging into smoothing splines and finding good resources, but no one talks about how to actually calculate the penalty $$\int \hat{f}^{"}(x)^2 dx$$ in the standard smoothing spline loss:

Since f is a linear function of B-splines (bs in R's stats library), I'd wonder if there's some built-in differentiator, but it doesn't look like it, and smooth.spline quickly hits Fortran code.

My main question: What is the "smart way" to do this computationally?

Optional: Are there less efficient methods in R that would get the job done with very few lines of code and use properties of the B-splines via bs?

Argument for the uniqueness of this question. The answer to this question could have been submitted in How to measure smoothness of a time series in R?, but it wasn't (and wasn't necessary to answer the question). In How to measure smoothness/roughness of a time series (duplicate), @David Waterworth suggests using standard derivative and integral approximations, but I would consider that a last resort.