From 16c7a7f741930b41f0badfcdca1dece3c5ae4f60 Mon Sep 17 00:00:00 2001 From: Jacob Hass Date: Sun, 6 Sep 2026 18:41:39 -0700 Subject: [PATCH 1/4] First pass at hermite polynomials --- include/boost/math/quadrature/hermite.hpp | 572 ++++++++++++++++++++++ 1 file changed, 572 insertions(+) create mode 100644 include/boost/math/quadrature/hermite.hpp diff --git a/include/boost/math/quadrature/hermite.hpp b/include/boost/math/quadrature/hermite.hpp new file mode 100644 index 0000000000..b08c2ecc11 --- /dev/null +++ b/include/boost/math/quadrature/hermite.hpp @@ -0,0 +1,572 @@ +// Copyright Jacob Hass 2026. +// Use, modification and distribution are subject to the +// Boost Software License, Version 1.0. (See accompanying file +// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) + +#ifndef BOOST_MATH_QUADRATURE_HERMITE_HPP +#define BOOST_MATH_QUADRATURE_HERMITE_HPP + +#include +#include + +#if __has_include() + #include + #include + #include + #include + #include + + #define EIGEN_SUPPORT +#endif + +namespace boost { namespace math{ namespace quadrature{ namespace detail { + + +#if defined(EIGEN_SUPPORT) && !defined(BOOST_MATH_GAUSS_NO_COMPUTE_ON_DEMAND) + +template +Real a(unsigned& n) +{ + return Real(0); +} + +template +Real b(unsigned& n) +{ + return Real(n) / Real(2); +} + +template +Real factorial(unsigned n) +{ + Real i = 1; + Real factorial = 1; + + while (i <= n) + { + factorial *= i; + i++; + } + return factorial; +} + +template +class hermite_detail +{ +public: + static std::pair, std::vector > calculate_values() + { + namespace mp = boost::multiprecision; + + Eigen::Matrix P = Eigen::Matrix::Zero(N, N); + for (unsigned n = 0; n < N; ++n) + { + P(n, n) = a(n); + + if (n > 0) + { + P(n, n-1) = sqrt(b(n)); + P(n-1, n) = sqrt(b(n)); + } + } + + Eigen::SelfAdjointEigenSolver > solver(P); + Eigen::Vector roots = solver.eigenvalues(); + + std::vector abscissa_vals(N); + std::vector weight_vals(N); + + for (unsigned i=0; i < roots.size(); i++) + { + abscissa_vals[i] = roots(i); + Real hermite_val = boost::math::hermite(N-1, roots(i)); + Real weight = mp::pow(Real(2), N-1) * factorial(N) * boost::math::constants::root_pi() / mp::pow(Real(N), 2) / mp::pow(hermite_val, 2.0); + weight_vals[i] = weight; + } + + // If N is odd, middle abscissa is 0 + if (N % 2 != 0) + { + abscissa_vals[N / 2] = Real(0); + } + + return std::make_pair(abscissa_vals, weight_vals); + } + + static const std::vector& weights() + { + static std::pair, std::vector > data = calculate_values(); + return data.second; + } + + static const std::vector& abscissa() + { + static std::pair, std::vector > data = calculate_values(); + return data.first; + } +}; + +#else + +template +class hermite_detail; + +#endif + +template +struct hermite_constant_category +{ + static const unsigned value = + (std::numeric_limits::is_specialized == 0) ? 999 : + (std::numeric_limits::radix == 2) ? + ( +#ifdef BOOST_HAS_FLOAT128 + (std::numeric_limits::digits <= 113) && std::is_constructible::value ? 0 : +#else + (std::numeric_limits::digits <= std::numeric_limits::digits) && std::is_constructible::value ? 0 : +#endif + (std::numeric_limits::digits10 <= 110) && std::is_constructible::value ? 4 : 999 + ) : (std::numeric_limits::digits10 <= 110) && std::is_constructible::value ? 4 : 999; + + using storage_type = + std::conditional_t<(std::numeric_limits::is_specialized == 0), T, + std::conditional_t<(std::numeric_limits::radix == 2), + std::conditional_t< ((std::numeric_limits::digits <= std::numeric_limits::digits) && std::is_constructible::value), + float, + std::conditional_t<((std::numeric_limits::digits <= std::numeric_limits::digits) && std::is_constructible::value), + double, + std::conditional_t<((std::numeric_limits::digits <= std::numeric_limits::digits) && std::is_constructible::value), + long double, +#ifdef BOOST_HAS_FLOAT128 + std::conditional_t<((std::numeric_limits::digits <= 113) && std::is_constructible::value), + __float128, + T + > + > +#else + T + > +#endif + > + >, T + > + >; +}; + +#ifndef BOOST_HAS_FLOAT128 +template +class hermite_detail +{ + using storage_type = typename hermite_constant_category::storage_type; + public: + static std::array const & abscissa() + { + static std::array data = { + static_cast(-2.65196135683523349244708200651661611e+00L), + static_cast(-1.67355162876747144503180139830359482e+00L), + static_cast(-8.16287882858964663038710959027145817e-01L), + static_cast(0.00000000000000000000000000000000000e+00L), + static_cast(8.16287882858964663038710959027145817e-01L), + static_cast(1.67355162876747144503180139830359482e+00L), + static_cast(2.65196135683523349244708200651661611e+00L), +}; + return data; + } + static std::array const & weights() + { + static std::array data = { + static_cast(9.71781245099519154149424255938959644e-04L), + static_cast(5.45155828191270305921785688416951260e-02L), + static_cast(4.25607252610127800520317466666391035e-01L), + static_cast(8.10264617556807326764876563813094941e-01L), + static_cast(4.25607252610127800520317466666391035e-01L), + static_cast(5.45155828191270305921785688416951260e-02L), + static_cast(9.71781245099519154149424255938959644e-04L), + }; + return data; + } + }; + +#else +template +class hermite_detail +{ + using storage_type = typename hermite_constant_category::storage_type; + public: + static std::array const & abscissa() + { + static std::array data = { + static_cast(-2.65196135683523349244708200651661611e+00Q), + static_cast(-1.67355162876747144503180139830359482e+00Q), + static_cast(-8.16287882858964663038710959027145817e-01Q), + static_cast(0.00000000000000000000000000000000000e+00Q), + static_cast(8.16287882858964663038710959027145817e-01Q), + static_cast(1.67355162876747144503180139830359482e+00Q), + static_cast(2.65196135683523349244708200651661611e+00Q), +}; + return data; + } + static std::array const & weights() + { + static std::array data = { + static_cast(9.71781245099519154149424255938959644e-04Q), + static_cast(5.45155828191270305921785688416951260e-02Q), + static_cast(4.25607252610127800520317466666391035e-01Q), + static_cast(8.10264617556807326764876563813094941e-01Q), + static_cast(4.25607252610127800520317466666391035e-01Q), + static_cast(5.45155828191270305921785688416951260e-02Q), + static_cast(9.71781245099519154149424255938959644e-04Q), + }; + return data; + } + }; + +#endif +template +class hermite_detail +{ + using storage_type = typename hermite_constant_category::storage_type; + public: + static std::array const & abscissa() + { + static std::array data = { + BOOST_MATH_HUGE_CONSTANT(T, 0, -2.6519613568352334924470820065166161144381584786255294172031003071471600949016031668052906818697878641942881494964276e+00), + BOOST_MATH_HUGE_CONSTANT(T, 0, -1.6735516287674714450318013983035948191078100577354089269242175099937138333690347986612540761793805055242734953861690e+00), + BOOST_MATH_HUGE_CONSTANT(T, 0, -8.1628788285896466303871095902714581674288940037863615684472203343594907048766511668519794976704116666704491757953733e-01), + BOOST_MATH_HUGE_CONSTANT(T, 0, 0.0000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000e+00), + BOOST_MATH_HUGE_CONSTANT(T, 0, 8.1628788285896466303871095902714581674288940037863615684472203343594907048766511668519794976704116666704491757953733e-01), + BOOST_MATH_HUGE_CONSTANT(T, 0, 1.6735516287674714450318013983035948191078100577354089269242175099937138333690347986612540761793805055242734953861690e+00), + BOOST_MATH_HUGE_CONSTANT(T, 0, 2.6519613568352334924470820065166161144381584786255294172031003071471600949016031668052906818697878641942881494964276e+00), +}; + return data; + } + static std::array const & weights() + { + static std::array data = { + BOOST_MATH_HUGE_CONSTANT(T, 0, 9.7178124509951915414942425593895964444240121701420164980001433798334142698580146031198718271051220413580750089473744e-04), + BOOST_MATH_HUGE_CONSTANT(T, 0, 5.4515582819127030592178568841695125960899915859346492986586846631369448562711195465546454603039043598893480308577392e-02), + BOOST_MATH_HUGE_CONSTANT(T, 0, 4.2560725261012780052031746666639103543970682101400152587878953913358027297072602373305673647185542480994656431711885e-01), + BOOST_MATH_HUGE_CONSTANT(T, 0, 8.1026461755680732676487656381309494070745117994166268718345498964704515867018614005712030022333470419027511937110060e-01), + BOOST_MATH_HUGE_CONSTANT(T, 0, 4.2560725261012780052031746666639103543970682101400152587878953913358027297072602373305673647185542480994656431711885e-01), + BOOST_MATH_HUGE_CONSTANT(T, 0, 5.4515582819127030592178568841695125960899915859346492986586846631369448562711195465546454603039043598893480308577392e-02), + BOOST_MATH_HUGE_CONSTANT(T, 0, 9.7178124509951915414942425593895964444240121701420164980001433798334142698580146031198718271051220413580750089473744e-04), + }; + return data; + } + }; + +#ifndef BOOST_HAS_FLOAT128 +template +class hermite_detail +{ + using storage_type = typename hermite_constant_category::storage_type; + public: + static std::array const & abscissa() + { + static std::array data = { + static_cast(-3.43615911883773760332672549431912138e+00L), + static_cast(-2.53273167423278979640896079775479348e+00L), + static_cast(-1.75668364929988177345140122010615676e+00L), + static_cast(-1.03661082978951365417749191675920902e+00L), + static_cast(-3.42901327223704608789165025557258031e-01L), + static_cast(3.42901327223704608789165025557258031e-01L), + static_cast(1.03661082978951365417749191675920902e+00L), + static_cast(1.75668364929988177345140122010615676e+00L), + static_cast(2.53273167423278979640896079775479348e+00L), + static_cast(3.43615911883773760332672549431912138e+00L), +}; + return data; + } + static std::array const & weights() + { + static std::array data = { + static_cast(7.64043285523262062915936785959522211e-06L), + static_cast(1.34364574678123269220156558584591387e-03L), + static_cast(3.38743944554810631361647312775859737e-02L), + static_cast(2.40138611082314686416523295005861395e-01L), + static_cast(6.10862633735325798783564990433419713e-01L), + static_cast(6.10862633735325798783564990433419713e-01L), + static_cast(2.40138611082314686416523295005861395e-01L), + static_cast(3.38743944554810631361647312775859737e-02L), + static_cast(1.34364574678123269220156558584591387e-03L), + static_cast(7.64043285523262062915936785959522211e-06L), + }; + return data; + } + }; + +#else +template +class hermite_detail +{ + using storage_type = typename hermite_constant_category::storage_type; + public: + static std::array const & abscissa() + { + static std::array data = { + static_cast(-3.43615911883773760332672549431912138e+00Q), + static_cast(-2.53273167423278979640896079775479348e+00Q), + static_cast(-1.75668364929988177345140122010615676e+00Q), + static_cast(-1.03661082978951365417749191675920902e+00Q), + static_cast(-3.42901327223704608789165025557258031e-01Q), + static_cast(3.42901327223704608789165025557258031e-01Q), + static_cast(1.03661082978951365417749191675920902e+00Q), + static_cast(1.75668364929988177345140122010615676e+00Q), + static_cast(2.53273167423278979640896079775479348e+00Q), + static_cast(3.43615911883773760332672549431912138e+00Q), +}; + return data; + } + static std::array const & weights() + { + static std::array data = { + static_cast(7.64043285523262062915936785959522211e-06Q), + static_cast(1.34364574678123269220156558584591387e-03Q), + static_cast(3.38743944554810631361647312775859737e-02Q), + static_cast(2.40138611082314686416523295005861395e-01Q), + static_cast(6.10862633735325798783564990433419713e-01Q), + static_cast(6.10862633735325798783564990433419713e-01Q), + static_cast(2.40138611082314686416523295005861395e-01Q), + static_cast(3.38743944554810631361647312775859737e-02Q), + static_cast(1.34364574678123269220156558584591387e-03Q), + static_cast(7.64043285523262062915936785959522211e-06Q), + }; + return data; + } + }; + +#endif +template +class hermite_detail +{ + using storage_type = typename hermite_constant_category::storage_type; + public: + static std::array const & abscissa() + { + static std::array data = { + BOOST_MATH_HUGE_CONSTANT(T, 0, -3.4361591188377376033267254943191213848406783093901772906871762944002361145649511153468627231618843292611320162500620e+00), + BOOST_MATH_HUGE_CONSTANT(T, 0, -2.5327316742327897964089607977547934803078465081567249459332874134930859380473865340227292325941771646240824732515703e+00), + BOOST_MATH_HUGE_CONSTANT(T, 0, -1.7566836492998817734514012201061567632954744937388471000872849560961705796098190771278266711469838314640623161517203e+00), + BOOST_MATH_HUGE_CONSTANT(T, 0, -1.0366108297895136541774919167592090162982561106572157244372326874036404187243878992980306028243900733070850998719429e+00), + BOOST_MATH_HUGE_CONSTANT(T, 0, -3.4290132722370460878916502555725803120830265867773062820577103907292018192312101560713520305443248063627958675235231e-01), + BOOST_MATH_HUGE_CONSTANT(T, 0, 3.4290132722370460878916502555725803120830265867773062820577103907292018192312101560713520305443248063627958675235231e-01), + BOOST_MATH_HUGE_CONSTANT(T, 0, 1.0366108297895136541774919167592090162982561106572157244372326874036404187243878992980306028243900733070850998719429e+00), + BOOST_MATH_HUGE_CONSTANT(T, 0, 1.7566836492998817734514012201061567632954744937388471000872849560961705796098190771278266711469838314640623161517203e+00), + BOOST_MATH_HUGE_CONSTANT(T, 0, 2.5327316742327897964089607977547934803078465081567249459332874134930859380473865340227292325941771646240824732515703e+00), + BOOST_MATH_HUGE_CONSTANT(T, 0, 3.4361591188377376033267254943191213848406783093901772906871762944002361145649511153468627231618843292611320162500620e+00), +}; + return data; + } + static std::array const & weights() + { + static std::array data = { + BOOST_MATH_HUGE_CONSTANT(T, 0, 7.6404328552326206291593678595952221082891147358323631952885551280911633480005783690893185993195845416600950144409630e-06), + BOOST_MATH_HUGE_CONSTANT(T, 0, 1.3436457467812326922015655858459138698692732786257703646162498020762483818413155201782152234375992775169973716992552e-03), + BOOST_MATH_HUGE_CONSTANT(T, 0, 3.3874394455481063136164731277585973698161942481867016739715157699771381277813573402914401584449869233310982466662099e-02), + BOOST_MATH_HUGE_CONSTANT(T, 0, 2.4013861108231468641652329500586139537004249390192739899066056147994315237169314350600241149865696523148172530499056e-01), + BOOST_MATH_HUGE_CONSTANT(T, 0, 6.1086263373532579878356499043341971323859272928403754564871663738953676910082005768001121074412857938126204657377493e-01), + BOOST_MATH_HUGE_CONSTANT(T, 0, 6.1086263373532579878356499043341971323859272928403754564871663738953676910082005768001121074412857938126204657377493e-01), + BOOST_MATH_HUGE_CONSTANT(T, 0, 2.4013861108231468641652329500586139537004249390192739899066056147994315237169314350600241149865696523148172530499056e-01), + BOOST_MATH_HUGE_CONSTANT(T, 0, 3.3874394455481063136164731277585973698161942481867016739715157699771381277813573402914401584449869233310982466662099e-02), + BOOST_MATH_HUGE_CONSTANT(T, 0, 1.3436457467812326922015655858459138698692732786257703646162498020762483818413155201782152234375992775169973716992552e-03), + BOOST_MATH_HUGE_CONSTANT(T, 0, 7.6404328552326206291593678595952221082891147358323631952885551280911633480005783690893185993195845416600950144409630e-06), + }; + return data; + } + }; + +#ifndef BOOST_HAS_FLOAT128 +template +class hermite_detail +{ + using storage_type = typename hermite_constant_category::storage_type; + public: + static std::array const & abscissa() + { + static std::array data = { + static_cast(-4.49999070730939155366438053053482422e+00L), + static_cast(-3.66995037340445253472922383311568148e+00L), + static_cast(-2.96716692790560324848896036354980632e+00L), + static_cast(-2.32573248617385774545404479448534701e+00L), + static_cast(-1.71999257518648893241583152515256151e+00L), + static_cast(-1.13611558521092066631913490555610213e+00L), + static_cast(-5.65069583255575748526020337198188867e-01L), + static_cast(0.00000000000000000000000000000000000e+00L), + static_cast(5.65069583255575748526020337198188867e-01L), + static_cast(1.13611558521092066631913490555610213e+00L), + static_cast(1.71999257518648893241583152515256151e+00L), + static_cast(2.32573248617385774545404479448534701e+00L), + static_cast(2.96716692790560324848896036354980632e+00L), + static_cast(3.66995037340445253472922383311568148e+00L), + static_cast(4.49999070730939155366438053053482422e+00L), +}; + return data; + } + static std::array const & weights() + { + static std::array data = { + static_cast(1.52247580425351702016062666964826618e-09L), + static_cast(1.05911554771106663577520791055016329e-06L), + static_cast(1.00004441232499868127296736176978509e-04L), + static_cast(2.77806884291277589607887049229213490e-03L), + static_cast(3.07800338725460822286814158757801456e-02L), + static_cast(1.58488915795935746883839384959994027e-01L), + static_cast(4.12028687498898627025891079567810319e-01L), + static_cast(5.64100308726417532852625797339963533e-01L), + static_cast(4.12028687498898627025891079567810319e-01L), + static_cast(1.58488915795935746883839384959994027e-01L), + static_cast(3.07800338725460822286814158757801456e-02L), + static_cast(2.77806884291277589607887049229213490e-03L), + static_cast(1.00004441232499868127296736176978509e-04L), + static_cast(1.05911554771106663577520791055016329e-06L), + static_cast(1.52247580425351702016062666964826618e-09L), + }; + return data; + } + }; + +#else +template +class hermite_detail +{ + using storage_type = typename hermite_constant_category::storage_type; + public: + static std::array const & abscissa() + { + static std::array data = { + static_cast(-4.49999070730939155366438053053482422e+00Q), + static_cast(-3.66995037340445253472922383311568148e+00Q), + static_cast(-2.96716692790560324848896036354980632e+00Q), + static_cast(-2.32573248617385774545404479448534701e+00Q), + static_cast(-1.71999257518648893241583152515256151e+00Q), + static_cast(-1.13611558521092066631913490555610213e+00Q), + static_cast(-5.65069583255575748526020337198188867e-01Q), + static_cast(0.00000000000000000000000000000000000e+00Q), + static_cast(5.65069583255575748526020337198188867e-01Q), + static_cast(1.13611558521092066631913490555610213e+00Q), + static_cast(1.71999257518648893241583152515256151e+00Q), + static_cast(2.32573248617385774545404479448534701e+00Q), + static_cast(2.96716692790560324848896036354980632e+00Q), + static_cast(3.66995037340445253472922383311568148e+00Q), + static_cast(4.49999070730939155366438053053482422e+00Q), +}; + return data; + } + static std::array const & weights() + { + static std::array data = { + static_cast(1.52247580425351702016062666964826618e-09Q), + static_cast(1.05911554771106663577520791055016329e-06Q), + static_cast(1.00004441232499868127296736176978509e-04Q), + static_cast(2.77806884291277589607887049229213490e-03Q), + static_cast(3.07800338725460822286814158757801456e-02Q), + static_cast(1.58488915795935746883839384959994027e-01Q), + static_cast(4.12028687498898627025891079567810319e-01Q), + static_cast(5.64100308726417532852625797339963533e-01Q), + static_cast(4.12028687498898627025891079567810319e-01Q), + static_cast(1.58488915795935746883839384959994027e-01Q), + static_cast(3.07800338725460822286814158757801456e-02Q), + static_cast(2.77806884291277589607887049229213490e-03Q), + static_cast(1.00004441232499868127296736176978509e-04Q), + static_cast(1.05911554771106663577520791055016329e-06Q), + static_cast(1.52247580425351702016062666964826618e-09Q), + }; + return data; + } + }; + +#endif +template +class hermite_detail +{ + using storage_type = typename hermite_constant_category::storage_type; + public: + static std::array const & abscissa() + { + static std::array data = { + BOOST_MATH_HUGE_CONSTANT(T, 0, -4.4999907073093915536643805305348242199319830534949905938271742846019366228525630586264389556919589962713446071341276e+00), + BOOST_MATH_HUGE_CONSTANT(T, 0, -3.6699503734044525347292238331156814846958556260108268200827250197712929612325341525137596165272567610075980577793224e+00), + BOOST_MATH_HUGE_CONSTANT(T, 0, -2.9671669279056032484889603635498063155705164076381381507212727522118956916046294049934561649283587290854852693064039e+00), + BOOST_MATH_HUGE_CONSTANT(T, 0, -2.3257324861738577454540447944853470123070950095634828957721625948740238048583608727061841554643103134570736658336089e+00), + BOOST_MATH_HUGE_CONSTANT(T, 0, -1.7199925751864889324158315251525615148429195630613888498938216300095132434761832428439718888176122448729491067241320e+00), + BOOST_MATH_HUGE_CONSTANT(T, 0, -1.1361155852109206663191349055561021252835292947846650850648392599800715967101004409631966601830482707422651868052475e+00), + BOOST_MATH_HUGE_CONSTANT(T, 0, -5.6506958325557574852602033719818886681552078800505788680930739682745058470915375159973526283034610292989738379345761e-01), + BOOST_MATH_HUGE_CONSTANT(T, 0, 0.0000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000e+00), + BOOST_MATH_HUGE_CONSTANT(T, 0, 5.6506958325557574852602033719818886681552078800505788680930739682745058470915375159973526283034610292989738379345761e-01), + BOOST_MATH_HUGE_CONSTANT(T, 0, 1.1361155852109206663191349055561021252835292947846650850648392599800715967101004409631966601830482707422651868052475e+00), + BOOST_MATH_HUGE_CONSTANT(T, 0, 1.7199925751864889324158315251525615148429195630613888498938216300095132434761832428439718888176122448729491067241320e+00), + BOOST_MATH_HUGE_CONSTANT(T, 0, 2.3257324861738577454540447944853470123070950095634828957721625948740238048583608727061841554643103134570736658336089e+00), + BOOST_MATH_HUGE_CONSTANT(T, 0, 2.9671669279056032484889603635498063155705164076381381507212727522118956916046294049934561649283587290854852693064039e+00), + BOOST_MATH_HUGE_CONSTANT(T, 0, 3.6699503734044525347292238331156814846958556260108268200827250197712929612325341525137596165272567610075980577793224e+00), + BOOST_MATH_HUGE_CONSTANT(T, 0, 4.4999907073093915536643805305348242199319830534949905938271742846019366228525630586264389556919589962713446071341276e+00), +}; + return data; + } + static std::array const & weights() + { + static std::array data = { + BOOST_MATH_HUGE_CONSTANT(T, 0, 1.5224758042535170201606266696482661827962160316014537930014349599535491582697947140375830838937903702636262186650310e-09), + BOOST_MATH_HUGE_CONSTANT(T, 0, 1.0591155477110666357752079105501632872963524018315499158761728455091303265515893684325618169279641299843625141139696e-06), + BOOST_MATH_HUGE_CONSTANT(T, 0, 1.0000444123249986812729673617697850851699387810876387789234107161752152460966503439378917351614005951603805751399265e-04), + BOOST_MATH_HUGE_CONSTANT(T, 0, 2.7780688429127758960788704922921349030124251134822666685717061210098979501372205289327949738880183107056484267658214e-03), + BOOST_MATH_HUGE_CONSTANT(T, 0, 3.0780033872546082228681415875780145569645078095357831670126204320943426296907483543543105075772562113024847790201415e-02), + BOOST_MATH_HUGE_CONSTANT(T, 0, 1.5848891579593574688383938495999402721002703620834420484147417162556491543397074116784210832986285797844312792251258e-01), + BOOST_MATH_HUGE_CONSTANT(T, 0, 4.1202868749889862702589107956781031893375232741692141877844273550979447720862335855236503495494156330750630603931556e-01), + BOOST_MATH_HUGE_CONSTANT(T, 0, 5.6410030872641753285262579733996353292453477640072243023805413420648986959478382400247069852378235816199417789938317e-01), + BOOST_MATH_HUGE_CONSTANT(T, 0, 4.1202868749889862702589107956781031893375232741692141877844273550979447720862335855236503495494156330750630603931556e-01), + BOOST_MATH_HUGE_CONSTANT(T, 0, 1.5848891579593574688383938495999402721002703620834420484147417162556491543397074116784210832986285797844312792251258e-01), + BOOST_MATH_HUGE_CONSTANT(T, 0, 3.0780033872546082228681415875780145569645078095357831670126204320943426296907483543543105075772562113024847790201415e-02), + BOOST_MATH_HUGE_CONSTANT(T, 0, 2.7780688429127758960788704922921349030124251134822666685717061210098979501372205289327949738880183107056484267658214e-03), + BOOST_MATH_HUGE_CONSTANT(T, 0, 1.0000444123249986812729673617697850851699387810876387789234107161752152460966503439378917351614005951603805751399265e-04), + BOOST_MATH_HUGE_CONSTANT(T, 0, 1.0591155477110666357752079105501632872963524018315499158761728455091303265515893684325618169279641299843625141139696e-06), + BOOST_MATH_HUGE_CONSTANT(T, 0, 1.5224758042535170201606266696482661827962160316014537930014349599535491582697947140375830838937903702636262186650310e-09), + }; + return data; + } + }; + +} // namespace detail + +template > +class hermite : public detail::hermite_detail::value> +{ + typedef detail::hermite_detail::value> base; +public: + + template + static auto integrate(F f, Real* pL1 = nullptr)->decltype(std::declval()(std::declval())) + { + // In many math texts, K represents the field of real or complex numbers. + // Too bad we can't put blackboard bold into C++ source! + typedef decltype(f(Real(0))) K; + static_assert(!std::is_integral::value, + "The return type cannot be integral, it must be either a real or complex floating point type."); + using std::abs; + unsigned non_zero_start = 1; + K result = Real(0); + if (N & 1) { + result = f(Real(0)) * static_cast(base::weights()[0]); + } + else { + result = 0; + non_zero_start = 0; + } + Real L1 = abs(result); + Real weight_total; + for (unsigned i = non_zero_start; i < base::abscissa().size(); ++i) + { + K fp = f(static_cast(base::abscissa()[i])); + result += fp * static_cast(base::weights()[i]); + L1 += abs(fp) * static_cast(base::weights()[i]); + } + if (pL1) + *pL1 = L1; + return result; + } +}; + +} // namespace quadrature +} // namespace math +} // namespace boost + +#endif // BOOST_MATH_QUADRATURE_HERMITE_HPP \ No newline at end of file From d3a5ddfaf4087cb9acae6a6c97fe7cfcc8feacd4 Mon Sep 17 00:00:00 2001 From: Jacob Hass Date: Mon, 7 Sep 2026 11:54:21 -0700 Subject: [PATCH 2/4] Excluded negative roots. Some numbers moved - need to check --- include/boost/math/quadrature/hermite.hpp | 297 ++++++++-------------- 1 file changed, 108 insertions(+), 189 deletions(-) diff --git a/include/boost/math/quadrature/hermite.hpp b/include/boost/math/quadrature/hermite.hpp index b08c2ecc11..7c7d13f3f2 100644 --- a/include/boost/math/quadrature/hermite.hpp +++ b/include/boost/math/quadrature/hermite.hpp @@ -6,8 +6,15 @@ #ifndef BOOST_MATH_QUADRATURE_HERMITE_HPP #define BOOST_MATH_QUADRATURE_HERMITE_HPP +#include +#include +#include +#include +#include +#include #include #include +#include #if __has_include() #include @@ -15,6 +22,7 @@ #include #include #include + #include #define EIGEN_SUPPORT #endif @@ -24,18 +32,6 @@ namespace boost { namespace math{ namespace quadrature{ namespace detail { #if defined(EIGEN_SUPPORT) && !defined(BOOST_MATH_GAUSS_NO_COMPUTE_ON_DEMAND) -template -Real a(unsigned& n) -{ - return Real(0); -} - -template -Real b(unsigned& n) -{ - return Real(n) / Real(2); -} - template Real factorial(unsigned n) { @@ -61,33 +57,45 @@ class hermite_detail Eigen::Matrix P = Eigen::Matrix::Zero(N, N); for (unsigned n = 0; n < N; ++n) { - P(n, n) = a(n); + P(n, n) = Real(0); if (n > 0) { - P(n, n-1) = sqrt(b(n)); - P(n-1, n) = sqrt(b(n)); + P(n, n-1) = sqrt(Real(n) / Real(2)); + P(n-1, n) = sqrt(Real(n) / Real(2)); } } Eigen::SelfAdjointEigenSolver > solver(P); Eigen::Vector roots = solver.eigenvalues(); - std::vector abscissa_vals(N); - std::vector weight_vals(N); - + // Filter out negative roots + std::vector pos_roots; + pos_roots.reserve(N); for (unsigned i=0; i < roots.size(); i++) { - abscissa_vals[i] = roots(i); - Real hermite_val = boost::math::hermite(N-1, roots(i)); + // Need to account for root near 0 + if (roots(i) >= -std::numeric_limits::epsilon() * 1000) + { + pos_roots.push_back(roots(i)); + } + } + + std::vector abscissa_vals(pos_roots.size()); + std::vector weight_vals(pos_roots.size()); + + for (unsigned i=0; i < pos_roots.size(); i++) + { + abscissa_vals[i] = pos_roots[i]; + Real hermite_val = boost::math::hermite(N-1, pos_roots[i]); Real weight = mp::pow(Real(2), N-1) * factorial(N) * boost::math::constants::root_pi() / mp::pow(Real(N), 2) / mp::pow(hermite_val, 2.0); weight_vals[i] = weight; } - // If N is odd, middle abscissa is 0 + // If N is odd, first abscissa is 0 if (N % 2 != 0) { - abscissa_vals[N / 2] = Real(0); + abscissa_vals[0] = Real(0); } return std::make_pair(abscissa_vals, weight_vals); @@ -155,16 +163,13 @@ struct hermite_constant_category #ifndef BOOST_HAS_FLOAT128 template -class hermite_detail +class hermite_detail { using storage_type = typename hermite_constant_category::storage_type; public: - static std::array const & abscissa() + static std::array const & abscissa() { - static std::array data = { - static_cast(-2.65196135683523349244708200651661611e+00L), - static_cast(-1.67355162876747144503180139830359482e+00L), - static_cast(-8.16287882858964663038710959027145817e-01L), + static std::array data = { static_cast(0.00000000000000000000000000000000000e+00L), static_cast(8.16287882858964663038710959027145817e-01L), static_cast(1.67355162876747144503180139830359482e+00L), @@ -172,12 +177,9 @@ class hermite_detail }; return data; } - static std::array const & weights() + static std::array const & weights() { - static std::array data = { - static_cast(9.71781245099519154149424255938959644e-04L), - static_cast(5.45155828191270305921785688416951260e-02L), - static_cast(4.25607252610127800520317466666391035e-01L), + static std::array data = { static_cast(8.10264617556807326764876563813094941e-01L), static_cast(4.25607252610127800520317466666391035e-01L), static_cast(5.45155828191270305921785688416951260e-02L), @@ -189,16 +191,13 @@ class hermite_detail #else template -class hermite_detail +class hermite_detail { using storage_type = typename hermite_constant_category::storage_type; public: - static std::array const & abscissa() + static std::array const & abscissa() { - static std::array data = { - static_cast(-2.65196135683523349244708200651661611e+00Q), - static_cast(-1.67355162876747144503180139830359482e+00Q), - static_cast(-8.16287882858964663038710959027145817e-01Q), + static std::array data = { static_cast(0.00000000000000000000000000000000000e+00Q), static_cast(8.16287882858964663038710959027145817e-01Q), static_cast(1.67355162876747144503180139830359482e+00Q), @@ -206,12 +205,9 @@ class hermite_detail }; return data; } - static std::array const & weights() + static std::array const & weights() { - static std::array data = { - static_cast(9.71781245099519154149424255938959644e-04Q), - static_cast(5.45155828191270305921785688416951260e-02Q), - static_cast(4.25607252610127800520317466666391035e-01Q), + static std::array data = { static_cast(8.10264617556807326764876563813094941e-01Q), static_cast(4.25607252610127800520317466666391035e-01Q), static_cast(5.45155828191270305921785688416951260e-02Q), @@ -223,16 +219,13 @@ class hermite_detail #endif template -class hermite_detail +class hermite_detail { using storage_type = typename hermite_constant_category::storage_type; public: - static std::array const & abscissa() + static std::array const & abscissa() { - static std::array data = { - BOOST_MATH_HUGE_CONSTANT(T, 0, -2.6519613568352334924470820065166161144381584786255294172031003071471600949016031668052906818697878641942881494964276e+00), - BOOST_MATH_HUGE_CONSTANT(T, 0, -1.6735516287674714450318013983035948191078100577354089269242175099937138333690347986612540761793805055242734953861690e+00), - BOOST_MATH_HUGE_CONSTANT(T, 0, -8.1628788285896466303871095902714581674288940037863615684472203343594907048766511668519794976704116666704491757953733e-01), + static std::array data = { BOOST_MATH_HUGE_CONSTANT(T, 0, 0.0000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000e+00), BOOST_MATH_HUGE_CONSTANT(T, 0, 8.1628788285896466303871095902714581674288940037863615684472203343594907048766511668519794976704116666704491757953733e-01), BOOST_MATH_HUGE_CONSTANT(T, 0, 1.6735516287674714450318013983035948191078100577354089269242175099937138333690347986612540761793805055242734953861690e+00), @@ -240,12 +233,9 @@ class hermite_detail }; return data; } - static std::array const & weights() + static std::array const & weights() { - static std::array data = { - BOOST_MATH_HUGE_CONSTANT(T, 0, 9.7178124509951915414942425593895964444240121701420164980001433798334142698580146031198718271051220413580750089473744e-04), - BOOST_MATH_HUGE_CONSTANT(T, 0, 5.4515582819127030592178568841695125960899915859346492986586846631369448562711195465546454603039043598893480308577392e-02), - BOOST_MATH_HUGE_CONSTANT(T, 0, 4.2560725261012780052031746666639103543970682101400152587878953913358027297072602373305673647185542480994656431711885e-01), + static std::array data = { BOOST_MATH_HUGE_CONSTANT(T, 0, 8.1026461755680732676487656381309494070745117994166268718345498964704515867018614005712030022333470419027511937110060e-01), BOOST_MATH_HUGE_CONSTANT(T, 0, 4.2560725261012780052031746666639103543970682101400152587878953913358027297072602373305673647185542480994656431711885e-01), BOOST_MATH_HUGE_CONSTANT(T, 0, 5.4515582819127030592178568841695125960899915859346492986586846631369448562711195465546454603039043598893480308577392e-02), @@ -257,39 +247,29 @@ class hermite_detail #ifndef BOOST_HAS_FLOAT128 template -class hermite_detail +class hermite_detail { using storage_type = typename hermite_constant_category::storage_type; public: - static std::array const & abscissa() + static std::array const & abscissa() { - static std::array data = { - static_cast(-3.43615911883773760332672549431912138e+00L), - static_cast(-2.53273167423278979640896079775479348e+00L), - static_cast(-1.75668364929988177345140122010615676e+00L), - static_cast(-1.03661082978951365417749191675920902e+00L), - static_cast(-3.42901327223704608789165025557258031e-01L), - static_cast(3.42901327223704608789165025557258031e-01L), - static_cast(1.03661082978951365417749191675920902e+00L), - static_cast(1.75668364929988177345140122010615676e+00L), - static_cast(2.53273167423278979640896079775479348e+00L), - static_cast(3.43615911883773760332672549431912138e+00L), + static std::array data = { + static_cast(3.42901327223704608789165025557258046e-01L), + static_cast(1.03661082978951365417749191675920910e+00L), + static_cast(1.75668364929988177345140122010615672e+00L), + static_cast(2.53273167423278979640896079775479347e+00L), + static_cast(3.43615911883773760332672549431912143e+00L), }; return data; } - static std::array const & weights() + static std::array const & weights() { - static std::array data = { - static_cast(7.64043285523262062915936785959522211e-06L), - static_cast(1.34364574678123269220156558584591387e-03L), - static_cast(3.38743944554810631361647312775859737e-02L), - static_cast(2.40138611082314686416523295005861395e-01L), - static_cast(6.10862633735325798783564990433419713e-01L), - static_cast(6.10862633735325798783564990433419713e-01L), - static_cast(2.40138611082314686416523295005861395e-01L), - static_cast(3.38743944554810631361647312775859737e-02L), - static_cast(1.34364574678123269220156558584591387e-03L), - static_cast(7.64043285523262062915936785959522211e-06L), + static std::array data = { + static_cast(6.10862633735325798783564990433419732e-01L), + static_cast(2.40138611082314686416523295005861392e-01L), + static_cast(3.38743944554810631361647312775859719e-02L), + static_cast(1.34364574678123269220156558584591379e-03L), + static_cast(7.64043285523262062915936785959522150e-06L), }; return data; } @@ -297,39 +277,29 @@ class hermite_detail #else template -class hermite_detail +class hermite_detail { using storage_type = typename hermite_constant_category::storage_type; public: - static std::array const & abscissa() + static std::array const & abscissa() { - static std::array data = { - static_cast(-3.43615911883773760332672549431912138e+00Q), - static_cast(-2.53273167423278979640896079775479348e+00Q), - static_cast(-1.75668364929988177345140122010615676e+00Q), - static_cast(-1.03661082978951365417749191675920902e+00Q), - static_cast(-3.42901327223704608789165025557258031e-01Q), - static_cast(3.42901327223704608789165025557258031e-01Q), - static_cast(1.03661082978951365417749191675920902e+00Q), - static_cast(1.75668364929988177345140122010615676e+00Q), - static_cast(2.53273167423278979640896079775479348e+00Q), - static_cast(3.43615911883773760332672549431912138e+00Q), + static std::array data = { + static_cast(3.42901327223704608789165025557258046e-01Q), + static_cast(1.03661082978951365417749191675920910e+00Q), + static_cast(1.75668364929988177345140122010615672e+00Q), + static_cast(2.53273167423278979640896079775479347e+00Q), + static_cast(3.43615911883773760332672549431912143e+00Q), }; return data; } - static std::array const & weights() + static std::array const & weights() { - static std::array data = { - static_cast(7.64043285523262062915936785959522211e-06Q), - static_cast(1.34364574678123269220156558584591387e-03Q), - static_cast(3.38743944554810631361647312775859737e-02Q), - static_cast(2.40138611082314686416523295005861395e-01Q), - static_cast(6.10862633735325798783564990433419713e-01Q), - static_cast(6.10862633735325798783564990433419713e-01Q), - static_cast(2.40138611082314686416523295005861395e-01Q), - static_cast(3.38743944554810631361647312775859737e-02Q), - static_cast(1.34364574678123269220156558584591387e-03Q), - static_cast(7.64043285523262062915936785959522211e-06Q), + static std::array data = { + static_cast(6.10862633735325798783564990433419732e-01Q), + static_cast(2.40138611082314686416523295005861392e-01Q), + static_cast(3.38743944554810631361647312775859719e-02Q), + static_cast(1.34364574678123269220156558584591379e-03Q), + static_cast(7.64043285523262062915936785959522150e-06Q), }; return data; } @@ -337,39 +307,29 @@ class hermite_detail #endif template -class hermite_detail +class hermite_detail { using storage_type = typename hermite_constant_category::storage_type; public: - static std::array const & abscissa() + static std::array const & abscissa() { - static std::array data = { - BOOST_MATH_HUGE_CONSTANT(T, 0, -3.4361591188377376033267254943191213848406783093901772906871762944002361145649511153468627231618843292611320162500620e+00), - BOOST_MATH_HUGE_CONSTANT(T, 0, -2.5327316742327897964089607977547934803078465081567249459332874134930859380473865340227292325941771646240824732515703e+00), - BOOST_MATH_HUGE_CONSTANT(T, 0, -1.7566836492998817734514012201061567632954744937388471000872849560961705796098190771278266711469838314640623161517203e+00), - BOOST_MATH_HUGE_CONSTANT(T, 0, -1.0366108297895136541774919167592090162982561106572157244372326874036404187243878992980306028243900733070850998719429e+00), - BOOST_MATH_HUGE_CONSTANT(T, 0, -3.4290132722370460878916502555725803120830265867773062820577103907292018192312101560713520305443248063627958675235231e-01), - BOOST_MATH_HUGE_CONSTANT(T, 0, 3.4290132722370460878916502555725803120830265867773062820577103907292018192312101560713520305443248063627958675235231e-01), - BOOST_MATH_HUGE_CONSTANT(T, 0, 1.0366108297895136541774919167592090162982561106572157244372326874036404187243878992980306028243900733070850998719429e+00), - BOOST_MATH_HUGE_CONSTANT(T, 0, 1.7566836492998817734514012201061567632954744937388471000872849560961705796098190771278266711469838314640623161517203e+00), - BOOST_MATH_HUGE_CONSTANT(T, 0, 2.5327316742327897964089607977547934803078465081567249459332874134930859380473865340227292325941771646240824732515703e+00), - BOOST_MATH_HUGE_CONSTANT(T, 0, 3.4361591188377376033267254943191213848406783093901772906871762944002361145649511153468627231618843292611320162500620e+00), + static std::array data = { + BOOST_MATH_HUGE_CONSTANT(T, 0, 3.4290132722370460878916502555725804574577169343357617326210300863769016299631431365924072451889514923095703125000000e-01), + BOOST_MATH_HUGE_CONSTANT(T, 0, 1.0366108297895136541774919167592091040917543343984240793575756300019937924949964269671909278258681297302246093750000e+00), + BOOST_MATH_HUGE_CONSTANT(T, 0, 1.7566836492998817734514012201061567205847985641023178720560650082445006137099596799089340493083000183105468750000000e+00), + BOOST_MATH_HUGE_CONSTANT(T, 0, 2.5327316742327897964089607977547934701687784149754958919647942962061508073712268185317952884361147880554199218750000e+00), + BOOST_MATH_HUGE_CONSTANT(T, 0, 3.4361591188377376033267254943191214342191179952800433354144214837346850421884170145858661271631717681884765625000000e+00), }; return data; } - static std::array const & weights() + static std::array const & weights() { - static std::array data = { - BOOST_MATH_HUGE_CONSTANT(T, 0, 7.6404328552326206291593678595952221082891147358323631952885551280911633480005783690893185993195845416600950144409630e-06), - BOOST_MATH_HUGE_CONSTANT(T, 0, 1.3436457467812326922015655858459138698692732786257703646162498020762483818413155201782152234375992775169973716992552e-03), - BOOST_MATH_HUGE_CONSTANT(T, 0, 3.3874394455481063136164731277585973698161942481867016739715157699771381277813573402914401584449869233310982466662099e-02), - BOOST_MATH_HUGE_CONSTANT(T, 0, 2.4013861108231468641652329500586139537004249390192739899066056147994315237169314350600241149865696523148172530499056e-01), - BOOST_MATH_HUGE_CONSTANT(T, 0, 6.1086263373532579878356499043341971323859272928403754564871663738953676910082005768001121074412857938126204657377493e-01), - BOOST_MATH_HUGE_CONSTANT(T, 0, 6.1086263373532579878356499043341971323859272928403754564871663738953676910082005768001121074412857938126204657377493e-01), - BOOST_MATH_HUGE_CONSTANT(T, 0, 2.4013861108231468641652329500586139537004249390192739899066056147994315237169314350600241149865696523148172530499056e-01), - BOOST_MATH_HUGE_CONSTANT(T, 0, 3.3874394455481063136164731277585973698161942481867016739715157699771381277813573402914401584449869233310982466662099e-02), - BOOST_MATH_HUGE_CONSTANT(T, 0, 1.3436457467812326922015655858459138698692732786257703646162498020762483818413155201782152234375992775169973716992552e-03), - BOOST_MATH_HUGE_CONSTANT(T, 0, 7.6404328552326206291593678595952221082891147358323631952885551280911633480005783690893185993195845416600950144409630e-06), + static std::array data = { + BOOST_MATH_HUGE_CONSTANT(T, 0, 6.1086263373532579878356499043341973238818374059758277936763908978968140237042483420282223960384726524353027343750000e-01), + BOOST_MATH_HUGE_CONSTANT(T, 0, 2.4013861108231468641652329500586139164344096899966477839326207760887354125331483167826718272408470511436462402343750e-01), + BOOST_MATH_HUGE_CONSTANT(T, 0, 3.3874394455481063136164731277585971926477457178246590999848008213597529213016723570461863346281461417675018310546875e-02), + BOOST_MATH_HUGE_CONSTANT(T, 0, 1.3436457467812326922015655858459137918851865251396942830563543358449043010113127527560550333873834460973739624023438e-03), + BOOST_MATH_HUGE_CONSTANT(T, 0, 7.6404328552326206291593678595952215021994824892912048011683019177853182368166190746930355182087168941507115960121155e-06), }; return data; } @@ -377,20 +337,13 @@ class hermite_detail #ifndef BOOST_HAS_FLOAT128 template -class hermite_detail +class hermite_detail { using storage_type = typename hermite_constant_category::storage_type; public: - static std::array const & abscissa() + static std::array const & abscissa() { - static std::array data = { - static_cast(-4.49999070730939155366438053053482422e+00L), - static_cast(-3.66995037340445253472922383311568148e+00L), - static_cast(-2.96716692790560324848896036354980632e+00L), - static_cast(-2.32573248617385774545404479448534701e+00L), - static_cast(-1.71999257518648893241583152515256151e+00L), - static_cast(-1.13611558521092066631913490555610213e+00L), - static_cast(-5.65069583255575748526020337198188867e-01L), + static std::array data = { static_cast(0.00000000000000000000000000000000000e+00L), static_cast(5.65069583255575748526020337198188867e-01L), static_cast(1.13611558521092066631913490555610213e+00L), @@ -402,16 +355,9 @@ class hermite_detail }; return data; } - static std::array const & weights() + static std::array const & weights() { - static std::array data = { - static_cast(1.52247580425351702016062666964826618e-09L), - static_cast(1.05911554771106663577520791055016329e-06L), - static_cast(1.00004441232499868127296736176978509e-04L), - static_cast(2.77806884291277589607887049229213490e-03L), - static_cast(3.07800338725460822286814158757801456e-02L), - static_cast(1.58488915795935746883839384959994027e-01L), - static_cast(4.12028687498898627025891079567810319e-01L), + static std::array data = { static_cast(5.64100308726417532852625797339963533e-01L), static_cast(4.12028687498898627025891079567810319e-01L), static_cast(1.58488915795935746883839384959994027e-01L), @@ -427,20 +373,13 @@ class hermite_detail #else template -class hermite_detail +class hermite_detail { using storage_type = typename hermite_constant_category::storage_type; public: - static std::array const & abscissa() + static std::array const & abscissa() { - static std::array data = { - static_cast(-4.49999070730939155366438053053482422e+00Q), - static_cast(-3.66995037340445253472922383311568148e+00Q), - static_cast(-2.96716692790560324848896036354980632e+00Q), - static_cast(-2.32573248617385774545404479448534701e+00Q), - static_cast(-1.71999257518648893241583152515256151e+00Q), - static_cast(-1.13611558521092066631913490555610213e+00Q), - static_cast(-5.65069583255575748526020337198188867e-01Q), + static std::array data = { static_cast(0.00000000000000000000000000000000000e+00Q), static_cast(5.65069583255575748526020337198188867e-01Q), static_cast(1.13611558521092066631913490555610213e+00Q), @@ -452,16 +391,9 @@ class hermite_detail }; return data; } - static std::array const & weights() + static std::array const & weights() { - static std::array data = { - static_cast(1.52247580425351702016062666964826618e-09Q), - static_cast(1.05911554771106663577520791055016329e-06Q), - static_cast(1.00004441232499868127296736176978509e-04Q), - static_cast(2.77806884291277589607887049229213490e-03Q), - static_cast(3.07800338725460822286814158757801456e-02Q), - static_cast(1.58488915795935746883839384959994027e-01Q), - static_cast(4.12028687498898627025891079567810319e-01Q), + static std::array data = { static_cast(5.64100308726417532852625797339963533e-01Q), static_cast(4.12028687498898627025891079567810319e-01Q), static_cast(1.58488915795935746883839384959994027e-01Q), @@ -477,20 +409,13 @@ class hermite_detail #endif template -class hermite_detail +class hermite_detail { using storage_type = typename hermite_constant_category::storage_type; public: - static std::array const & abscissa() + static std::array const & abscissa() { - static std::array data = { - BOOST_MATH_HUGE_CONSTANT(T, 0, -4.4999907073093915536643805305348242199319830534949905938271742846019366228525630586264389556919589962713446071341276e+00), - BOOST_MATH_HUGE_CONSTANT(T, 0, -3.6699503734044525347292238331156814846958556260108268200827250197712929612325341525137596165272567610075980577793224e+00), - BOOST_MATH_HUGE_CONSTANT(T, 0, -2.9671669279056032484889603635498063155705164076381381507212727522118956916046294049934561649283587290854852693064039e+00), - BOOST_MATH_HUGE_CONSTANT(T, 0, -2.3257324861738577454540447944853470123070950095634828957721625948740238048583608727061841554643103134570736658336089e+00), - BOOST_MATH_HUGE_CONSTANT(T, 0, -1.7199925751864889324158315251525615148429195630613888498938216300095132434761832428439718888176122448729491067241320e+00), - BOOST_MATH_HUGE_CONSTANT(T, 0, -1.1361155852109206663191349055561021252835292947846650850648392599800715967101004409631966601830482707422651868052475e+00), - BOOST_MATH_HUGE_CONSTANT(T, 0, -5.6506958325557574852602033719818886681552078800505788680930739682745058470915375159973526283034610292989738379345761e-01), + static std::array data = { BOOST_MATH_HUGE_CONSTANT(T, 0, 0.0000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000e+00), BOOST_MATH_HUGE_CONSTANT(T, 0, 5.6506958325557574852602033719818886681552078800505788680930739682745058470915375159973526283034610292989738379345761e-01), BOOST_MATH_HUGE_CONSTANT(T, 0, 1.1361155852109206663191349055561021252835292947846650850648392599800715967101004409631966601830482707422651868052475e+00), @@ -502,16 +427,9 @@ class hermite_detail }; return data; } - static std::array const & weights() + static std::array const & weights() { - static std::array data = { - BOOST_MATH_HUGE_CONSTANT(T, 0, 1.5224758042535170201606266696482661827962160316014537930014349599535491582697947140375830838937903702636262186650310e-09), - BOOST_MATH_HUGE_CONSTANT(T, 0, 1.0591155477110666357752079105501632872963524018315499158761728455091303265515893684325618169279641299843625141139696e-06), - BOOST_MATH_HUGE_CONSTANT(T, 0, 1.0000444123249986812729673617697850851699387810876387789234107161752152460966503439378917351614005951603805751399265e-04), - BOOST_MATH_HUGE_CONSTANT(T, 0, 2.7780688429127758960788704922921349030124251134822666685717061210098979501372205289327949738880183107056484267658214e-03), - BOOST_MATH_HUGE_CONSTANT(T, 0, 3.0780033872546082228681415875780145569645078095357831670126204320943426296907483543543105075772562113024847790201415e-02), - BOOST_MATH_HUGE_CONSTANT(T, 0, 1.5848891579593574688383938495999402721002703620834420484147417162556491543397074116784210832986285797844312792251258e-01), - BOOST_MATH_HUGE_CONSTANT(T, 0, 4.1202868749889862702589107956781031893375232741692141877844273550979447720862335855236503495494156330750630603931556e-01), + static std::array data = { BOOST_MATH_HUGE_CONSTANT(T, 0, 5.6410030872641753285262579733996353292453477640072243023805413420648986959478382400247069852378235816199417789938317e-01), BOOST_MATH_HUGE_CONSTANT(T, 0, 4.1202868749889862702589107956781031893375232741692141877844273550979447720862335855236503495494156330750630603931556e-01), BOOST_MATH_HUGE_CONSTANT(T, 0, 1.5848891579593574688383938495999402721002703620834420484147417162556491543397074116784210832986285797844312792251258e-01), @@ -556,8 +474,9 @@ class hermite : public detail::hermite_detail(base::abscissa()[i])); - result += fp * static_cast(base::weights()[i]); - L1 += abs(fp) * static_cast(base::weights()[i]); + K fm = f(-static_cast(base::abscissa()[i])); + result += (fp + fm) * static_cast(base::weights()[i]); + L1 += (abs(fp) + abs(fm)) * static_cast(base::weights()[i]); } if (pL1) *pL1 = L1; From b8a9e8109ecbbae2f35055e4440ecd359763d214 Mon Sep 17 00:00:00 2001 From: Jacob Hass Date: Mon, 7 Sep 2026 11:55:41 -0700 Subject: [PATCH 3/4] Added generating hermite code --- tools/hermite.cpp | 173 ++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 173 insertions(+) create mode 100644 tools/hermite.cpp diff --git a/tools/hermite.cpp b/tools/hermite.cpp new file mode 100644 index 0000000000..1189b8b5e9 --- /dev/null +++ b/tools/hermite.cpp @@ -0,0 +1,173 @@ +// Copyright (c) 2017 John Maddock +// Use, modification and distribution are subject to the +// Boost Software License, Version 1.0. (See accompanying file +// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) + +#include +#include +#include +#include +#include +#include +#include +#include +#include + +template +void print_gauss_constants(const char* suffix, int prec, int tag) +{ + std::string suff = std::string(suffix) + ")"; + const char* suffix_append = suff.c_str(); + + auto ab = T::abscissa(); + auto w = T::weights(); + std::cout << std::setprecision(prec) << std::scientific; + std::size_t order = (ab[0] == 0) ? (ab.size() * 2) - 1 : ab.size() * 2; + std::cout << + "template \n" + "class hermite_detail\n" + "{\n" + " using storage_type = typename hermite_constant_category::storage_type; \n" + " public:\n" + " static std::array<" << (prec > 40 ? "T" : "storage_type") << ", " << ab.size() << "> const & abscissa()\n" + " {\n" + " static std::array<" << (prec > 40 ? "T" : "storage_type") << ", " << ab.size() << "> data = {\n"; + for (unsigned i = 0; i < ab.size(); ++i) + std::cout << " " << (prec > 40 ? "BOOST_MATH_HUGE_CONSTANT(T, 0, " : "static_cast(") << ab[i] << (prec > 40 ? ")" : suffix_append) << ",\n"; + std::cout << + "};\n" + " return data;\n" + " }\n" + " static std::array<" << (prec > 40 ? "T" : "storage_type") << ", " << w.size() << "> const & weights()\n" + " {\n" + " static std::array<" << (prec > 40 ? "T" : "storage_type") << ", " << w.size() << "> data = {\n"; + for (unsigned i = 0; i < w.size(); ++i) + std::cout << " " << (prec > 40 ? "BOOST_MATH_HUGE_CONSTANT(T, 0, " : "static_cast(") << w[i] << (prec > 40 ? ")" : suffix_append) << ",\n"; + + std::cout << " };\n" + " return data;\n" + " }\n" + " };\n\n"; +} + +template +void print_gauss_kronrod_constants(const char* suffix, int prec, int tag) +{ + auto ab = T::abscissa(); + auto w = T::weights(); + std::cout << std::setprecision(prec) << std::scientific; + std::size_t order = (ab.size() * 2) - 1; + std::cout << + " template \n" + " class gauss_kronrod_detail\n" + " {\n" + " public:\n" + " static " << (prec > 40 ? " " : "constexpr ") << "std::array const & abscissa()\n" + " {\n" + " static " << (prec > 40 ? " " : "constexpr ") << "std::array data = {\n"; + + for (unsigned i = 0; i < ab.size(); ++i) + std::cout << " " << (prec > 40 ? "BOOST_MATH_HUGE_CONSTANT(T, 0, " : "") << ab[i] << (prec > 40 ? ")" : suffix) << ",\n"; + + std::cout << " };\n" + " return data;\n" + " }\n" + " static " << (prec > 40 ? " " : "constexpr ") << "std::array const & weights()\n" + " {\n" + " static " << (prec > 40 ? " " : "constexpr ") << "std::array data = {\n"; + + for (unsigned i = 0; i < w.size(); ++i) + std::cout << " " << (prec > 40 ? "BOOST_MATH_HUGE_CONSTANT(T, 0, " : "") << w[i] << (prec > 40 ? ")" : suffix) << ",\n"; + + std::cout << " };\n" + " return data;\n" + " }\n" + " };\n\n"; +} + + +int main() +{ + using boost::math::quadrature::detail::hermite_detail; + typedef boost::multiprecision::number > mp_type; + + std::cout << "#ifndef BOOST_HAS_FLOAT128\n"; + print_gauss_constants >("L", 35, 0); + std::cout << "#else\n"; + print_gauss_constants >("Q", 35, 0); + std::cout << "#endif\n"; + print_gauss_constants >("", 115, 4); + + std::cout << "#ifndef BOOST_HAS_FLOAT128\n"; + print_gauss_constants >("L", 35, 0); + std::cout << "#else\n"; + print_gauss_constants >("Q", 35, 0); + std::cout << "#endif\n"; + print_gauss_constants >("", 115, 4); + + std::cout << "#ifndef BOOST_HAS_FLOAT128\n"; + print_gauss_constants >("L", 35, 0); + std::cout << "#else\n"; + print_gauss_constants >("Q", 35, 0); + std::cout << "#endif\n"; + print_gauss_constants >("", 115, 4); + + // std::cout << "#ifndef BOOST_HAS_FLOAT128\n"; + // print_gauss_constants >("L", 35, 0); + // std::cout << "#else\n"; + // print_gauss_constants >("Q", 35, 0); + // std::cout << "#endif\n"; + // print_gauss_constants >("", 115, 4); + + // std::cout << "#ifndef BOOST_HAS_FLOAT128\n"; + // print_gauss_constants >("L", 35, 0); + // std::cout << "#else\n"; + // print_gauss_constants >("Q", 35, 0); + // std::cout << "#endif\n"; + // print_gauss_constants >("", 115, 4); + + // std::cout << "#ifndef BOOST_HAS_FLOAT128\n"; + // print_gauss_constants >("L", 35, 0); + // std::cout << "#else\n"; + // print_gauss_constants >("Q", 35, 0); + // std::cout << "#endif\n"; + // print_gauss_constants >("", 115, 4); + + // print_gauss_kronrod_constants >("f", 9, 0); + // print_gauss_kronrod_constants >("", 17, 1); + // print_gauss_kronrod_constants >("L", 35, 2); + // print_gauss_kronrod_constants >("Q", 35, 3); + // print_gauss_kronrod_constants >("", 115, 4); + + // print_gauss_kronrod_constants >("f", 9, 0); + // print_gauss_kronrod_constants >("", 17, 1); + // print_gauss_kronrod_constants >("L", 35, 2); + // print_gauss_kronrod_constants >("Q", 35, 3); + // print_gauss_kronrod_constants >("", 115, 4); + + // print_gauss_kronrod_constants >("f", 9, 0); + // print_gauss_kronrod_constants >("", 17, 1); + // print_gauss_kronrod_constants >("L", 35, 2); + // print_gauss_kronrod_constants >("Q", 35, 3); + // print_gauss_kronrod_constants >("", 115, 4); + + // print_gauss_kronrod_constants >("f", 9, 0); + // print_gauss_kronrod_constants >("", 17, 1); + // print_gauss_kronrod_constants >("L", 35, 2); + // print_gauss_kronrod_constants >("Q", 35, 3); + // print_gauss_kronrod_constants >("", 115, 4); + + // print_gauss_kronrod_constants >("f", 9, 0); + // print_gauss_kronrod_constants >("", 17, 1); + // print_gauss_kronrod_constants >("L", 35, 2); + // print_gauss_kronrod_constants >("Q", 35, 3); + // print_gauss_kronrod_constants >("", 115, 4); + + // print_gauss_kronrod_constants >("f", 9, 0); + // print_gauss_kronrod_constants >("", 17, 1); + // print_gauss_kronrod_constants >("L", 35, 2); + // print_gauss_kronrod_constants >("Q", 35, 3); + // print_gauss_kronrod_constants >("", 115, 4); + + return 0; +} From e402ba846d3c8eed50f1a7f443abcae866fc4f11 Mon Sep 17 00:00:00 2001 From: Jacob Hass Date: Tue, 8 Sep 2026 20:06:17 -0700 Subject: [PATCH 4/4] Added hermite root finding --- include/boost/math/quadrature/hermite.hpp | 109 +++++++++++++++------- 1 file changed, 76 insertions(+), 33 deletions(-) diff --git a/include/boost/math/quadrature/hermite.hpp b/include/boost/math/quadrature/hermite.hpp index 7c7d13f3f2..42a88dec7e 100644 --- a/include/boost/math/quadrature/hermite.hpp +++ b/include/boost/math/quadrature/hermite.hpp @@ -23,6 +23,7 @@ #include #include #include + #include #define EIGEN_SUPPORT #endif @@ -32,11 +33,52 @@ namespace boost { namespace math{ namespace quadrature{ namespace detail { #if defined(EIGEN_SUPPORT) && !defined(BOOST_MATH_GAUSS_NO_COMPUTE_ON_DEMAND) +template +struct hermite_functor +{ + hermite_functor(unsigned const& N_) : N(N_) {}; + + std::pair operator()(Real const& x) + { + Real fx = boost::math::hermite(N, x); + Real fdx = Real(2) * Real(N) * boost::math::hermite(N-1, x); + return std::make_pair(fx, fdx); + } + +private: + unsigned N; +}; + +template +Real hermite_minimizer(unsigned N, Real guess) +{ + Real min; + Real max; + if (guess > 0) + { + min = guess / 1.1; // Minimum possible value is half our guess. + max = guess * 1.1; // Maximum possible value is twice our guess. + } + else + { + min = guess * 1.1; // Minimum possible value is half our guess. + max = guess / 1.1; // Maximum possible value is twice our guess. + } + + const int digits = std::numeric_limits::digits; + int get_digits = static_cast(digits); // Accuracy doubles with each step, so stop when we have + // just over half the digits correct. + const std::uintmax_t maxit = 10000; + std::uintmax_t it = maxit; + Real result = boost::math::tools::newton_raphson_iterate(hermite_functor(N), guess, min, max, get_digits, it); + return result; +} + template Real factorial(unsigned n) { - Real i = 1; - Real factorial = 1; + Real i = Real(1); + Real factorial = Real(1); while (i <= n) { @@ -77,7 +119,8 @@ class hermite_detail // Need to account for root near 0 if (roots(i) >= -std::numeric_limits::epsilon() * 1000) { - pos_roots.push_back(roots(i)); + Real root = hermite_minimizer(N, roots(i)); + pos_roots.push_back(root); } } @@ -254,22 +297,22 @@ class hermite_detail static std::array const & abscissa() { static std::array data = { - static_cast(3.42901327223704608789165025557258046e-01L), - static_cast(1.03661082978951365417749191675920910e+00L), - static_cast(1.75668364929988177345140122010615672e+00L), - static_cast(2.53273167423278979640896079775479347e+00L), - static_cast(3.43615911883773760332672549431912143e+00L), + static_cast(3.42901327223704608789165025557258031e-01L), + static_cast(1.03661082978951365417749191675920902e+00L), + static_cast(1.75668364929988177345140122010615676e+00L), + static_cast(2.53273167423278979640896079775479348e+00L), + static_cast(3.43615911883773760332672549431912138e+00L), }; return data; } static std::array const & weights() { static std::array data = { - static_cast(6.10862633735325798783564990433419732e-01L), - static_cast(2.40138611082314686416523295005861392e-01L), - static_cast(3.38743944554810631361647312775859719e-02L), - static_cast(1.34364574678123269220156558584591379e-03L), - static_cast(7.64043285523262062915936785959522150e-06L), + static_cast(6.10862633735325798783564990433419713e-01L), + static_cast(2.40138611082314686416523295005861395e-01L), + static_cast(3.38743944554810631361647312775859737e-02L), + static_cast(1.34364574678123269220156558584591387e-03L), + static_cast(7.64043285523262062915936785959522211e-06L), }; return data; } @@ -284,22 +327,22 @@ class hermite_detail static std::array const & abscissa() { static std::array data = { - static_cast(3.42901327223704608789165025557258046e-01Q), - static_cast(1.03661082978951365417749191675920910e+00Q), - static_cast(1.75668364929988177345140122010615672e+00Q), - static_cast(2.53273167423278979640896079775479347e+00Q), - static_cast(3.43615911883773760332672549431912143e+00Q), + static_cast(3.42901327223704608789165025557258031e-01Q), + static_cast(1.03661082978951365417749191675920902e+00Q), + static_cast(1.75668364929988177345140122010615676e+00Q), + static_cast(2.53273167423278979640896079775479348e+00Q), + static_cast(3.43615911883773760332672549431912138e+00Q), }; return data; } static std::array const & weights() { static std::array data = { - static_cast(6.10862633735325798783564990433419732e-01Q), - static_cast(2.40138611082314686416523295005861392e-01Q), - static_cast(3.38743944554810631361647312775859719e-02Q), - static_cast(1.34364574678123269220156558584591379e-03Q), - static_cast(7.64043285523262062915936785959522150e-06Q), + static_cast(6.10862633735325798783564990433419713e-01Q), + static_cast(2.40138611082314686416523295005861395e-01Q), + static_cast(3.38743944554810631361647312775859737e-02Q), + static_cast(1.34364574678123269220156558584591387e-03Q), + static_cast(7.64043285523262062915936785959522211e-06Q), }; return data; } @@ -314,22 +357,22 @@ class hermite_detail static std::array const & abscissa() { static std::array data = { - BOOST_MATH_HUGE_CONSTANT(T, 0, 3.4290132722370460878916502555725804574577169343357617326210300863769016299631431365924072451889514923095703125000000e-01), - BOOST_MATH_HUGE_CONSTANT(T, 0, 1.0366108297895136541774919167592091040917543343984240793575756300019937924949964269671909278258681297302246093750000e+00), - BOOST_MATH_HUGE_CONSTANT(T, 0, 1.7566836492998817734514012201061567205847985641023178720560650082445006137099596799089340493083000183105468750000000e+00), - BOOST_MATH_HUGE_CONSTANT(T, 0, 2.5327316742327897964089607977547934701687784149754958919647942962061508073712268185317952884361147880554199218750000e+00), - BOOST_MATH_HUGE_CONSTANT(T, 0, 3.4361591188377376033267254943191214342191179952800433354144214837346850421884170145858661271631717681884765625000000e+00), + BOOST_MATH_HUGE_CONSTANT(T, 0, 3.4290132722370460878916502555725803120830265867773062820577103907292018192312101560713520305443248063627958675235231e-01), + BOOST_MATH_HUGE_CONSTANT(T, 0, 1.0366108297895136541774919167592090162982561106572157244372326874036404187243878992980306028243900733070850998719429e+00), + BOOST_MATH_HUGE_CONSTANT(T, 0, 1.7566836492998817734514012201061567632954744937388471000872849560961705796098190771278266711469838314640623161517203e+00), + BOOST_MATH_HUGE_CONSTANT(T, 0, 2.5327316742327897964089607977547934803078465081567249459332874134930859380473865340227292325941771646240824732515703e+00), + BOOST_MATH_HUGE_CONSTANT(T, 0, 3.4361591188377376033267254943191213848406783093901772906871762944002361145649511153468627231618843292611320162500620e+00), }; return data; } static std::array const & weights() { static std::array data = { - BOOST_MATH_HUGE_CONSTANT(T, 0, 6.1086263373532579878356499043341973238818374059758277936763908978968140237042483420282223960384726524353027343750000e-01), - BOOST_MATH_HUGE_CONSTANT(T, 0, 2.4013861108231468641652329500586139164344096899966477839326207760887354125331483167826718272408470511436462402343750e-01), - BOOST_MATH_HUGE_CONSTANT(T, 0, 3.3874394455481063136164731277585971926477457178246590999848008213597529213016723570461863346281461417675018310546875e-02), - BOOST_MATH_HUGE_CONSTANT(T, 0, 1.3436457467812326922015655858459137918851865251396942830563543358449043010113127527560550333873834460973739624023438e-03), - BOOST_MATH_HUGE_CONSTANT(T, 0, 7.6404328552326206291593678595952215021994824892912048011683019177853182368166190746930355182087168941507115960121155e-06), + BOOST_MATH_HUGE_CONSTANT(T, 0, 6.1086263373532579878356499043341971323859272928403754564871663738953676910082005768001121074412857938126204657377493e-01), + BOOST_MATH_HUGE_CONSTANT(T, 0, 2.4013861108231468641652329500586139537004249390192739899066056147994315237169314350600241149865696523148172530499056e-01), + BOOST_MATH_HUGE_CONSTANT(T, 0, 3.3874394455481063136164731277585973698161942481867016739715157699771381277813573402914401584449869233310982466662099e-02), + BOOST_MATH_HUGE_CONSTANT(T, 0, 1.3436457467812326922015655858459138698692732786257703646162498020762483818413155201782152234375992775169973716992552e-03), + BOOST_MATH_HUGE_CONSTANT(T, 0, 7.6404328552326206291593678595952221082891147358323631952885551280911633480005783690893185993195845416600950144409630e-06), }; return data; }