[NAThompson]

Adds tools::multivariate_newton and tools::multivariate_halley in boost/math/tools/multivariate_roots.hpp for systems F(x) = 0. f(x) returns (F, J), or (F, J, d2) for Halley, where d2(v) = F″(x)[v, v], so the n³ tensor is never formed. Halley is Cuyt's Padé form: with a = −J⁻¹F and b = J⁻¹F″[a, a], xᵢ += aᵢ² / (aᵢ + bᵢ/2), with one LU factorization serving both solves (Cuyt & Van der Cruyssen 1983; Cuyt & Rall, ACM TOMS 1985). Both use a backtracking line search on ½‖F‖²; a singular J or a failed line search raises an evaluation error. Also bound in the nanobind layer; example/multivariate_roots_ellipse.cpp gets J and F″[v, v] from forward-mode autodiff.
Closes #1616

Newton doubles the digits per step and Halley triples them (cpp_bin_float_50):

Basins of z³ = 1 as a real 2-D system. Componentwise Halley is not complex Halley: undamped, its basin boundaries are wilder than Newton's; the line search tames both, and Halley needs 28% fewer iterations.

Every pixel solves a 3×3 ray–torus system:

google/benchmark (reporting/performance/multivariate_roots_performance.cpp, Apple clang -O2 arm64):

system	Newton	Halley
Broyden tridiagonal, n = 10	3.57 µs	3.05 µs
Broyden tridiagonal, n = 100	629 µs	466 µs
Rosenbrock gradient from (−1.2, 1)	63.3 µs	14.4 µs

Limitations: the LU is dense; the componentwise quotient is not affine-invariant; a line search along the Newton direction creeps through curved valleys (253 Newton iterations on the Rosenbrock gradient).

🤖 Generated with Claude Code