diff --git a/export/CMakeLists.txt.export b/export/CMakeLists.txt.export index 972f3c2f2..4a94cf30e 100644 --- a/export/CMakeLists.txt.export +++ b/export/CMakeLists.txt.export @@ -237,6 +237,28 @@ configure_file( libint2_install_path(_libint2_installed_datadir "${LIBINT2_INSTALL_DATADIR}" "\$\{_IMPORT_PREFIX\}") include(srclist.cmake) +# Keep the truncated-Coulomb numerical kernel in its own translation unit, +# compiled with IEEE floating-point semantics: its error bounds rely on +# infinities, compensated sums and exact products, which value-unsafe modes +# (-ffast-math, and icx's default -fp-model=fast) break. Contraction of a*b+c +# into a fused multiply-add is disabled as well: GCC's default contracts +# across statements, which makes compensated sums inexact; the kernel calls +# std::fma itself where it wants one. The options come after +# CMAKE_CXX_FLAGS on the command line, so they override them. +if (MSVC) # includes clang-cl and icx-cl + set(_libint2_ieee_fp_options "/fp:precise") +elseif (CMAKE_CXX_COMPILER_ID STREQUAL "IntelLLVM") + set(_libint2_ieee_fp_options "-fp-model=precise;-ffp-contract=off") +elseif (CMAKE_CXX_COMPILER_ID MATCHES "GNU|Clang") + set(_libint2_ieee_fp_options "-fno-fast-math;-ffp-contract=off") +endif() +foreach(_source IN LISTS LIBINT2_LIBRARY_CXX_SRC) + if (_source MATCHES "truncated_coulomb") + set_source_files_properties(${_source} PROPERTIES + SKIP_UNITY_BUILD_INCLUSION ON + COMPILE_OPTIONS "${_libint2_ieee_fp_options}") + endif() +endforeach() list(LENGTH LIBINT2_LIBRARY_CXX_SRC _source_count) message(STATUS "Loading ${_source_count} library source files from LIBINT2_LIBRARY_CXX_SRC") message(DEBUG "LIBINT2_LIBRARY_CXX_SRC=${LIBINT2_LIBRARY_CXX_SRC}") diff --git a/export/tests/CMakeLists.txt b/export/tests/CMakeLists.txt index 94a0b74f5..23b3fb274 100644 --- a/export/tests/CMakeLists.txt +++ b/export/tests/CMakeLists.txt @@ -85,6 +85,7 @@ if (LIBINT2_REQUIRE_CXX_API) unit/test-core.cc unit/test-core-ints.cc unit/test-molden.cc + unit/test-truncated-coulomb.cc unit/test-permute.cc unit/test-precision.cc unit/test-shell-order.cc diff --git a/export/tests/unit/test-core-ints.cc b/export/tests/unit/test-core-ints.cc index 4c9f8e05f..ce9264403 100644 --- a/export/tests/unit/test-core-ints.cc +++ b/export/tests/unit/test-core-ints.cc @@ -92,6 +92,25 @@ TEST_CASE("Boys reference values", "[core-ints][precision]") { #endif // LIBINT_HAS_MPFR +TEST_CASE("Boys interpolation endpoint", "[core-ints][precision]") { + const auto evaluator = libint2::FmEval_Chebyshev7::instance(40); + // The table covers [0,117); its right endpoint uses the asymptotic branch. + for (double T : + {std::nextafter(117., 0.), 117., std::nextafter(117., 118.)}) { + double values[41]; + long double reference[41]; + libint2::FmEval_Reference2::eval(reference, T, 40); + for (int last : {0, 40}) { + evaluator->eval(values, T, last); + for (int m = 0; m <= last; ++m) { + INFO("T=" << T << " m=" << m); + CHECK(std::abs(values[m] - reference[m]) <= + std::numeric_limits::epsilon()); + } + } + } +} + TEST_CASE("Slater/Yukawa core integral values", "[core-ints]") { using scalar_type = libint2::scalar_type; diff --git a/export/tests/unit/test-truncated-coulomb.cc b/export/tests/unit/test-truncated-coulomb.cc new file mode 100644 index 000000000..1e36fc821 --- /dev/null +++ b/export/tests/unit/test-truncated-coulomb.cc @@ -0,0 +1,460 @@ +/* Independent truncated-Coulomb tests. LGPL-3.0-or-later; see COPYING.LESSER. + */ +#include + +#include +#include +#include +#include +#include +#include + +#include "catch.hpp" +#include "truncated-coulomb-reference.h" +#include "truncated-coulomb-vrr-reference.h" + +namespace { +template +void check_truncated_coulomb_radial(Evaluator& eval) { + const double eps = std::numeric_limits::epsilon(); + std::array values; + for (const auto& ref : truncated_coulomb_reference) { + // Check that the result is independent of the requested prefix, including + // the high-order region where direct upward recurrence loses accuracy. + for (int mmax : {0, 1, 2, 8, 9, 21, 40, 64}) { + eval.eval(values.data(), 1., ref.t, mmax, ref.r); + for (int m = 0; m <= mmax; ++m) { + INFO("R=" << ref.r << " T=" << ref.t << " mmax=" << mmax << " m=" << m); + CHECK(std::isfinite(values[m])); + CHECK(std::abs(values[m] - ref.g[m]) <= eps); + } + // Equal dimensionless R and T at a different physical scale. + eval.eval(values.data(), 4., ref.t, mmax, ref.r / 2); + for (int m = 0; m <= mmax; ++m) + CHECK(std::abs(values[m] - ref.g[m]) <= eps); + } + } + for (double t : {0., 0.1, 10., 100., 1e5}) { + eval.eval(values.data(), 1., t, 64, 0.); + for (auto value : values) CHECK(value == 0.); + eval.eval(values.data(), 1., t, 64, + std::numeric_limits::infinity()); + std::array ref; + libint2::FmEval_Reference2::eval(ref.data(), t, 64); + for (int m = 0; m <= 64; ++m) CHECK(std::abs(values[m] - ref[m]) <= eps); + } + CHECK_THROWS_AS(eval.eval(values.data(), 1., 0., 1, -1.), + std::invalid_argument); + CHECK_THROWS_AS(eval.eval(values.data(), 1., 0., 1, + std::numeric_limits::quiet_NaN()), + std::invalid_argument); + CHECK_THROWS_AS(eval.eval(values.data(), 0., 0., 1, 1.), + std::invalid_argument); +} +} // namespace + +TEST_CASE("Truncated Coulomb auxiliary integrals", + "[core-ints][truncated-coulomb]") { + const double eps = std::numeric_limits::epsilon(); + libint2::os_core_ints::truncated_coulomb_gm_eval tc(192, eps); + check_truncated_coulomb_radial(tc); + + // Independent defining-erf references, generated at 100 decimal digits and + // checked at 140 digits. These inputs exposed cancellation at m=40 and a + // parasitic high-order solution at m=192; the criterion is absolute error. + std::array values; + decltype(tc)::scratch_type scratch(192); + tc.eval(values.data(), 1., 15., 40, 1., scratch); + CHECK(std::abs(values[40] - 9.1989225488094219e-09) <= eps); + tc.eval(values.data(), 1., 59., 192, .1, scratch); + CHECK(std::abs(values[192] - 1.4877535827085252e-28) <= eps); + // Reusing a larger workspace must not retain the preceding hierarchy. + tc.eval(values.data(), 1., 15., 40, 1., scratch); + CHECK(std::abs(values[40] - 9.1989225488094219e-09) <= eps); + // A reserve beyond the native Boys table must still use its supported + // orders for later, lower-order requests, even if the singleton has grown. + std::array coulomb; + libint2::FmEval_Chebyshev7::instance(8)->eval(coulomb.data(), .24, 8); + tc.eval(values.data(), 1., .24, 8, std::numeric_limits::infinity(), + scratch); + for (int m = 0; m <= 8; ++m) CHECK(values[m] == coulomb[m]); + tc.eval(values.data(), 1., .24, 8, 20., scratch); + for (int m = 0; m <= 8; ++m) CHECK(std::abs(values[m] - coulomb[m]) <= eps); +} + +// 0 < T < min(n + 1/2, 12) uses the double-word origin expansion: its top +// of range at a high order, and R well beyond sqrt(T), where U_k and V_k +// grow with the dominant solution and the error estimate needs its nested +// form. Defining-erf references at 50 digits; the criterion is absolute. +// Requests with T >= n + 1/2 take the other paths and must agree. +TEST_CASE("Truncated Coulomb auxiliary integrals by the origin expansion", + "[core-ints][truncated-coulomb]") { + const double eps = std::numeric_limits::epsilon(); + struct Reference { + double r, t; + int n; + std::vector g; + }; + const Reference refs[] = { + {10.17, + 4.39, + 8, + {0.42168486463540115, 0.046615505172546444, 0.014515465409703409, + 0.0068538266318990484, 0.0040519427338097875, 0.0027410883125107839, + 0.0020217815737101615, 0.0015811425099987123, 0.0012888847869632435}}, + {2.5, + 0.75, + 16, + {0.78684342664308601, 0.23274225857973911, 0.10702959412734653, + 0.079677306351208520, 0.064820966424317426, 0.050486149778876075, + 0.039215987596959884, 0.032395090517522593, 0.028933570819893606, + 0.026978128157491320, 0.025178016413891936, 0.023023311904602649, + 0.020622817168460760, 0.018330621074997705, 0.016456474963152301, + 0.015131173955875089, 0.014302465367542966}}, + {4.0, + 6.0, + 24, + {0.35648410865939795, 0.033264962713753418, 0.0057267648960145567, + 0.0031756712700991371, 0.0015358745407417276, 0.00079818953623173328, + 0.00057881633889808070, 0.00045803099002024138, 0.00034837190263657074, + 0.00027074722046168703, 0.00022641815596211149, 0.00019844919959376965, + 0.00017484643544159914, 0.00015342899987390319, 0.00013583485881713887, + 0.00012261888641286325, 0.00011258053758827510, 0.00010414278242677521, + 9.6430909906661019e-5, 8.9356055604776100e-5, 8.3142386188761031e-5, + 7.7896557283820035e-5, 7.3481447750625323e-5, 6.9627983537232975e-5, + 6.6101460978547776e-5}}, + {1.3, + 11.9, + 32, + {0.00030356804207075164, 0.00021926625544449961, 0.00015621835795788411, + 0.00010956135496208932, 7.5444668313430262e-5, 5.0837636044396094e-5, + 3.3370582399328925e-5, 2.1204065307623918e-5, 1.2921709858899920e-5, + 7.4427271828251431e-6, 3.9508125676048421e-6, 1.8366272052554936e-6, + 6.5150563597032163e-7, 7.0405436508649165e-8, -1.3756457694838029e-7, + -1.3242897395803016e-7, -2.1602980738434684e-8, 1.2591410098219692e-7, + 2.6861185930436844e-7, 3.8415069602994759e-7, 4.6311306128179859e-7, + 5.0438309904018694e-7, 5.1179632360591788e-7, 4.9176966149080371e-7, + 4.5167919355463708e-7, 3.9880001985529897e-7, 3.3966139346921833e-7, + 2.7970195820779460e-7, 2.2313571606541575e-7, 1.7296021462410085e-7, + 1.3105520979003631e-7, 9.8333427916515778e-8, 7.4915618858296132e-8}}, + }; + libint2::os_core_ints::truncated_coulomb_gm_eval tc(32, eps); + decltype(tc)::scratch_type scratch(32); + std::array g; + for (const auto& ref : refs) { + for (int n : {ref.n, static_cast(ref.t - .5)}) { + tc.eval(g.data(), 1., ref.t, n, ref.r, scratch); + for (int m = 0; m <= n; ++m) { + INFO("R=" << ref.r << " T=" << ref.t << " n=" << n << " m=" << m); + CHECK(std::abs(g[m] - ref.g[m]) <= eps); + } + } + } +} + +// Just beyond the distance |R - sqrt(T)| at which an absolute criterion +// would switch to F_m or zero; the high orders are far below eps here, and +// G_16 even exceeds F_16 by 1e8 at T = 1e5. Defining-erf references at 50 +// digits; the criterion is relative to the larger of |G_m| and F_m. +TEST_CASE("Truncated Coulomb auxiliary integrals near the Boys and zero limits", + "[core-ints][truncated-coulomb]") { + constexpr int mmax = 16; + struct Reference { + double r, t, g[mmax + 1]; + }; + const Reference refs[] = { + {108.6, + 1e4, + {0.0088622692545275801, 4.4311346272637901e-7, 6.6467019408956851e-11, + 1.6616754852239213e-14, 5.8158641982837245e-18, 2.617138889227676e-21, + 1.4394263890752218e-24, 9.3562715289889424e-28, 7.0172036467416507e-31, + 5.9646230997349737e-34, 5.66639194437939e-37, 5.9497115711008598e-40, + 6.8421659681667526e-43, 8.5528911216056172e-46, 1.1532116269244219e-48, + 1.7822050813808079e-51, -5.6290925324490396e-54}}, + {91.4, + 1e4, + {2.1884352861136341e-36, 1.8957055345478377e-37, 1.6322092908399323e-38, + 1.3966766794916191e-39, 1.1876156724097307e-40, 1.0033637183153641e-41, + 8.4213758114545613e-43, 7.0208024419088905e-44, 5.8130393957925464e-45, + 4.7792858020065692e-46, 3.9011509332364499e-47, 3.1609337260052282e-48, + 2.5418476958424058e-49, 2.0281923028827299e-50, 1.6054734130845148e-51, + 1.2604767988555799e-52, 9.8129961037524071e-54}}, + {324.8277660168379, + 1e5, + {0.0028024956081989643, 1.4012478040994822e-8, 2.1018717061492233e-13, + 5.2546792653730582e-18, 1.8391377428805704e-22, 8.2761198429625666e-27, + 4.551865913629409e-31, 2.9587128438598042e-35, 2.2190346327157234e-39, + 1.8861794840531705e-43, 1.7918586678269412e-47, 1.8844588437658945e-51, + 1.4099508446526933e-55, 1.9074151329311718e-57, -4.648509778688959e-59, + 1.138728452152227e-60, -2.7627424580754614e-62}}, + {307.62776601683794, + 1e5, + {6.9204400160012852e-37, 1.8949573389055659e-38, 5.1553492922574712e-40, + 1.3933279360641161e-41, 3.7404711744624385e-43, 9.9727609630811688e-45, + 2.6403173634543515e-46, 6.9403414288238827e-48, 1.8110001939567021e-49, + 4.6902331769513531e-51, 1.2053967397591896e-52, 3.0735642877976119e-54, + 7.7739850384888677e-56, 1.9500309936894007e-57, 4.8499523714967992e-59, + 1.1957136764392447e-60, 2.9214574907968413e-62}}, + }; + libint2::os_core_ints::truncated_coulomb_gm_eval tc( + mmax, std::numeric_limits::epsilon()); + decltype(tc)::scratch_type scratch(mmax); + for (const auto& ref : refs) { + double g[mmax + 1]; + long double f[mmax + 1]; + tc.eval(g, 1., ref.t, mmax, ref.r, scratch); + libint2::FmEval_Reference2::eval(f, ref.t, mmax); + for (int m = 0; m <= mmax; ++m) { + INFO("R=" << ref.r << " T=" << ref.t << " m=" << m); + const double scale = + std::max(std::abs(ref.g[m]), static_cast(f[m])); + CHECK(std::abs(g[m] - ref.g[m]) <= 1e-12 * scale); + } + } +} + +// (L0|00)^(0), L = 0..mmax, by the Obara-Saika VRR with P = A, from G_m. +template +std::array vrr_l000(const double* g, double zeta, double rho, + double wp) { + std::array, N> a{}; + for (std::size_t m = 0; m < N; ++m) a[0][m] = g[m]; + for (std::size_t l = 0; l + 1 < N; ++l) + for (std::size_t m = 0; m + l + 1 < N; ++m) + a[l + 1][m] = + wp * a[l][m + 1] + + (l ? l / (2 * zeta) * (a[l - 1][m] - rho / zeta * a[l - 1][m + 1]) + : 0.); + std::array result; + for (std::size_t l = 0; l < N; ++l) result[l] = a[l][0]; + return result; +} + +// What must hold is the accuracy of the integrals, which neither a bound +// relative to each F_m nor an absolute one is the right measure of: far +// outside the cutoff band the high-order G_m enter integrals multiplied by +// small powers of the geometry, while near it an error of eps F_0 in G_4 can +// already exceed eps (L0|00) by 1e8 for tight primitives (the third +// reference geometry). Compare (L0|00) from evaluated and reference G_m, +// relative to the larger of the truncated and the Coulomb integral, against +// the double-precision VRR floor (~2^(L/2) eps; the tolerance leaves a factor +// of 8). +TEST_CASE("Truncated Coulomb integrals across the cutoff transition", + "[core-ints][truncated-coulomb]") { + constexpr int mmax = 16; + constexpr std::size_t N = mmax + 1; + const double eps = std::numeric_limits::epsilon(); + libint2::os_core_ints::truncated_coulomb_gm_eval tc(mmax, eps); + decltype(tc)::scratch_type scratch(mmax); + const auto boys = libint2::FmEval_Chebyshev7::instance(mmax); + for (const auto& ref : truncated_coulomb_vrr_reference) { + double g[N], f[N]; + tc.eval(g, 1., ref.t, mmax, ref.r, scratch); + boys->eval(f, ref.t, mmax); + const auto computed = vrr_l000(g, ref.zeta, ref.rho, ref.wp); + const auto exact = vrr_l000(ref.g, ref.zeta, ref.rho, ref.wp); + const auto coulomb = vrr_l000(f, ref.zeta, ref.rho, ref.wp); + for (int l = 0; l <= mmax; ++l) { + INFO("zeta=" << ref.zeta << " rho=" << ref.rho << " R=" << ref.r + << " T=" << ref.t << " L=" << l); + const double scale = std::max(std::abs(exact[l]), std::abs(coulomb[l])); + CHECK(std::abs(computed[l] - exact[l]) <= + 8 * eps * std::ldexp(1., l / 2) * scale); + } + } +} + +#ifdef LIBINT_INCLUDE_ERI +TEST_CASE("Truncated Coulomb Engine limits and derivatives", + "[engine][2-body][truncated-coulomb]") { + using namespace libint2; + for (const auto op : {Operator::truncated_coulomb}) { + REQUIRE(rank(op) == 2); + const double infinity = std::numeric_limits::infinity(); + Shell a{{0.8}, {{0, false, {1.0}}}, {{0.1, -0.2, 0.3}}}; + Shell b{{1.1}, {{0, false, {1.0}}}, {{-0.4, 0.3, 0.2}}}; + Shell c{{0.6}, {{0, false, {1.0}}}, {{0.5, 0.1, -0.1}}}; + Shell d{{1.3}, {{0, false, {1.0}}}, {{-0.1, -0.3, 0.4}}}; + + std::vector> brakets{ + {BraKet::xx_xx, LIBINT_INCLUDE_ERI}}; +#ifdef LIBINT_INCLUDE_ERI3 + brakets.emplace_back(BraKet::xs_xx, LIBINT_INCLUDE_ERI3); +#endif +#ifdef LIBINT_INCLUDE_ERI2 + brakets.emplace_back(BraKet::xs_xs, LIBINT_INCLUDE_ERI2); +#endif + for (const auto& entry : brakets) { + const auto braket = entry.first; + for (int deriv = 0; deriv <= std::min(2, entry.second); ++deriv) { + Engine trunc(op, 1, 1, deriv); + Engine coul(Operator::coulomb, 1, 1, deriv); + trunc.set(braket).set_params(infinity); + coul.set(braket); + auto compute = [&](Engine& e) { + if (braket == BraKet::xx_xx) + e.compute(a, b, c, d); + else if (braket == BraKet::xs_xx) + e.compute(a, c, d); + else + e.compute(a, c); + }; + compute(trunc); + compute(coul); + REQUIRE(trunc.results().size() == coul.results().size()); + for (std::size_t k = 0; k < trunc.results().size(); ++k) + CHECK(trunc.results()[k][0] == + Approx(coul.results()[k][0]).margin(2e-14)); + trunc.set_params(0.); + compute(trunc); + if (trunc.results()[0]) + for (auto target : trunc.results()) CHECK(target[0] == 0.); + } + } + + Engine e(op, 1, 0, 0); + CHECK_THROWS_AS(e.set_params(-1.), std::invalid_argument); + e.set_params(1.2); + e.compute(a, b, c, d); + const auto value = e.results()[0][0]; + e.compute(c, d, a, b); + CHECK(e.results()[0][0] == Approx(value).margin(2e-15)); + e.compute(b, a, d, c); + CHECK(e.results()[0][0] == Approx(value).margin(2e-15)); + + // Changing the operator must restore its default (zero) cutoff, rather + // than retaining the old core parameters. + e.set(Operator::coulomb); + e.set(op); + e.compute(a, b, c, d); + CHECK((!e.results()[0] || e.results()[0][0] == 0.)); + e.set_params(1.2); + e.compute(a, b, c, d); + CHECK(e.results()[0][0] == value); + +#if LIBINT_INCLUDE_ERI >= 1 + Engine grad(op, 1, 0, 1); + grad.set_params(1.2); + grad.compute(a, b, c, d); + std::vector analytic; + for (auto target : grad.results()) analytic.push_back(target[0]); + const double step = 1e-4; + for (int center = 0; center < 4; ++center) { + for (int xyz = 0; xyz < 3; ++xyz) { + std::array shells{{a, b, c, d}}; + auto sample = [&](double displacement) { + shells[center].O[xyz] += displacement; + e.compute(shells[0], shells[1], shells[2], shells[3]); + const auto v = e.results()[0][0]; + shells[center].O[xyz] -= displacement; + return v; + }; + const auto fd = (sample(-2 * step) - 8 * sample(-step) + + 8 * sample(step) - sample(2 * step)) / + (12 * step); + CHECK(analytic[3 * center + xyz] == Approx(fd).margin(3e-11)); + } + } +#endif +#if LIBINT_INCLUDE_ERI >= 2 + Engine hessian(op, 1, 0, 2); + hessian.set_params(1.2); + hessian.compute(a, b, c, d); + const double hxx = hessian.results()[0][0]; + const double h = 1e-3; + auto displaced = [&](double dx) { + Shell moved = a; + moved.O[0] += dx; + e.compute(moved, b, c, d); + return e.results()[0][0]; + }; + const auto fd2 = + (-displaced(2 * h) + 16 * displaced(h) - 30 * displaced(0) + + 16 * displaced(-h) - displaced(-2 * h)) / + (12 * h * h); + CHECK(hxx == Approx(fd2).margin(2e-9)); +#endif + } +} + +// (aa|bb) of normalized s primitives at distance D equals +// [erf(sqrt(mu)(Rc-D)) - erf(sqrt(mu)(Rc+D)) + 2 erf(sqrt(mu) D)] / (2D), +// mu = 4ab/(2a+2b): the radial average of theta(Rc-r)/r over the Gaussian +// distribution of r12. This checks the cutoff units and its sqrt(rho) +// scaling independently of the radial evaluator. mpmath values, 50 digits. +TEST_CASE("Truncated Coulomb Engine against closed-form (ss|ss)", + "[engine][2-body][truncated-coulomb]") { + using namespace libint2; + Shell a{{0.8}, {{0, false, {1.0}}}, {{0, 0, 0}}}; + Shell b{{1.3}, {{0, false, {1.0}}}, {{0, 0, 1.5}}}; + const std::pair refs[] = {{0.5, 0.031552595422039718}, + {1.2, 0.2011613560258891}, + {1.5, 0.31017100129777589}, + {2.5, 0.59039957931514569}, + {4.0, 0.64335169847410041}}; + Engine e(Operator::truncated_coulomb, 1, 0); + for (const auto& ref : refs) { + INFO("Rc=" << ref.first); + e.set_params(ref.first); + e.compute(a, a, b, b); + CHECK(e.results()[0][0] == Approx(ref.second).epsilon(1e-14)); + } +} + +TEST_CASE("Truncated Coulomb contractions and reduced centers", + "[engine][2-body][truncated-coulomb]") { + using namespace libint2; + for (const auto op : {Operator::truncated_coulomb}) { + const double inf = std::numeric_limits::infinity(); + for (bool pure : {false, true}) { + Shell a{{0.8, 1.7}, {{1, pure, {0.6, 0.4}}}, {{0.1, -0.2, 0.3}}}; + Shell b{{1.1}, {{0, pure, {1.0}}}, {{-0.4, 0.3, 0.2}}}; + Shell c{{0.6, 1.3}, {{1, pure, {0.7, 0.3}}}, {{0.5, 0.1, -0.1}}}; + Shell d{{1.3}, {{0, pure, {1.0}}}, {{-0.1, -0.3, 0.4}}}; + Engine e(op, 2, 0), coul(Operator::coulomb, 2, 1); + // Exercise resizing the radial evaluator after constructing an s engine. + e.set_max_l(1).set_params(inf); + e.compute(a, b, c, d); + coul.compute(a, b, c, d); + const auto size = a.size() * b.size() * c.size() * d.size(); + for (std::size_t i = 0; i < size; ++i) + CHECK(e.results()[0][i] == Approx(coul.results()[0][i]).margin(3e-14)); + e.set_params(1.2); + e.compute(a, b, c, d); + std::vector ref(e.results()[0], e.results()[0] + size); + e.compute(c, d, a, b); + for (std::size_t i = 0; i < a.size(); ++i) + for (std::size_t j = 0; j < c.size(); ++j) + CHECK(e.results()[0][j * a.size() + i] == + Approx(ref[i * c.size() + j]).margin(2e-14)); + +#ifdef LIBINT_INCLUDE_ERI3 + Engine three(op, 2, 1); + three.set(BraKet::xs_xx).set_params(1.2); + e.compute(a, Shell::unit(), c, d); + ref.assign(e.results()[0], + e.results()[0] + a.size() * c.size() * d.size()); + three.compute(a, c, d); + for (std::size_t i = 0; i < ref.size(); ++i) + CHECK(three.results()[0][i] == Approx(ref[i]).margin(2e-14)); +#endif +#ifdef LIBINT_INCLUDE_ERI2 + Engine two(op, 2, 1); + two.set(BraKet::xs_xs).set_params(1.2); + e.compute(a, Shell::unit(), c, Shell::unit()); + ref.assign(e.results()[0], e.results()[0] + a.size() * c.size()); + two.compute(a, c); + for (std::size_t i = 0; i < ref.size(); ++i) + CHECK(two.results()[0][i] == Approx(ref[i]).margin(2e-14)); + // Increase center count after creating a two-center core evaluator. + two.set(BraKet::xx_xx); + two.compute(a, b, c, d); + e.compute(a, b, c, d); + for (std::size_t i = 0; i < size; ++i) + CHECK(two.results()[0][i] == Approx(e.results()[0][i]).margin(2e-14)); +#endif + } + } +} + +#endif // LIBINT_INCLUDE_ERI diff --git a/export/tests/unit/truncated-coulomb-reference.h b/export/tests/unit/truncated-coulomb-reference.h new file mode 100644 index 000000000..8dea9092e --- /dev/null +++ b/export/tests/unit/truncated-coulomb-reference.h @@ -0,0 +1,391 @@ +// Generated by tools/generate_truncated_coulomb_reference.py +// mpmath 1.4.1; decimal precision 100; mmax 64. +// Independent mathematical reference; rounded once to double. +// clang-format off +struct TruncatedCoulombReference { + double r, t; + double g[65]; +}; +static const TruncatedCoulombReference truncated_coulomb_reference[] = { + {0.00000001, 0, { + 9.9999999999999998e-17, 9.9999999999999998e-17, 9.9999999999999998e-17, 9.9999999999999986e-17, + 9.9999999999999986e-17, 9.9999999999999986e-17, 9.9999999999999973e-17, 9.9999999999999973e-17, + 9.9999999999999973e-17, 9.9999999999999973e-17, 9.9999999999999961e-17, 9.9999999999999961e-17, + 9.9999999999999961e-17, 9.9999999999999961e-17, 9.9999999999999949e-17, 9.9999999999999949e-17, + 9.9999999999999949e-17, 9.9999999999999949e-17, 9.9999999999999936e-17, 9.9999999999999936e-17, + 9.9999999999999936e-17, 9.9999999999999924e-17, 9.9999999999999924e-17, 9.9999999999999924e-17, + 9.9999999999999924e-17, 9.9999999999999912e-17, 9.9999999999999912e-17, 9.9999999999999912e-17, + 9.9999999999999912e-17, 9.9999999999999899e-17, 9.9999999999999899e-17, 9.9999999999999899e-17, + 9.9999999999999887e-17, 9.9999999999999887e-17, 9.9999999999999887e-17, 9.9999999999999887e-17, + 9.9999999999999875e-17, 9.9999999999999875e-17, 9.9999999999999875e-17, 9.9999999999999875e-17, + 9.9999999999999862e-17, 9.9999999999999862e-17, 9.9999999999999862e-17, 9.999999999999985e-17, + 9.999999999999985e-17, 9.999999999999985e-17, 9.999999999999985e-17, 9.9999999999999838e-17, + 9.9999999999999838e-17, 9.9999999999999838e-17, 9.9999999999999838e-17, 9.9999999999999825e-17, + 9.9999999999999825e-17, 9.9999999999999825e-17, 9.9999999999999825e-17, 9.9999999999999813e-17, + 9.9999999999999813e-17, 9.9999999999999813e-17, 9.9999999999999801e-17, 9.9999999999999801e-17, + 9.9999999999999801e-17, 9.9999999999999801e-17, 9.9999999999999788e-17, 9.9999999999999788e-17, + 9.9999999999999788e-17, + }}, + {0.1, 0, { + 0.0099501662508319471, 0.0099170543086051031, 0.0098840305912597567, 0.0098510949098371964, + 0.0098182470757145219, 0.0097854869006041446, 0.0097528141965532687, 0.009720228775943392, + 0.0096877304514897955, 0.0096553190362410352, 0.0096229943435784428, 0.0095907561872156169, + 0.0095586043811979205, 0.0095265387399019741, 0.0094945590780351596, 0.0094626652106351121, + 0.0094308569530692221, 0.0093991341210341341, 0.009367496530555245, 0.0093359439979862069, + 0.0093044763400084284, 0.0092730933736305723, 0.0092417949161880657, 0.009210580785342596, + 0.009179450799081619, 0.0091484047757178638, 0.0091174425338888324, 0.0090865638925563142, + 0.0090557686710058868, 0.0090250566888464244, 0.0089944277660096023, 0.0089638817227494158, + 0.0089334183796416777, 0.008903037557583534, 0.0088727390777929727, 0.0088425227618083337, + 0.0088123884314878272, 0.0087823359090090369, 0.0087523650168684387, 0.0087224755778809119, + 0.008692667415179255, 0.0086629403522137018, 0.0086332942127514325, 0.0086037288208760943, + 0.0085742440009873178, 0.0085448395778002333, 0.008515515376344986, 0.0084862712219662649, + 0.0084571069403228064, 0.0084280223573869311, 0.0083990172994440546, 0.0083700915930922104, + 0.0083412450652415714, 0.0083124775431139743, 0.0082837888542424461, 0.0082551788264707202, + 0.0082266472879527659, 0.0081981940671523139, 0.0081698189928423783, 0.0081415218941047855, + 0.0081133026003297071, 0.0080851609412151745, 0.0080570967467666187, 0.0080291098472963941, + 0.0080012000734233073, + }}, + {1, 0, { + 0.63212055882855767, 0.45595981372381411, 0.32262648039048075, 0.22344025854231594, + 0.15120802374672862, 0.10000433058018862, 0.064980300377998915, 0.042197631674034751, + 0.02848480846740499, 0.021312886874780463, 0.018688517244869876, 0.019062108925738827, + 0.021249270887697073, 0.024363865431164783, 0.027761196196310099, 0.030990017427312607, + 0.033752200570640246, 0.035869028310361005, 0.037253206445975448, 0.037885791865920053, + 0.037797331422623308, 0.037052592834396839, 0.035738345796488352, 0.033953720165755596, + 0.031802729201131481, 0.029388600134982963, 0.026809602493065795, 0.024156107187537343, + 0.021508647037799988, 0.018936782536588016, 0.016498605836005265, 0.014240741499880463, + 0.012198724935482737, 0.010397658924200222, 0.0088530656293498607, 0.0075718661520183172, + 0.0065534323875186089, 0.0057906668354683324, 0.0052710753426255998, 0.0049778056956334632, + 0.0048906316979563059, 0.0049868680115388638, 0.0052422057534065912, 0.0056314627306501153, + 0.006129245381180367, 0.0067105220577882378, 0.0073511093342854619, 0.0080280746001585909, + 0.0087200594109054647, 0.0094075289339164009, 0.010072953426261503, 0.010700928046588253, + 0.011278237478388142, 0.011793871860944629, 0.012239000417567253, 0.01260690896442784, + 0.012892907200003589, 0.013094211334138016, 0.013209807233535984, 0.01324029885106485, + 0.013187746281257778, 0.013055497353657705, 0.012848016247125062, 0.012570712188472376, + 0.012229770892981396, + }}, + {3, 0, { + 0.99987659019591335, 0.3340326555564912, 0.19819821686033448, 0.14493395356020156, + 0.11081531935845891, 0.089761219458611943, 0.076844346152957294, 0.06740694631478121, + 0.059358726047453618, 0.052504049617811709, 0.047093212384734885, 0.04305275926927226, + 0.03996295545125314, 0.037342470160122807, 0.034882136983399127, 0.032503659051358803, + 0.030283400136758603, 0.02833413318888715, 0.026716951960401861, 0.025410942079986124, + 0.024331796916649818, 0.02337377299542151, 0.022449863415678752, 0.021514957496811402, + 0.020568392490579755, 0.01964085763034322, 0.018774394653097468, 0.018003835097688575, + 0.017345096384933451, 0.016792115618059419, 0.016321161065989184, 0.015899496666583322, + 0.015494963584338724, 0.015083672032454636, 0.01465418542791474, 0.014207860128341418, + 0.013756034513597803, 0.013315363470225371, 0.012902744182252629, 0.012531062632620658, + 0.012206550587622779, 0.011928032071734655, 0.011687882396065772, 0.011474202968136603, + 0.011273564084990481, 0.011073676681398503, 0.010865484803782374, 0.010644371625452833, + 0.010410390548108059, 0.010167625776320365, 0.0099229246756533628, 0.0096843142934278093, + 0.0094594181781906695, 0.0092541394617191974, 0.0090717907524229518, 0.0089127512679708417, + 0.0087746355098170924, 0.0086528799882045482, 0.0085416037366272436, 0.0084345775301450585, + 0.0083261437478460067, 0.0082119579128977598, 0.0080894662035977619, 0.0079580821758342517, + 0.0078190728090770014, + }}, + {0.1, 0.1, { + 0.0090062792310367599, 0.0089763103229537327, 0.0089464212611438977, 0.0089166118745999529, + 0.0088868819926185005, 0.0088572314447995844, 0.0088276600610462386, 0.0087981676715640101, + 0.008768754106860524, 0.0087394191977450044, 0.0087101627753278302, 0.0086809846710200738, + 0.008651884716533045, 0.0086228627438778328, 0.0085939185853648598, 0.0085650520736034168, + 0.0085362630415012179, 0.0085075513222639381, 0.0084789167493947697, 0.0084503591566939696, + 0.0084218783782583993, 0.0083934742484810827, 0.0083651466020507551, 0.0083368952739514106, + 0.008308720099461854, 0.0082806209141552539, 0.0082525975538986929, 0.0082246498548527221, + 0.0081967776534709168, 0.0081689807864994239, 0.0081412590909765228, 0.0081136124042321765, + 0.0080860405638875903, 0.0080585434078547669, 0.008031120774336065, 0.008003772501823752, + 0.0079764984290995669, 0.0079492983952342798, 0.0079221722395872475, 0.0078951198018059731, + 0.0078681409218256705, 0.0078412354398688221, 0.0078144031964447452, 0.0077876440323491457, + 0.0077609577886636888, 0.0077343443067555582, 0.0077078034282770224, 0.0076813349951649987, + 0.0076549388496406148, 0.0076286148342087825, 0.0076023627916577538, 0.0075761825650586968, + 0.0075500739977652576, 0.0075240369334131293, 0.0074980712159196244, 0.0074721766894832375, + 0.0074463531985832195, 0.0074206005879791474, 0.0073949187027104943, 0.0073693073880962004, + 0.007343766489734246, 0.0073182958535012244, 0.0072928953255519133, 0.0072675647523188511, + 0.0072423039805119081, + }}, + {1, 1, { + 0.30578343743121617, 0.22353090785207047, 0.16051446063675165, 0.1129515247559839, + 0.077692805295879158, 0.052130471820707389, 0.034118052808226947, 0.021900706193489234, + 0.014054672192720866, 0.0094348380428787727, 0.0071294561658082839, 0.0064211585445520742, + 0.0067535017003224644, 0.0077023594259370345, 0.0089515551414149543, 0.010272193103561093, + 0.011505208379362811, 0.012546710084893093, 0.01333574144957218, 0.013844124296688954, + 0.014068094999305518, 0.014021474301280706, 0.013730144975322444, 0.01322763947964649, + 0.012551664897273601, 0.011741414794309514, 0.010835537486808294, 0.0098706478071970208, + 0.0088802850359291433, 0.0078942334170194074, 0.0069381337938213914, 0.0060333255532312949, + 0.0051968674062068659, 0.0044416936995249363, 0.0037768700743938586, 0.0032079184762107181, + 0.0027371868797487651, 0.0023642437186837954, 0.0020862809817610906, 0.0018985133358478568, + 0.0017945635267689855, 0.0017668267534037942, 0.0018068087639292491, 0.0019054341344902294, + 0.0020533226039047238, 0.0022410324924863299, 0.0024592711636167728, 0.0026990732243764993, + 0.0029519477338836629, 0.0032099961193763488, 0.0034660028120197025, 0.0037135008258777657, + 0.0039468146310966679, 0.0041610827306636806, 0.0043522623518378357, 0.0045171186195308681, + 0.0046532004991045438, 0.0047588056884702796, 0.0048329365110794751, 0.004875248718383512, + 0.0048859949577008321, 0.0048659645034183257, 0.0048164206896196965, 0.0047390373234881171, + 0.0046358352035360088, + }}, + {3, 1.000001, { + 0.74475117928399825, 0.19301480699478943, 0.096424859064270632, 0.068570314812382363, + 0.050335251975315243, 0.038575123355758877, 0.032018382289449539, 0.027957865424385699, + 0.024622415096584287, 0.021581727961463334, 0.019017764465842378, 0.017083912238223874, + 0.015702562949911389, 0.014653851764227066, 0.013728917355293153, 0.012819003679224939, + 0.011918990549527974, 0.0110804082397495, 0.010357480024085333, 0.0097733886020108934, + 0.0093133543820348736, 0.0089371694493338236, 0.0085989756523870507, 0.0082639769621387915, + 0.0079167235552967986, 0.0075605312755838063, 0.0072108232643770505, 0.0068863105486156534, + 0.0066014068014224559, 0.0063618725366416558, 0.0061641393105867819, 0.0059975732240861454, + 0.0058483155997312167, 0.0057032714852204751, 0.0055531467143859464, 0.0053939551125154576, + 0.0052269420220770581, 0.0050572699000465734, 0.0048920268531955184, 0.004738150641274033, + 0.0046007509261207719, 0.0044821224929189919, 0.004381533130242499, 0.0042956904272023599, + 0.0042196715271262796, 0.0041480488969140002, 0.0040759573256559673, 0.0039999059484499668, + 0.003918222736415118, 0.0038311067192148634, 0.0037403384585543479, 0.0036487512829798526, + 0.0035595902797163704, 0.003475884394345762, 0.0033999346614104332, 0.0033329862644148026, + 0.0032751119803810656, 0.0032252968676545735, 0.0031816842274141732, 0.0031419240530708176, + 0.0031035582516268593, 0.0030643808239227948, 0.0030227234933132108, 0.0029776347878496673, + 0.0029289399679381933, + }}, + {1, 5, { + 0.015323581325245571, 0.011709912761402501, 0.0088389237715927015, 0.0065789165404490082, + 0.0048185222236491777, 0.0034639376189624141, 0.0024364963261597427, 0.0016705379774350413, + 0.0011115427180641362, 0.00071450139004990176, 0.00044249484649506858, 0.00026545852636324942, + 0.00015911087115166192, 0.0001040263880320796, 8.4836177782035159e-05, 8.9540568346182996e-05, + 0.00010892014272833827, 0.00013603293838391364, 0.00016578693839959662, 0.00019457818541612778, + 0.00021998593932444848, 0.00024051728012729556, 0.00025539443799149663, 0.00026437892257500549, + 0.00026762723158511056, 0.00026557355189097927, 0.0002588354324030241, 0.00024813891277010566, + 0.00023426004160550569, 0.00021798011779541575, 0.00020005234336081685, 0.00018117789079401717, + 0.0001619896658360851, 0.00014304229199689421, 0.00012480705910591165, 0.00010767076787172304, + 9.1937568591750612e-05, 7.7833037302065501e-05, 6.5509859065027582e-05, 5.5054597817256674e-05, + 4.6495127101404103e-05, 3.980837776043646e-05, 3.4928128795343045e-05, 3.1752627438144033e-05, + 3.015187529611744e-05, 2.9974460279393279e-05, 3.1053849917710194e-05, 3.3214091484923198e-05, + 3.627488886956309e-05, 4.0056046058942939e-05, 4.4381283068362072e-05, 4.9081442701414562e-05, + 5.3997116164695239e-05, 5.8980722715780594e-05, 6.3898083581373028e-05, 6.8629533680856543e-05, + 7.307061652561614e-05, 7.7132408295455503e-05, 8.0741516746069361e-05, 8.3839799471700289e-05, + 8.6383844304115502e-05, 8.8344252418432734e-05, 8.970476216257709e-05, 9.0461248836040776e-05, + 9.06206327043874e-05, + }}, + {3, 15, { + 0.024825552765284854, 0.0086055981295167245, 0.0026137569867761669, 0.00062996533088244321, + 8.0170526500387667e-05, -2.3000182738204848e-05, -1.895628984519167e-05, -4.746993166210417e-06, + 1.7301975625923541e-06, 2.4866043932213371e-06, 1.3561028457293986e-06, 3.2333863933313845e-07, + -1.3329857542613918e-07, -1.7419671599229458e-07, -6.0029855132187398e-08, 4.558595549065798e-08, + 8.9388294665288953e-08, 8.1867466710259292e-08, 5.2429114354357394e-08, 2.4571652538752982e-08, + 8.7762418797338745e-09, 5.0973419060463973e-09, 8.6112686155830899e-09, 1.3919756746787893e-08, + 1.7493445333155809e-08, 1.812890696022463e-08, 1.6325235445588284e-08, 1.3340942295716348e-08, + 1.0412818111822095e-08, 8.3258237201989098e-09, 7.317736412602483e-09, 7.2030170275981545e-09, + 7.5840301486243498e-09, 8.0486648283151522e-09, 8.3001320726679266e-09, 8.2066664591920535e-09, + 7.7865049969848852e-09, 7.1555724310782889e-09, 6.4651706336556706e-09, 5.849712483578083e-09, + 5.3948676941648352e-09, 5.1277561928641797e-09, 5.0248079884205049e-09, 5.0300390362036111e-09, + 5.0763053780493083e-09, 5.1036910785165806e-09, 5.0715895179967044e-09, 4.9634511354848375e-09, + 4.7850694269876039e-09, 4.5584232511112165e-09, 4.313483923143874e-09, 4.0801809972910155e-09, + 3.8821241799174158e-09, 3.7329287733689945e-09, 3.635276698562842e-09, 3.5822901814422457e-09, + 3.560460020792999e-09, 3.5532575065299183e-09, 3.5446306595168017e-09, 3.5217816600226331e-09, + 3.4768778516367797e-09, 3.4076054105567484e-09, 3.3166893418567223e-09, 3.2106504876143373e-09, + 3.0981415256036535e-09, + }}, + {5, 24, { + 0.10072555410123596, 0.012409352704075571, 0.00056296622236436567, -0.00015621338307643143, + -2.3011363528942111e-05, 8.8158251003798074e-06, 3.3918660059148895e-06, -2.2269502490613155e-07, + -4.2119851481707799e-07, -5.7834728690593434e-08, 5.3223046137289364e-08, 2.6283198091137186e-08, + -1.4709795326070311e-09, -5.769768158061231e-09, -1.6754726730322385e-09, 8.1320275161801728e-10, + 8.253236618078459e-10, 1.7398300047048713e-10, -1.4125333219119302e-10, -1.121653282743051e-10, + -3.7060111397033879e-12, 4.3671491521497327e-11, 3.2465383541295835e-11, 8.0263481442827883e-12, + -3.7236229361545075e-12, -2.4904983782810129e-12, 2.5606688617776819e-12, 5.1987581686731589e-12, + 4.6646747899654135e-12, 2.8994403262862976e-12, 1.6269699809224784e-12, 1.3099848543230559e-12, + 1.5428398220789131e-12, 1.7882832533547392e-12, 1.7921333457284085e-12, 1.5950611869142371e-12, + 1.3543823484446354e-12, 1.1865103205295047e-12, 1.114330424704578e-12, 1.0970250513705906e-12, + 1.0833989608709318e-12, 1.046170513891622e-12, 9.8669244467746728e-13, 9.2130901776895499e-13, + 8.6539310117267845e-13, 8.2488630224721649e-13, 7.9662500788927483e-13, 7.734816450886559e-13, + 7.4956176332130987e-13, 7.2271457147648854e-13, 6.9414571143152456e-13, 6.664738865423048e-13, + 6.4180821356092562e-13, 6.2078938852127571e-13, 6.0270376671156828e-13, 5.8622327636063605e-13, + 5.7019423379922266e-13, 5.5408598827663063e-13, 5.3800675339030525e-13, 5.2242392758028607e-13, + 5.0781211420730335e-13, 4.9440899690540852e-13, 4.8215237060000173e-13, 4.7076938444299488e-13, + 4.5993314998739351e-13, + }}, + {5, 24.000001, { + 0.10072554169188352, 0.012409352141109269, 0.00056296637857773733, -0.00015621336006506348, + -2.3011372344765525e-05, 8.8158217085136873e-06, 3.391866228609704e-06, -2.2269460370764522e-07, + -4.2119845698232263e-07, -5.7834781913626487e-08, 5.3223019854090514e-08, 2.6283199562113836e-08, + -1.470973762839705e-09, -5.7697664825881498e-09, -1.6754734862345781e-09, 8.1320192629444164e-10, + 8.2532348782477466e-10, 1.7398314172376339e-10, -1.4125322002586648e-10, -1.1216532456827213e-10, + -3.7060548111787214e-12, 4.3671459056117766e-11, 3.2465375514945824e-11, 8.0263518679044822e-12, + -3.7236204456548466e-12, -2.4905009389472778e-12, 2.5606636630218405e-12, 5.1987535039998139e-12, + 4.6646718905258978e-12, 2.8994386993169698e-12, 1.6269686709383943e-12, 1.3099833114841264e-12, + 1.5428380337965539e-12, 1.7882814612221893e-12, 1.7921317506678972e-12, 1.5950598325324805e-12, + 1.3543811619348708e-12, 1.1865092061996275e-12, 1.1143293276800672e-12, 1.0970239679721517e-12, + 1.0833979147009101e-12, 1.0461695271996371e-12, 9.866915233688812e-13, 9.2130815237626552e-13, + 8.6539227628677366e-13, 8.2488550562259458e-13, 7.9662423440800373e-13, 7.7348089552725315e-13, + 7.4956104060708464e-13, 7.2271387733110965e-13, 6.9414504495795824e-13, 6.6647324473440101e-13, + 6.4180759277183781e-13, 6.2078878581780143e-13, 6.0270318048857637e-13, 5.8622270616667867e-13, + 5.7019367971350274e-13, 5.540854502701379e-13, 5.3800623096663106e-13, 5.2242341976841853e-13, + 5.0781161979854701e-13, 4.9440851475327277e-13, 4.8215189983084682e-13, 4.7076892451006907e-13, + 4.5993270058770454e-13, + }}, + {5, 24.001, { + 0.10071314503004102, 0.01240878965975035, 0.00056312242424029133, -0.00015619036730555525, + -2.30201776580724e-05, 8.8124331230965733e-06, 3.39208849035018e-06, -2.222738553175477e-07, + -4.2114065348124485e-07, -5.7887938594886818e-08, 5.3196762204669987e-08, 2.6284666186064994e-08, + -1.4652106023207989e-09, -5.7680922789243673e-09, -1.6762854631510296e-09, 8.123775149712462e-10, + 8.2514960819940318e-10, 1.7412419772063383e-10, -1.4114116872320139e-10, -1.121616004328302e-10, + -3.7496663998710475e-12, 4.3639030151750568e-11, 3.2457355331755268e-11, 8.0300705215431959e-12, + -3.7211311583080581e-12, -2.4930564485410342e-12, 2.5554724354632253e-12, 5.1940949433322439e-12, + 4.6617761629058482e-12, 2.8978140110407347e-12, 1.6256607671900923e-12, 1.3084429083439795e-12, + 1.541052434626442e-12, 1.7864919173139212e-12, 1.7905389615349613e-12, 1.593707397635275e-12, + 1.3531963951065244e-12, 1.1853965384368016e-12, 1.1132339411783659e-12, 1.095942175330565e-12, + 1.0823532835497459e-12, 1.0451842819572544e-12, 9.8577156821879984e-13, 9.204440369781938e-13, + 8.645686130540517e-13, 8.2409006385525194e-13, 7.9585190090464344e-13, 7.7273244456695625e-13, + 7.4883939571163592e-13, 7.2202075889505735e-13, 6.9347955834564849e-13, 6.6583238862301114e-13, + 6.4118772542660691e-13, 6.2018697777119225e-13, 6.0211782843999862e-13, 5.8565335908018428e-13, + 5.6964041672727251e-13, 5.5354824265058871e-13, 5.3748458328639999e-13, 5.2191636259023749e-13, + 5.0731794620814019e-13, 4.9392707984286332e-13, 4.8168183110725144e-13, 4.7030967591979965e-13, + 4.5948396954522283e-13, + }}, + {8, 64, { + 0.055389182840797378, 0.0043389779909437298, 0.00010169479666274366, -2.6545125130361575e-05, + -2.1669422678713361e-06, 5.3495293185651982e-07, 1.0032392095936986e-07, -1.3274109780739722e-08, + -5.7189094049510687e-09, 1.6902945968810223e-10, 3.6801664960312614e-10, 3.1959107824725838e-11, + -2.3859676020200577e-11, -6.0363004419889661e-12, 1.2496332105594129e-12, 7.9810047866650672e-13, + 7.6927849592875231e-15, -8.8699917269487178e-14, -1.8435453579080286e-14, 7.4587014165571689e-15, + 4.1203831570330227e-15, -1.0698857689069266e-16, -6.3949174160087049e-16, -1.4490977666988444e-16, + 6.6086041771890336e-17, 4.2794551478423329e-17, 5.1010250292708728e-19, -7.7255498628710451e-18, + -2.4692687770265948e-18, 7.4399527452710886e-19, 7.4987770664038065e-19, 9.6927000606964107e-20, + -1.2980528536093528e-19, -6.8586967750567336e-20, 4.6464584483181796e-21, 1.8125442553337923e-20, + 5.9636774795498607e-21, -2.1625238572504615e-21, -2.4323027019900292e-21, -4.5627883229748643e-22, + 4.6203901210643128e-22, 3.3163149488896649e-22, 1.9210373033293267e-23, -8.7417382748169941e-23, + -4.7139320613577721e-23, 3.6772637125966488e-24, 1.6401191798142489e-23, 7.1082222823021329e-24, + -1.6232068500284555e-24, -3.1814293115343314e-24, -1.1560864988255874e-24, 4.5361043360799693e-25, + 6.5129250048961442e-25, 2.07103658095163e-25, -1.1454597499046347e-25, -1.4207462550641311e-25, + -4.1718305732752036e-26, 2.8309117652578871e-26, 3.3237916002898881e-26, 9.6798707008087584e-27, + -6.9496778603681839e-27, -8.288733568387073e-27, -2.5540741707536651e-27, 1.7284237458332444e-27, + 2.2238405072811642e-27, + }}, + {8.000001, 64, { + 0.055389245340797331, 0.0043389784792210728, 0.00010169431982549426, -2.6545147571422341e-05, + -2.1669321404298799e-06, 5.349541768880977e-07, 1.0032361216030378e-07, -1.3274191251369284e-08, + -5.7189000345747676e-09, 1.6903540120386708e-10, 3.6801666853860837e-10, 3.1958656206618212e-11, + -2.3859738405842464e-11, -6.0362683149649828e-12, 1.2496444197777138e-12, 7.9809896003151928e-13, + 7.6912409337383887e-15, -8.8700028776164938e-14, -1.8435278540609537e-14, 7.4587538101165156e-15, + 4.1203699849791334e-15, -1.0699951226934033e-16, -6.3949237172258678e-16, -1.449081478818577e-16, + 6.6086586359216132e-17, 4.2794413812008783e-17, 5.0996190995426995e-19, -7.7255657552446085e-18, + -2.4692457905882517e-18, 7.4400595255742494e-19, 7.4987637981685232e-19, 9.6924203145469345e-20, + -1.2980597683847522e-19, -6.8586568920754843e-20, 4.6467738236341091e-21, 1.8125466127101485e-20, + 5.9636071036441225e-21, -2.1625586107558116e-21, -2.4322984898203181e-21, -4.5626756187372936e-22, + 4.6204277335107345e-22, 3.3162996435369192e-22, 1.9208573811621098e-23, -8.741776956975721e-23, + -4.7138957940846093e-23, 3.6775612344731395e-24, 1.6401225246410479e-23, 7.1081431220891479e-24, + -1.6232589441511102e-24, -3.1814302919147399e-24, -1.1560692258756466e-24, 4.5362024856732771e-25, + 6.5129186864428061e-25, 2.0709976784753812e-25, -1.1454799127766748e-25, -1.4207436392248916e-25, + -4.1717388377892555e-26, 2.8309573739608995e-26, 3.3237839957189753e-26, 9.6796430235260638e-27, + -6.9497919396266815e-27, -8.2887145559433986e-27, -2.5540147512207089e-27, 1.7284551919210294e-27, + 2.2238365001937889e-27, + }}, + {5, 32, { + 0.027645307962757153, 0.0055066584449269601, 0.0008473811954847855, 6.4539020166722513e-05, + -1.2695408753157689e-05, -5.4112836970132262e-06, -4.0064067693526899e-07, 3.3896177618252072e-07, + 1.362844444801762e-07, 4.5478852911193296e-09, -1.5613549831478602e-08, -6.1015868421838021e-09, + 1.7563854040510194e-10, 1.1221789456871636e-09, 4.3949925961484919e-10, -3.0901193446422056e-11, + -1.0859539640306354e-10, -4.5303684006691563e-11, 3.6689741693832831e-12, 1.3666497373249963e-11, + 6.523426165179652e-12, -1.6774853886183765e-13, -2.0340568812895196e-12, -1.1791723866479867e-12, + -7.5605783129107564e-14, 3.5376066973529932e-13, 2.6671187354832973e-13, 5.782487975458467e-14, + -5.72758023363757e-14, -6.2557902728975533e-14, -2.2520963171372971e-14, 9.6085613878234001e-15, + 1.7995949709257967e-14, 1.1050439766061909e-14, 1.7162138219867535e-15, -3.0210294600700822e-15, + -2.8623972870578468e-15, -5.7965290341774917e-16, 1.3123748755964364e-15, 1.8675605758064843e-15, + 1.4076397501851233e-15, 6.7291425173688845e-16, 1.9165194884643074e-16, 9.4500572104776868e-17, + 2.3976527600881866e-16, 4.2219714555259875e-16, 5.1286206216416351e-16, 4.9144488077416567e-16, + 4.0686730736554155e-16, 3.2067381530596695e-16, 2.6984723178490705e-16, 2.5862037809749753e-16, + 2.7016773325840794e-16, 2.8361880685614872e-16, 2.858565681035134e-16, 2.7470731727441504e-16, + 2.5564140919800899e-16, 2.3620646102229406e-16, 2.2168143529809372e-16, 2.1344751439481705e-16, + 2.0974731001847857e-16, 2.075755033907829e-16, 2.044392323076919e-16, 1.9926774176968305e-16, + 1.9238779527914303e-16, + }}, + {5, 32.000001, { + 0.027645302456099145, 0.0055066575975457992, 0.00084738113094575914, 6.4539032862128519e-05, + -1.2695403341874207e-05, -5.4112832963723808e-06, -4.0064101589697618e-07, 3.3896163989807889e-07, + 1.3628443993228312e-07, 4.5479009046660707e-09, -1.5613543729891689e-08, -6.1015870178217805e-09, + 1.7563741822637883e-10, 1.1221785061878897e-09, 4.3949929051598826e-10, -3.0901084851048583e-11, + -1.0859535109937781e-10, -4.5303687675658893e-11, 3.6689605028892058e-12, 1.3666490849823731e-11, + 6.5234263329271732e-12, -1.677465048055511e-13, -2.0340557021171739e-12, -1.1791723110420268e-12, + -7.5606136889643043e-14, 3.5376040302345533e-13, 2.6671181572342149e-13, 5.7824937030355583e-14, + -5.7275739778484393e-14, -6.2557880208007604e-14, -2.2520972779925339e-14, 9.6085433918792607e-15, + 1.7995938658819089e-14, 1.1050438049846582e-14, 1.7162168430147746e-15, -3.0210265976730922e-15, + -2.8623967074042886e-15, -5.7965421579168768e-16, 1.3123730080365691e-15, 1.8675591681670743e-15, + 1.4076390772709691e-15, 6.7291406008498732e-16, 1.9165185434597875e-16, 9.4500332339712557e-17, + 2.397648538119306e-16, 4.221966326907836e-16, 5.1286157071948742e-16, 4.9144447390701964e-16, + 4.0686698669186199e-16, 3.2067354545886515e-16, 2.6984697316466466e-16, 2.5862010792990674e-16, + 2.701674496397447e-16, 2.8361852099971867e-16, 2.8585629339632467e-16, 2.7470706163312456e-16, + 2.5564117299165937e-16, 2.3620623934096606e-16, 2.2168122185068473e-16, 2.1344730464761136e-16, + 2.0974710244307793e-16, 2.0757529895165074e-16, 2.0443903304004683e-16, 1.9926754938198071e-16, + 1.9238761037039748e-16, + }}, + {6, 36.000001, { + 0.073852235817559966, 0.0079701697205169005, 0.00033209049224741346, -7.3388771487195562e-05, + -1.1153798600177903e-05, 2.3454686517231338e-06, 8.8936602136609653e-07, -4.0494405145732222e-08, + -7.6786850114945888e-08, -8.8497413194098428e-09, 6.5913353472540105e-09, 2.4409696226520867e-09, + -3.1362861919441577e-10, -4.3018856555824285e-10, -6.4740565609355443e-11, 5.4813602021949807e-11, + 2.8598241889882776e-11, -1.3622333135283962e-12, -6.3620922550928813e-12, -1.9973626962325359e-12, + 7.1000309457297723e-13, 7.6930704036749994e-13, 1.4216967133933919e-13, -1.4148819785234376e-13, + -1.0313199319805671e-13, -7.6201675396057555e-15, 2.733063706347003e-14, 1.637892483027302e-14, + -7.4052073616096013e-17, -5.4009073021146672e-15, -2.9391954497002736e-15, 2.5853087214934302e-16, + 1.2600478591261999e-15, 7.037942928353755e-16, -1.6815367274668533e-17, -2.666943373523852e-16, + -1.5086540633157496e-16, 2.434433846370766e-17, 9.4836086348222002e-17, 6.9092724363270385e-17, + 1.9147592118533829e-17, -7.122930515608575e-18, -5.2359042581644136e-18, 7.475456971078103e-18, + 1.5847680612603873e-17, 1.5901757975811551e-17, 1.1425795223045836e-17, 7.3104078253282104e-18, + 5.7495073795305701e-18, 6.1983014505083363e-18, 7.1109564160246186e-18, 7.4247823586332069e-18, + 6.9920117192347971e-18, 6.2265892718439936e-18, 5.5765330053504795e-18, 5.2295799090046812e-18, + 5.1166430566280881e-18, 5.0685701113304228e-18, 4.9609624082419095e-18, 4.7673689397869361e-18, + 4.5316147538192538e-18, 4.3115386286819162e-18, 4.1389670394904547e-18, 4.0115477888324276e-18, + 3.9073702421200368e-18, + }}, + {10, 40, { + 0.14012476628594236, 0.0017515680712017309, 6.5678867274483488e-05, 4.1076295043814153e-06, + 3.5803712686362272e-07, 4.0931884900705807e-08, 5.3468266172469953e-09, 9.7691335064855129e-10, + 1.4760714733222201e-10, 4.0577161509274016e-11, 8.2406408579642278e-12, 1.9773086740497482e-12, + 7.7620892557488272e-13, 1.7991555764472048e-13, 6.1048385040411031e-14, 2.9604260106596645e-14, + 8.8553769016243945e-15, 3.4723868529185512e-15, 1.9617937300087705e-15, 8.0472078599305464e-16, + 3.4780441639372197e-16, 2.0770708231240679e-16, 1.1239409224585654e-16, 5.7354630489886901e-17, + 3.4734865457940712e-17, 2.2210219076350677e-17, 1.3688500749748412e-17, 8.8644325208838621e-18, + 6.1611821561112725e-18, 4.3565627341233443e-18, 3.1276494729031074e-18, 2.3300827795811149e-18, + 1.7904962808876594e-18, 1.4014607194185911e-18, 1.1182084052373084e-18, 9.1152312409628183e-19, + 7.5657859015165825e-19, 6.3720950184893892e-19, 5.4406039606385617e-19, 4.7057125325843754e-19, + 4.1165642511477293e-19, 3.6369486508526449e-19, 3.2420444974463385e-19, 2.9136205383836536e-19, + 2.6376043816805387e-19, 2.4033018224897548e-19, 2.2026810031455024e-19, 2.0295642943322802e-19, + 1.8790738542548335e-19, 1.7473375428430525e-19, 1.6312825588511401e-19, 1.5284474127602852e-19, + 1.436832645476049e-19, 1.3547998045038396e-19, 1.2810003586553726e-19, 1.2143181521288202e-19, + 1.1538220589495714e-19, 1.0987295964292659e-19, 1.0483796160571163e-19, 1.002210836253423e-19, + 9.5974426842541742e-20, 9.2056893541104833e-20, 8.8433049238610855e-20, 8.507221192335366e-20, + 8.1947706726315484e-20, + }}, + {20, 400, { + 0.02215567313631895, 0.00065269459142039865, 2.4476047178264951e-06, -7.6595247051358436e-07, + -9.631771616993863e-09, 2.8030195224338189e-09, 7.500239489830876e-11, -1.6612908546891741e-11, + -7.8495889650223648e-13, 1.3284575718033505e-13, 1.0119437302933083e-14, -1.2949614227517681e-15, + -1.5325082020574843e-16, 1.4213353093765645e-17, 2.6406585845180585e-18, -1.5935026027143987e-19, + -5.0515964749497694e-20, 1.4341248834671887e-21, 1.0506108405534542e-21, 5.6003086547007695e-24, + -2.3281198397874052e-23, -1.0140383204328361e-24, 5.3752090101276467e-25, 4.8780891159569323e-26, + -1.2545372225148782e-26, -1.9595547869073471e-27, 2.8071164146330743e-28, 7.4460753776151095e-29, + -5.2504232367923821e-30, -2.7662154221754204e-30, 3.0939155553025526e-32, 1.0100954557476259e-31, + 4.7408624903723565e-33, -3.5869359286369608e-33, -4.0993337431195713e-34, 1.2010414412543398e-34, + 2.5109025556979659e-35, -3.4954575214561228e-36, -1.3585227461304562e-36, 6.2844192569403284e-38, + 6.8285272504740454e-38, 2.0932177365703496e-39, -3.2070898913585012e-39, -3.5788231131558673e-40, + 1.371807463900023e-40, 3.0244000227878744e-41, -4.8881886902630141e-42, -2.0955390640020594e-42, + 9.2709205121118948e-44, 1.2946975724591333e-43, 6.3705997628129505e-45, -7.2046487507044771e-45, + -1.0769262040133407e-45, 3.4698053859350254e-46, 1.0306072082361817e-46, -1.220629821147018e-47, + -8.081253230046281e-48, -1.9163077142841656e-50, 5.5303266332597602e-49, 6.1041796805653616e-50, + -3.2614776359789437e-50, -8.0529860813504264e-51, 1.483532250442276e-51, 7.7250236600574824e-52, + -2.2145978617607986e-53, + }}, + {100, 10000, { + 0.0044311346272637902, 2.5221556731363191e-05, 3.7832335097044787e-09, -1.2490541916225739e-09, + -6.2466896706790087e-13, 1.8717202396481944e-13, 1.9667820693065069e-16, -4.6697947798307578e-17, + -8.4212683659277555e-20, 1.6290392643984634e-20, 4.4966533377188696e-23, -7.2950317476810499e-24, + -2.8639275331847727e-26, 3.9852266847868847e-27, 2.1146999923370076e-29, -2.5672401172709887e-30, + -1.7741814607089997e-32, 1.9033311661182e-33, 1.6657694136257103e-35, -1.5945793993913701e-36, + -1.7295050630227009e-38, 1.4881123043217804e-39, 1.9666699993001116e-41, -1.5291635607811572e-42, + -2.4298461966199657e-44, 1.7136869486751306e-45, 3.2399876499868141e-47, -2.0774517390305975e-48, + -4.6358181304776465e-50, 2.7051107291433755e-51, 7.0821023559211063e-53, -3.7598928599384549e-54, + -1.1501457860662292e-55, 5.5467165943648605e-57, 1.9779665133443172e-58, -8.6390974770700275e-60, + -3.5897229122015068e-61, 1.4133873333425913e-62, 6.8537695794275029e-64, -2.4165732881213056e-65, + -1.3727807354178211e-66, 4.2949018075151605e-68, 2.8771247638855043e-69, -7.8869093913345812e-71, + -6.2946949819778018e-72, 1.4855660904105821e-73, 1.4344912675786827e-74, -2.8423926363386674e-76, + -3.3981405042499157e-77, 5.4444993282579962e-79, 8.3516387627674988e-80, -1.0179832782881088e-81, + -2.1256904247760159e-82, 1.7603565845605231e-84, 5.5933580523239532e-85, -2.383979982146934e-87, + -1.5190165601916171e-87, 1.3783687115981443e-91, 4.2506997850282832e-90, 1.819122145504535e-92, + -1.2236729414721951e-92, -1.0895465507620822e-94, 3.6180420722116635e-95, 5.0106914011612727e-97, + -1.0968972155154744e-97, + }}, +}; +// clang-format on diff --git a/export/tests/unit/truncated-coulomb-vrr-reference.h b/export/tests/unit/truncated-coulomb-vrr-reference.h new file mode 100644 index 000000000..490d3dfd8 --- /dev/null +++ b/export/tests/unit/truncated-coulomb-vrr-reference.h @@ -0,0 +1,684 @@ +// Generated by tools/generate_truncated_coulomb_reference.py --vrr-output +// mpmath 1.4.1; decimal precision 100; mmax 16. +// G_m(R, T) for physical geometries across the cutoff transition; +// zeta, rho and W-P define the (L0|00) Obara-Saika recurrence. +// clang-format off +struct TruncatedCoulombVRRReference { + double zeta, rho, wp, r, t; + double g[17]; +}; +static const TruncatedCoulombVRRReference truncated_coulomb_vrr_reference[] = { + {0.0112, 0.011199672489995071, 17.938088173556768, 0.91541656904529334, 3.6039854696257025, { + 0.035018746558547219, 0.027501782223422509, 0.021374910603150221, 0.016419724182643074, + 0.012447325787716864, 0.0092948488656160056, 0.0068223407121181913, 0.0049099747646543469, + 0.0034555609398262883, 0.0023723256327073563, 0.001586935429910283, 0.0010377408331266997, + 0.000673218358255763, 0.00045059128002212649, 0.00033461104497243827, 0.0002964829879825896, + 0.00031292146921304247, + }}, + {0.14499999999999999, 0.14499983313511267, 2.0071569143838284, 0.82250302216781057, 0.58415910969484874, { + 0.30781308250732753, 0.24724430397201191, 0.19678488015820167, 0.15503974834390627, + 0.12077605428697882, 0.092906928494858085, 0.070476683159819753, 0.052647319252612922, + 0.03868624099435012, 0.02795508216302137, 0.019899555465512996, 0.014040242545505052, + 0.0099642481276182582, 0.0073176473428272128, 0.0057986604622857693, 0.0051514941079144451, + 0.0051607925287883646, + }}, + {9920, 8840.035118525022, 2.0443053897237324, 220.01021861494436, 46522.270406162497, { + 0.0041087951539969369, 4.415948581730991e-08, 1.4229724371916724e-12, 9.2932758586969676e-17, + -3.0427469925426849e-19, 5.6642068775739967e-21, -9.9752138356958289e-23, 1.6863889476801508e-24, + -2.7208142296214554e-26, 4.1567770592057315e-28, -5.9468365626026018e-30, 7.8309034584131307e-32, + -9.2133840843249903e-34, 9.104136573858645e-36, -6.2741882121171803e-38, -1.7262347144663126e-41, + 9.68895795078459e-42, + }}, + {24.730921310199538, 24.69604475591548, 3.8716670370604649, 12.246329291783365, 371.23523495816056, { + 7.1330272502272575e-25, 2.634763397507812e-25, 9.6726266589246385e-26, 3.5286124004784754e-26, + 1.278901504437877e-26, 4.6042139809748319e-27, 1.6461334736339572e-27, 5.8433953521994238e-28, + 2.0589635923734438e-28, 7.1994597658643328e-29, 2.4974393008807105e-29, 8.5921120042778571e-30, + 2.9306909512299634e-30, 9.9071310461517779e-31, 3.3178907456330934e-31, 1.1003346767366907e-31, + 3.6118491161631561e-32, + }}, + {107810.98505610731, 0.03771483168463792, 1.4637046315625918e-05, 12.492301679654465, 66.026875854371497, { + 0.10906478969231902, 0.00082591210700321287, 1.8763078530320738e-05, 7.1043887279774885e-07, + 3.765674249319856e-08, 2.5677331592748457e-09, 2.1332447122262322e-10, 2.1237480033157785e-11, + 2.3201143684512955e-12, 3.3170305517803748e-13, 3.7093975518469692e-14, 8.8640914762754514e-15, + 8.8407678947756675e-16, 2.5517579009622186e-16, 6.2249634323815307e-17, 2.961081793236907e-18, + 4.0694689229670935e-18, + }}, + {13.014253415611124, 12.999573398225996, 66.279604492826138, 236.5952207417518, 57235.994623261198, { + 3.3958544462562211e-07, 3.9968997642010422e-09, 4.4264851744601446e-11, 4.5569840766876586e-13, + 4.2811349821326098e-15, 3.5573620720693197e-17, 2.4542997751716164e-19, 1.1719763119131567e-21, + 1.3850308534797836e-25, -7.0497777840774629e-26, -7.9775523279922122e-28, -3.3554531521000664e-30, + 3.2470305897711216e-32, 6.8843061875224001e-34, 4.3093998585190208e-36, -3.0152089827918613e-38, + -8.6935288109740993e-40, + }}, + {0.046050947132961054, 0.046050942243674731, 156.91664484520962, 19.910104456403868, 1133.904920912677, { + 2.8964619643338446e-86, 1.1882591018108635e-86, 4.8672802425837672e-87, 1.9906192530615585e-87, + 8.1285180122720734e-88, 3.3139844222134163e-88, 1.3489628110606259e-88, 5.4821931409795653e-89, + 2.2243742146982273e-89, 9.0106219216756108e-90, 3.6440909828692682e-90, 1.4713134956174137e-90, + 5.9305738814812824e-91, 2.3864786053325014e-91, 9.586970980697693e-92, 3.8446953541948946e-92, + 1.5391913075501492e-92, + }}, + {22.131287134235002, 22.130293550806368, 6.4106292371011593, 45.523830163281723, 909.55191106843199, { + 0.029385372159625178, 1.6153763079397406e-05, 2.6640199777748659e-08, 7.3223417634445302e-11, + 2.8176727309550632e-13, 1.394041080558382e-15, 8.4296738314414275e-18, 6.0241619238648182e-20, + 4.9674146004390992e-22, 4.6421785925482601e-24, 4.8486178845366053e-26, 5.5973152459028807e-28, + 7.0770150163578946e-30, 9.7259635901988677e-32, 1.4435735538552021e-33, 2.3013328075263186e-35, + 3.9217836917913072e-37, + }}, + {3.9719192255035369, 3.0313997351418172, 7.6613705041382412, 14.685820952605123, 305.47138797197562, { + 1.9949098211212093e-06, 3.4034000645058655e-07, 5.5516010366400283e-08, 8.5919949964613877e-09, + 1.2479768415593792e-09, 1.6732643278158403e-10, 2.0139030673015689e-11, 2.0577721068891249e-12, + 1.5301293926617099e-13, 2.1999998059520793e-15, -1.7563011266817836e-15, -3.6042370065662181e-16, + -3.6720744787722292e-17, 5.0993889691313022e-21, 7.78687834776506e-19, 1.5210604988342819e-19, + 1.088478948371713e-20, + }}, + {0.6347663879312937, 0.62749966908641375, 2.2099772552729502, 16.145522266628806, 3.1361003496227196, { + 0.49429993402270878, 0.071880355942840243, 0.027452749181514618, 0.014956764808696852, + 0.0097645852519086043, 0.0070835348200046092, 0.0054951953001411528, 0.0044618507965133525, + 0.0037428404889067696, 0.0032167932856209838, 0.0028167397788959429, 0.0025030481155822052, + 0.002250915128413231, 0.0020440932312414777, 0.0018715284459754245, 0.0017254565400903314, + 0.0016002736532318975, + }}, + {492073.98955975013, 3.4184405614128002, 6.1774887631898073e-05, 19.653328171402912, 270.30663988431166, { + 0.053903303113252289, 9.9738213673480222e-05, 5.4750939497824262e-07, 6.1614162025985532e-09, + -1.0868500343622925e-10, 2.794865671907108e-11, -3.6568304411352711e-12, 4.3058878512635009e-13, + -3.8377356392871051e-14, 1.5798338905453184e-15, 2.4715186912240634e-16, -6.5311113182956678e-17, + 6.6306132624572185e-18, 1.2816955377340915e-19, -1.5376819068058946e-19, 2.1346671119919435e-20, + 3.0870741378560438e-22, + }}, + {287.99394733161313, 1.6678101226300761, 0.010880183769951441, 16.567973390103145, 5.8869802262400199, { + 0.3650376959548865, 0.03076809646572715, 0.007603978694271824, 0.0029934318740310628, + 0.0015439730940377118, 0.0009444921207182117, 0.00064668710153280176, 0.00047830866880617238, + 0.00037364530656647092, 0.00030377418513856624, 0.00025449079435946591, 0.00021819010880773661, + 0.00019050747981587618, 0.00016879129449208564, 0.00015135263866006743, 0.00013707210571425307, + 0.00012518220302899113, + }}, + {179305.1615856381, 3776.4954028499606, 0.029098942112227112, 88.76110865652808, 7208.6016010399562, { + 0.010438051087083089, 7.2401166677352281e-07, 1.501129915247211e-10, 7.5855940304160434e-14, + -9.6579269834435704e-16, 3.979973200713671e-17, -1.5166737776132083e-18, 5.4424269873154147e-20, + -1.8142583630462134e-21, 5.5028473028494902e-23, -1.462737639039628e-24, 3.1308016000645842e-26, + -3.937603807530197e-28, -5.9686550560976104e-30, 5.8347900844057897e-31, -2.0387157387861763e-32, + 3.4012525974303702e-34, + }}, + {0.61382450907256159, 0.6136300660825591, 31.827709424632271, 18.755122037972178, 622.0031564652088, { + 3.9022548370215755e-20, 9.8319527984652642e-21, 2.4541575549475982e-21, 6.0666607126027647e-22, + 1.4846194849934474e-22, 3.5951596398675605e-23, 8.6111676516651027e-24, 2.0390673654069417e-24, + 4.7707905666236596e-25, 1.1022429488398595e-25, 2.5130588373858056e-26, 5.649922760678324e-27, + 1.2515050061101675e-27, 2.7287253941805557e-28, 5.8499434947544633e-29, 1.231577934975111e-29, + 2.5424454256058541e-30, + }}, + {0.24841272195340627, 0.010071704849717988, 3.5464310196600137, 1.6083322374642537, 77.059770411980892, { + 1.8537259734597591e-25, 1.5405751373650349e-25, 1.2784818557095894e-25, 1.0594104855521189e-25, + 8.7654637123144476e-26, 7.2411669052215453e-26, 5.9723669168020748e-26, 4.917778444542466e-26, + 4.0425488831201013e-26, 3.3172916241852833e-26, 2.7172595159231726e-26, 2.2216388697317895e-26, + 1.8129470510438741e-26, 1.4765189952500973e-26, 1.2000699914472953e-26, 9.7332381462195042e-27, + 7.8769679448179284e-27, + }}, + {28.5843159165488, 12.322668230632033, 1.5115145964742114, 19.074939324150794, 151.48706143133148, { + 0.072004090668502552, 0.00023765755962314227, 2.3532461192819911e-06, 3.8835760906694812e-08, + 8.9727242636524236e-10, 2.6653932556967207e-11, 9.6771716130776477e-13, 4.1522764314893499e-14, + 2.0557579596934949e-15, 1.1534940679741891e-16, 7.233748853564842e-18, 5.0139177394949099e-19, + 3.8062681908300889e-20, 3.1407580535715554e-21, 2.7989151599163098e-22, 2.6791374620244796e-23, + 2.740950894422078e-24, + }}, + {1236.5347567559706, 729.26167551458434, 0.67221194748561242, 27.188800651697786, 947.41819157882014, { + 5.463757046961374e-09, 6.6347190734186701e-10, 7.812638341342609e-11, 8.8912514519809277e-12, + 9.7385447180825645e-13, 1.0210198819810553e-13, 1.0172295291855883e-14, 9.5318977628606224e-16, + 8.270759010980518e-17, 6.4736607262043181e-18, 4.3400648897763837e-19, 2.1676840462601828e-20, + 3.0113661636403272e-22, -9.5100825170001787e-23, -1.4129902825317119e-23, -1.1201364991049994e-24, + -3.3938099896751e-26, + }}, + {0.091952358633389991, 0.090957113643486479, 31.874291710306323, 11.153167646253713, 94.443084530083382, { + 0.089258451226073787, 0.00081021119644210458, -3.6990492522705118e-05, 4.3313568818192404e-06, + -5.0342306650514746e-08, -5.4490735430018629e-08, 7.0458863108571691e-09, 6.484044026584734e-10, + -2.6529776732200938e-10, 1.7610473031727013e-13, 1.0326891736797964e-11, -7.1446684947371672e-13, + -4.5750987486083299e-13, 6.2433106668182694e-14, 2.3705886793596043e-14, -4.821224323147876e-15, + -1.453346430085756e-15, + }}, + {260126.03021814057, 1902.0104903721103, 0.023710480549422808, 134.53076223310649, 20000.250294042595, { + 6.0107598103852626e-25, 2.9607537505911618e-26, 1.4442579624282542e-27, 6.9746986515457446e-29, + 3.3335254874103624e-30, 1.5762642427483163e-31, 7.3712143413682715e-33, 3.4076807071742916e-34, + 1.5566795770530945e-35, 7.0235116869408795e-37, 3.1282377113236677e-38, 1.3746423468747113e-39, + 5.9559854239631183e-41, 2.5426817550111308e-42, 1.0687385413018811e-43, 4.4189427444781398e-45, + 1.7956048991539559e-46, + }}, + {0.46181700256458885, 0.21071933870830395, 42.401446702081344, 42.269681430837309, 1819.6874122158338, { + 0.006056941055405397, 0.00011983904763529117, 1.1739733222043773e-06, -2.0780494685731272e-08, + -8.5572191834846762e-10, 7.9928077815915415e-12, 1.0333199641492096e-12, 1.0072903238947867e-15, + -1.6976808053069364e-15, -2.3146738725768596e-17, 3.4467002943678313e-18, 1.0337743680191992e-19, + -8.0838971521315114e-21, -4.2800579525259619e-22, 2.0510363577464991e-23, 1.8414469376042911e-24, + -5.0804722559953096e-26, + }}, + {1779.9339402769649, 0.1988054085416493, 0.0042730327530860436, 4.0764507592180435, 290.97254787636047, { + 7.3243011176065893e-77, 5.6029885718371053e-77, 4.2833088552833875e-77, 3.2722208255694126e-77, + 2.4980825541826325e-77, 1.9057648217852165e-77, 1.4528724545923387e-77, 1.106824136341099e-77, + 8.4259705872179538e-78, 6.4098613636954074e-78, 4.8726123460287092e-78, 3.7013206723680085e-78, + 2.8095077976865166e-78, 2.1309803688836653e-78, 1.6151069409001363e-78, 1.2231863962355457e-78, + 9.2565759387686043e-79, + }}, + {70820.339970237721, 9.7585525916386207, 0.00028648882726085866, 8.1045575706939204, 42.183827924586545, { + 0.13489242188436576, 0.0020430340191726517, -3.7432224244318993e-05, 1.8493494457407182e-05, + -5.7602378444573728e-07, -2.9758885160757278e-07, 9.1797814984239393e-08, 5.1861151686677843e-09, + -5.9059752978183736e-09, 3.4119592967632571e-10, 4.6817134526989181e-10, -6.2763004547560815e-11, + -4.3001811020089683e-11, 9.0469432045823464e-12, 5.0464767935423648e-12, -1.2073814554499183e-12, + -7.0728031268276791e-13, + }}, + {0.38379677102492221, 0.015243859405177003, 4.2497952867494542, 1.2352830412832705, 174.51940338590455, { + 8.2423263850508982e-66, 7.5210963209864403e-66, 6.86105818907997e-66, 6.2571775657377559e-66, + 5.7048217186151163e-66, 5.1997285190878027e-66, 4.7379776887577316e-66, 4.3159642099563372e-66, + 3.9303737422312316e-66, 3.5781598979990171e-66, 3.2565232409740632e-66, 2.9628918806913078e-66, + 2.6949035454786862e-66, 2.4503890246466911e-66, 2.2273568784908253e-66, 2.0239793219866331e-66, + 1.8385791948333277e-66, + }}, + {422252.30123975646, 156204.34669011415, 4.4639469026471934, 4777.5325033454919, 22745158.928382773, { + 0.00018582342170901889, 4.0849004901244566e-12, 2.6939142322459701e-19, 2.9609753890138158e-26, + 4.5563163107276387e-33, 9.014412018018014e-40, 2.1797712726710342e-46, 6.2455941346667943e-53, + -2.526492984778114e-58, 4.5276166537908513e-61, -7.4398159104909018e-64, 1.2118941216732105e-66, + -1.9564303064246606e-69, 3.1293146719604403e-72, -4.9579231802766873e-75, 7.7783669735395009e-78, + -1.2080319931027461e-80, + }}, + {0.14354460108586714, 0.1435413848513212, 34.795873920362318, 19.458964815138085, 173.80092785386418, { + 0.067223155604473867, 0.00019339124489886682, 1.6690754815314394e-06, 2.4008437442506224e-08, + 4.8348148704490997e-10, 1.251815349122518e-11, 3.9614198274090692e-13, 1.481535764265587e-14, + 6.393244073862222e-16, 3.1267153991104484e-17, 1.7090633849956721e-18, 1.0325385339365411e-19, + 6.8310913600117939e-21, 4.9165630998815954e-22, 3.8064569926271747e-23, 3.2179708606306189e-24, + 2.7325798793983886e-25, + }}, + {3.3360608344064429, 3.3359664188507332, 24.246003321232141, 41.639549223693599, 1961.2211748102284, { + 1.8252953929961126e-06, 1.1647314428517485e-07, 7.0207011088258396e-09, 3.9531455284005288e-10, + 2.0446881627450351e-11, 9.4480361794720914e-13, 3.6939219523336307e-14, 1.0580020894328595e-15, + 8.0664591490481276e-18, -1.4343360084791811e-18, -1.0772903062957433e-19, -3.4755431233759804e-21, + 5.9915527432659153e-23, 1.389502872100702e-23, 7.126597028434619e-25, -8.9642566846253977e-28, + -2.6769997754117268e-27, + }}, + {102.760104798028, 9.0672534124734145, 8.9233175623266519, 299.65680446815446, 92731.230178138401, { + 9.0323567506876158e-15, 1.4717181004286922e-16, 2.3510051353551495e-18, 3.6775431846160757e-20, + 5.6251148560670027e-22, 8.399914778155029e-24, 1.2222802154101246e-25, 1.7292201602463715e-27, + 2.3721697739300416e-29, 3.145003025084413e-31, 4.0129619080202551e-33, 4.901404908578739e-35, + 5.6884361424536875e-37, 6.207410917599015e-39, 6.2664166613860738e-41, 5.6904874889809075e-43, + 4.386753519733656e-45, + }}, + {266562.85249662452, 2268.3687379131502, 0.058143726499563009, 332.83269381952266, 105898.87013402081, { + 0.0027233237743573419, 1.2858134231795046e-08, 1.8212849035389672e-13, 4.299585305381522e-18, + 1.4210301346734458e-22, 6.0384361030682057e-27, 3.1361426792973383e-31, 1.9249824703152717e-35, + 1.3548241988864913e-39, 2.8880084520821539e-43, -3.7476252010589448e-45, 7.7767956992301496e-47, + -1.5895499294041593e-48, 3.2070950087276316e-50, -6.3836725175136551e-52, 1.2528227286828021e-53, + -2.4225951316970522e-55, + }}, + {18.574541764113235, 2.8452627359219984, 2.5447493538882737, 16.807518220586008, 785.24262552768801, { + 1.8980470742879676e-58, 7.6381593209763044e-59, 3.0665959439684412e-59, 1.2282834402674448e-59, + 4.907978282731935e-60, 1.956393305233317e-60, 7.7794053423681309e-61, 3.085746244114477e-61, + 1.2209113746565945e-61, 4.8184171515795978e-62, 1.8967357173618471e-62, 7.446923479170362e-63, + 2.9160816960820632e-63, 1.1388325970739863e-63, 4.4354915723386936e-64, 1.7227771906111453e-64, + 6.6727570819504467e-65, + }}, + {3954.0191792015307, 0.011951897705152077, 0.00015933179442177405, 0.27271316514932759, 33.208217912559356, { + 5.7794335841821602e-16, 5.6562996123957302e-16, 5.535335500260254e-16, 5.4165094765370645e-16, + 5.2997901646791096e-16, 5.1851465785890801e-16, 5.0725481184176128e-16, 4.9619645664012584e-16, + 4.8533660827398549e-16, 4.7467232015129865e-16, 4.6420068266351953e-16, 4.539188227849622e-16, + 4.4382390367597425e-16, 4.3391312428988818e-16, 4.2418371898371902e-16, 4.1463295713257571e-16, + 4.0525814274775487e-16, + }}, + {9456.4065587906916, 1.4598608532863884, 0.0068007467191039698, 46.689286193647163, 2833.0485002492505, { + 1.9609676547784176e-22, 2.4394169219664954e-23, 3.0046828504501259e-24, 3.6631999757150088e-25, + 4.4188897521803595e-26, 5.2721144976945909e-27, 6.218562587403565e-28, 7.2481466817432656e-29, + 8.3440257730756678e-30, 9.4818880998619351e-31, 1.0629650635915409e-31, 1.1747736170507022e-32, + 1.2790074947205375e-33, 1.3705939106377314e-34, 1.4442651394012095e-35, 1.4949114738425709e-36, + 1.5179991096215319e-37, + }}, + {590.20781870475173, 0.059835946087753319, 0.0064087056797393193, 0.42320654485810449, 239.1048628649425, { + 6.2349006417534871e-102, 6.090641393160258e-102, 5.949473946616319e-102, 5.8113368395506174e-102, + 5.6761697446689303e-102, 5.5439134504295013e-102, 5.414509841832014e-102, 5.2879018815152032e-102, + 5.1640335911584647e-102, 5.0428500331828894e-102, 4.924297292747221e-102, 4.8083224600342837e-102, + 4.6948736128235028e-102, 4.5838997993452068e-102, 4.4753510214124302e-102, 4.3691782178260493e-102, + 4.2653332480490858e-102, + }}, + {3.7454032880753609, 3.7446739509163036, 19.921262766086294, 42.77688592657244, 1486.6779259800808, { + 0.022984575837767297, 7.7301830827898129e-06, 7.7991119018126077e-09, 1.3151164247067287e-11, + 2.7256784244384341e-14, 4.4776707060169627e-16, -3.2867326982571741e-17, 2.9665895848925146e-18, + -2.5012610713433936e-19, 1.972083796316296e-20, -1.4288718893902186e-21, 9.2469369361307841e-23, + -5.0456304798024409e-24, 1.9716644427888135e-25, -9.8045237067537393e-28, -7.1636482052949244e-28, + 7.7331824582110157e-29, + }}, + {46.470198218183391, 37.037600168683603, 1.1679198500802255, 19.837054047150126, 79.53023953650019, { + 0.099375378414713861, 0.00062476473724881591, 1.1783531790359096e-05, 3.7041042058445817e-07, + 1.6301176503443166e-08, 9.2235726552575054e-10, 6.3786617391783459e-11, 5.2132750443472134e-12, + 4.9163139782396175e-13, 5.254437690440748e-14, 6.2765004041359339e-15, 8.2865655412971797e-16, + 1.198229808941817e-16, 1.8832927826023531e-17, 3.1968283653242717e-18, 5.8284762585089275e-19, + 1.1359375067068224e-19, + }}, + {0.051115323541494041, 0.051113180434664979, 127.24777570088449, 20.738259369241458, 827.69383846614062, { + 1.0408304174698977e-31, 2.9340810442187683e-32, 8.2264715162106387e-33, 2.2937759352194211e-33, + 6.3595773439741423e-34, 1.753018831260577e-34, 4.803570734601179e-35, 1.3082660451273919e-35, + 3.5409004557884817e-36, 9.5223653748633454e-37, 2.5439760834681547e-37, 6.7505218926485045e-38, + 1.7788199115526254e-38, 4.6537782296116107e-39, 1.2085406845099517e-39, 3.1145479034848575e-40, + 7.963351551209881e-41, + }}, + {987.65217519296687, 708.31765308050831, 1.7864182939148088, 60.90462258913886, 4394.8726843793074, { + 1.6760157669987933e-16, 1.3870974263990746e-17, 1.130794893630105e-18, 9.0736093483429488e-20, + 7.1602762531855429e-21, 5.5516598455343183e-22, 4.2247034510764666e-23, 3.1515383810101273e-24, + 2.3014184778224976e-25, 1.6425177244742777e-26, 1.1434972485890142e-27, 7.747763974996901e-29, + 5.0947319649657227e-30, 3.2400933921057636e-31, 1.9840115285970261e-32, 1.1627862613121253e-33, + 6.4688305450065146e-35, + }}, + {61.646006290724806, 11.489013779787623, 10.620237325776104, 187.24671706067107, 37307.451742168727, { + 1.5602108326441497e-19, 4.8382345365151024e-21, 1.4803582571771605e-22, 4.46660706024756e-24, + 1.328161524408942e-25, 3.8894691370240485e-27, 1.1209040876016327e-28, 3.1762922324287893e-30, + 8.8417320472242998e-32, 2.4152226330150202e-33, 6.4662660640635081e-35, 1.6944194046430553e-36, + 4.3386470557355469e-38, 1.0834835685832291e-39, 2.6328644685869118e-41, 6.2080575974732706e-43, + 1.4154105003341615e-44, + }}, + {0.014691180814874355, 0.009379130983252949, 69.907514130086611, 0.097669338714730727, 112.46011549522025, { + 1.934292918458527e-51, 1.9285490870972481e-51, 1.9228191042130128e-51, 1.9171029425500986e-51, + 1.9114005748977499e-51, 1.9057119740901142e-51, 1.9000371130061786e-51, 1.8943759645697059e-51, + 1.8887285017491716e-51, 1.8830946975577007e-51, 1.8774745250530041e-51, 1.8718679573373151e-51, + 1.8662749675573272e-51, 1.8606955289041304e-51, 1.8551296146131481e-51, 1.849577197964075e-51, + 1.8440382522808135e-51, + }}, + {243164.91718078626, 92710.266285385092, 0.77053145050767935, 617.7292274845164, 378664.63375759724, { + 0.0014396098477127955, 4.279783276671967e-09, -9.1524799350197112e-12, 3.2130063069824261e-14, + -9.9485790194912464e-17, 2.5556915511178893e-19, -4.5789772919271182e-22, 7.5006302998797479e-26, + 3.3331586405839146e-27, -1.3485468490626797e-29, 1.6589176691952975e-32, 9.6352464750093516e-35, + -5.8788002964561491e-37, 8.5487488018009916e-40, 6.024687457902115e-42, -3.7698465499769817e-44, + 3.3009012702108435e-47, + }}, + {0.17289071431378139, 0.17288940155607313, 8.4515176740184632, 8.2148507299693563, 12.349357816075688, { + 0.25218686475748436, 0.010210350048372594, 0.0012400125535750575, 0.00025085233338121036, + 7.0920025503341753e-05, 2.5667244888324202e-05, 1.1255923677393129e-05, 5.7490748429960715e-06, + 3.3160966994628979e-06, 2.107045441693348e-06, 1.4454671742546368e-06, 1.0535917620170552e-06, + 8.0571134012539203e-07, 6.4012477416765312e-07, 5.24352640133497e-07, 4.4025313912057534e-07, + 3.7715787309474538e-07, + }}, + {6018.263488794335, 5723.7745632323158, 3.3465159654464602, 267.50942633054387, 70867.265945482868, { + 0.0032193431859157788, 6.7303503232881732e-07, -3.1624037713552806e-09, 1.0794827959578473e-11, + -7.7950231795753708e-15, -1.9036326199741433e-16, 1.1391135385648778e-18, 1.210707461009259e-21, + -5.3933004841533438e-23, 1.9910325335995957e-25, 2.0838450578394491e-27, -2.2565767244958724e-29, + -3.9130112051077074e-32, 1.9331575662520048e-33, -5.8666884923938355e-36, -1.4910781892036265e-37, + 1.2823302589475358e-39, + }}, + {19.860834841883541, 19.859028062896897, 83.704452082133344, 379.60364653402274, 139166.31452306718, { + 0.0023756249488557225, 8.5352010542104159e-09, 9.1996411812662828e-14, 1.6526343342266725e-18, + 4.1563361787300572e-23, 1.3440037497962314e-27, 5.2525821584148472e-32, 1.2183799218081628e-35, + -1.5735077432780944e-37, 2.5176847315945519e-39, -3.961628589371911e-41, 6.1306404718678053e-43, + -9.3218111125163812e-45, 1.3912661796033864e-46, -2.0357612563184841e-48, 2.916521338217192e-50, + -4.084556141904497e-52, + }}, + {641872.87532913941, 3809.9938565182738, 0.053847671162906827, 552.1421266475005, 313550.5899054535, { + 1.7285749561626197e-31, 2.4318120923830719e-33, 3.3941862367366062e-35, 4.6991766381343078e-37, + 6.4520889669397511e-39, 8.7837195864359926e-41, 1.1853804652318628e-42, 1.5853822185658181e-44, + 2.1008562777906404e-46, 2.757572400277024e-48, 3.5842633115390591e-50, 4.6119089090375241e-52, + 5.8725274467381952e-54, 7.3973763871534414e-56, 9.2144898350157561e-58, 1.1345525058709156e-59, + 1.3801961307647577e-61, + }}, + {31.749625445361385, 31.566507898605266, 68.544337476203921, 390.4353626316481, 150035.43391273171, { + 0.0022879437596230118, 7.7426034878552362e-09, -8.6372248081337465e-13, 7.1003792465672502e-15, + -5.0296560667553335e-17, 3.2961243692586238e-19, -1.9538945573717341e-21, 1.0060767668252843e-23, + -4.100708878178949e-26, 9.1777458790333055e-29, 3.608453955220017e-31, -5.6033412041568162e-33, + 3.2425210042158384e-35, -5.2187328607769378e-38, -8.7857522912019486e-40, 9.2010455626660619e-42, + -3.1900155760637459e-44, + }}, + {1250.7705378495664, 1184.7682830459275, 89.40427171334585, 3239.7643381998382, 10554518.021943556, { + 4.9317339655024406e-41, 1.3757047379286675e-43, 3.814362470667139e-46, 1.0510942293354236e-48, + 2.878287340896251e-51, 7.8315400512648646e-54, 2.1170260122608838e-56, 5.6847878268694997e-59, + 1.5161837335672633e-61, 4.0158418273045379e-64, 1.0561443013844546e-66, 2.7575494280781714e-69, + 7.1466869272775034e-72, 1.8381962978708147e-74, 4.6914263391249435e-77, 1.1878445042576163e-79, + 2.9830940127366719e-82, + }}, + {0.028327008616249882, 0.022516292196339053, 115.4234492469364, 6.6813977894844907, 474.78025583797864, { + 5.6283116035590594e-103, 3.9169159368552933e-103, 2.7241195728476369e-103, 1.8933113629143001e-103, + 1.3150129024394183e-103, 9.1274234799294607e-104, 6.3310347131245936e-104, 4.3884137382116522e-104, + 3.039798947634391e-104, 2.1041882904684421e-104, 1.4555421300829754e-104, 1.0061514906052873e-104, + 6.9502164348233492e-105, 4.7976389607836666e-105, 3.309397763219521e-105, 2.2811845111973387e-105, + 1.5713022642630638e-105, + }}, + {0.028762826118416858, 0.028762806057089092, 91.243918117197509, 1.2450378762151741, 239.46371232019948, { + 1.3114105832093972e-91, 1.2116002004194737e-91, 1.1191902808569769e-91, 1.0336462454625914e-91, + 9.5447113929996517e-92, 8.8120305130787054e-92, 8.1341270649459786e-92, 7.5070121934233062e-92, + 6.9269799790242382e-92, 6.3905878873152212e-92, 5.894638534460008e-92, 5.436162682607779e-92, + 5.0124033843049901e-92, 4.6208012002866081e-92, 4.2589804198549626e-92, 3.9247362176036426e-92, + 3.6160226845081808e-92, + }}, + {504.57899132066751, 4.2387826131661166, 0.14762204423320782, 31.863922025319713, 1308.9396225315384, { + 1.2763375609800877e-11, 1.566178307382788e-12, 1.8802203550962336e-13, 2.20464434641588e-14, + 2.5198411971164315e-15, 2.8008509687741818e-16, 3.018890297422439e-17, 3.1441586797706448e-18, + 3.1499059565124538e-19, 3.0174283754691307e-20, 2.7412729539728768e-21, 2.3335495214267902e-22, + 1.8260027074452397e-23, 1.2685355615482992e-24, 7.2334065991652599e-26, 2.5474170708712085e-27, + -8.3816348310896705e-29, + }}, + {0.28241799623115577, 0.28241299403255005, 1.8011890034094586, 1.6207049344688518, 0.91625980133843921, { + 0.60203843540399038, 0.28628531104675348, 0.12956841877122094, 0.060224221732723564, + 0.035170262385194537, 0.02978328407529534, 0.030983568451079999, 0.032609765196263664, + 0.032406138674955773, 0.030125365758464355, 0.026387133936551729, 0.022035652092325265, + 0.01781561249866374, 0.014242305575257984, 0.011582341021726329, 0.0098905620088674158, + 0.009069197915218364, + }}, + {1821.9924515980522, 2.6062229116299851, 0.034300985020018021, 32.891835188063325, 1498.6312072959706, { + 2.1219057583336731e-18, 3.2431272981377124e-19, 4.8978171515604788e-20, 7.3049712572356634e-21, + 1.0753932815305137e-21, 1.5616366786723635e-22, 2.2354198375127705e-23, 3.1519098044884488e-24, + 4.3737474121956612e-25, 5.9673465647119496e-26, 7.9961642039109489e-27, 1.0510179909839586e-27, + 1.3531141041503353e-28, 1.7033806095735994e-29, 2.0924546148436522e-30, 2.5020261640049042e-31, + 2.9032518715911119e-32, + }}, + {3752.5975554047, 681.75768190662609, 2.6994789306718561, 375.14647095297249, 150519.86730517593, { + 1.9745661407452403e-76, 6.5463530873888349e-78, 2.1640148769461917e-79, 7.1325139903035132e-81, + 2.3438772430983256e-82, 7.6793588051094788e-84, 2.5084225938127382e-85, 8.1686374205801841e-87, + 2.6519111408440306e-88, 8.5825380454086115e-90, 2.7688959424901326e-91, 8.9046750585883303e-93, + 2.8545402774873817e-94, 9.1210898518651462e-96, 2.9049270520256148e-97, 9.2211865840098137e-99, + 2.9173347602617714e-100, + }}, + {0.016083171428850837, 0.016082445080071141, 112.0757523690494, 0.212641009127504, 202.02942654352773, { + 9.1449597973869443e-89, 9.0529735984044307e-89, 8.9617878846924829e-89, 8.871396597230328e-89, + 8.7817937172889513e-89, 8.6929732661932423e-89, 8.6049293050854056e-89, 8.5176559346896158e-89, + 8.4311472950779368e-89, 8.3453975654374605e-89, 8.2604009638386899e-89, 8.1761517470051584e-89, + 8.092644210084259e-89, 8.0098726864192931e-89, 7.9278315473227406e-89, 7.8465152018507215e-89, + 7.7659180965786645e-89, + }}, + {1.7475458847450525, 1.5844951504853213, 9.7707970067602954, 2.3873685971904619, 184.00356714716651, { + 9.0568650022475563e-58, 7.5171231866427131e-58, 6.2351292700808339e-58, 5.1683993660181793e-58, + 4.2813447030743203e-58, 3.5441696927274333e-58, 2.9319434069358988e-58, 2.4238175930811729e-58, + 2.0023684544879451e-58, 1.6530429093501941e-58, 1.3636930003879763e-58, 1.1241846398964556e-58, + 9.2606900655605882e-59, 7.6230671823316925e-59, 6.270364374556645e-59, 5.1538086460341074e-59, + 4.2328417331721952e-59, + }}, + {551.10009432393201, 0.011264865394650847, 0.003292572361745856, 9.7814919667294813, 292.28438939164681, { + 1.1459177817233248e-26, 4.9675239931579856e-27, 2.1425109673753214e-27, 9.1927009138493678e-28, + 3.9231720795687163e-28, 1.6650961984921905e-28, 7.027142523686178e-29, 2.9483642658267097e-29, + 1.2296097780709981e-29, 5.0962770182222616e-30, 2.09869073360089e-30, 8.5853402847605843e-31, + 3.488006670079647e-31, 1.4070101884533327e-31, 5.6337578138445256e-32, 2.2384702887935886e-32, + 8.8230166916960915e-33, + }}, + {15122.364049670126, 2.6688143337451766, 0.0026213787934185164, 19.170634830674874, 588.81732472827844, { + 1.0575301043278944e-14, 2.2706870578574565e-15, 4.8066434793404775e-16, 1.0023175341309439e-16, + 2.0571448235574361e-17, 4.1513398538459537e-18, 8.2277910992746411e-19, 1.5994975765147338e-19, + 3.045305930945348e-20, 5.6682607983811085e-21, 1.0292310154572557e-21, 1.818404688149759e-22, + 3.1158493522725949e-23, 5.1566880278109538e-24, 8.1976551550508837e-25, 1.242308387056327e-25, + 1.7747122842751251e-26, + }}, + {2708.029735278365, 358.1023357140445, 0.85742789905012606, 108.13383380280447, 15055.49104419138, { + 9.7429640352274025e-97, 1.1597166173289686e-97, 1.3775815957966503e-98, 1.6329810935763009e-99, + 1.9316878956713902e-100, 2.2802278395347198e-101, 2.6859530508649034e-102, 3.1571177490003869e-103, + 3.7029557852688122e-104, 4.3337588632275413e-105, 5.0609541543838335e-106, 5.8971797685133861e-107, + 6.8563562668860847e-108, 7.9537521234496457e-109, 9.2060407484351236e-110, 1.0631346397383128e-110, + 1.224927600410175e-111, + }}, + {39.838079010808684, 15.033208755925541, 2.0431659696564344, 14.343156169719318, 440.70946919524624, { + 1.1043634834663842e-22, 3.5494998131135849e-23, 1.1324782605939498e-23, 3.585878691136397e-24, + 1.1265552821268078e-24, 3.5105955934054195e-25, 1.0848051633436946e-25, 3.3229714737830523e-26, + 1.0086848026979429e-26, 3.0330231960186313e-27, 9.0304705295772037e-28, 2.6611293637889093e-28, + 7.757650697096774e-29, 2.2359848944462108e-29, 6.3682803208709085e-30, 1.7910196070317601e-30, + 4.9702456821792988e-31, + }}, + {0.069739456749215722, 0.0096700563195362271, 4.992856753035058, 6.7289572329531948, 12.537944090445187, { + 0.25028228991407103, 0.0099816194638970759, 0.0011933343726752766, 0.00023836626157271493, + 6.598691487451473e-05, 2.379900846272392e-05, 1.0162229397587723e-05, 5.1770369754023742e-06, + 2.94441028528062e-06, 1.8487135017324422e-06, 1.2617201938724397e-06, 9.1277157687285336e-07, + 6.934244613051576e-07, 5.4870095350569402e-07, 4.4779220800031746e-07, 3.7462348269170212e-07, + 3.2005596106290044e-07, + }}, + {0.43347740933360429, 0.4318786868543249, 22.543445664795509, 19.011488767133915, 221.11171831192931, { + 0.059599011939448034, 0.00013477130509151447, 9.1426400352120949e-07, 1.0340126294628892e-08, + 1.6290352830183913e-10, 3.5048562036990502e-12, 4.2979656897274825e-14, 1.0980342811671824e-14, + -1.6197038304933102e-15, 3.1278070702684773e-16, -4.9776427305123125e-17, 6.7281795494917482e-18, + -6.3922732783801102e-19, 6.7143486213922256e-21, 1.3553714803937278e-20, -3.3165664628802219e-21, + 3.8558729991867105e-22, + }}, + {45.931005276389364, 45.92376271013724, 5.8489077468931354, 35.117520186657821, 1571.5347192057443, { + 1.7434698821270893e-12, 2.0421281973749637e-13, 2.3441801216768422e-14, 2.633444594566941e-15, + 2.8904836135973719e-16, 3.0938119843545318e-17, 3.2217976854631669e-18, 3.2551838696403374e-19, + 3.1800134892088493e-20, 2.9905723580826127e-21, 2.6918210551489772e-22, 2.3007070639018622e-23, + 1.8457859830170973e-24, 1.3647739139099065e-25, 9.0001277059658611e-27, 4.9232012901086964e-28, + 1.7422207053362458e-29, + }}, + {7.5245236150149299, 7.5244122113634226, 43.568190516011299, 132.84574084586865, 14283.166067143111, { + 0.0074153678183698753, 2.595841770483973e-07, 2.7261201315044163e-11, 4.7715613588214918e-15, + 1.1692410966426307e-18, 3.6837665473872416e-22, 1.4185031466684007e-25, 6.4553407906912502e-29, + 3.3896585464729707e-32, 2.0172066549936284e-35, 1.3416817484551243e-38, 9.863120188167488e-42, + 7.9412282704497966e-45, 6.9498144118741116e-48, 6.5687463213166074e-51, 6.668467005938962e-54, + 7.2365775281311996e-57, + }}, + {0.018737146827942343, 0.018737142925348879, 23.281607535443719, 0.25747596801714617, 10.156158695615463, { + 3.095817655593249e-06, 3.030929668216741e-06, 2.9671448843245319e-06, 2.9044480498971928e-06, + 2.8428240872325206e-06, 2.7822580932000916e-06, 2.7227353375109006e-06, 2.664241261001976e-06, + 2.6067614739358487e-06, 2.5502817543147672e-06, 2.4947880462095406e-06, 2.4402664581029016e-06, + 2.386703261247272e-06, 2.3340848880368261e-06, 2.2823979303937368e-06, 2.2316291381684985e-06, + 2.181765417554216e-06, + }}, + {11049.078861563739, 93.227062805786119, 0.61809860872929301, 692.49714922935902, 500294.44078660791, { + 1.0464227787090027e-100, 2.197263929820586e-102, 4.6035691331781071e-104, 9.6236088167717713e-106, + 2.0072684121235159e-107, 4.177249183611295e-109, 8.6733306198613606e-111, 1.796739199162247e-112, + 3.7134769571883173e-114, 7.6571259081095238e-116, 1.5751898628771766e-117, 3.2327690223321233e-119, + 6.6188654631070438e-121, 1.3519227614244818e-122, 2.7546845103374728e-124, 5.5993239699166157e-126, + 1.1353637263347735e-127, + }}, + {85.9414338520602, 85.71746978788714, 0.27057062698758932, 12.867899322062678, 6.3080792875529061, { + 0.35272013457268081, 0.027813427040395414, 0.0064693822704733539, 0.00241954654078961, + 0.0011980903047173511, 0.0007103021236767133, 0.00047493021040101083, 0.00034499916875532061, + 0.00026580666352993726, 0.00021378811545811493, 0.00017758540893808222, 0.00015121600306671802, + 0.00013129508057125038, 0.00011579186260764652, 0.00010342702723036966, 9.3360849005856436e-05, + 8.5022536684168745e-05, + }}, + {16.264706933594681, 16.262759606631217, 6.2011175915236771, 29.187646570198222, 625.51564454860272, { + 0.035434462700264263, 2.8324212356643119e-05, 6.7920309906882957e-08, 2.7174172150847691e-10, + 1.4762538882470418e-12, 1.7226018850931785e-14, -7.9051728821683975e-16, 1.1927328279855993e-16, + -1.4805096712543525e-17, 1.7202630516612313e-18, -1.8117819513273843e-19, 1.6669241460348731e-20, + -1.2291008609393661e-21, 5.2577474442237143e-23, 2.8943392250541413e-24, -9.3813486309825711e-25, + 1.1011109627253892e-25, + }}, + {4982.2835515011284, 121.61089759262231, 0.52293517638147868, 233.20677996859757, 55818.749152719851, { + 2.9558260859910372e-08, 4.0096193772151936e-10, 5.1889025695931251e-12, 6.3605398376394355e-14, + 7.3128134625805859e-16, 7.7729239193020889e-18, 7.4639035064023993e-20, 6.2057716615286737e-22, + 4.0439743536441842e-24, 1.355501376733891e-26, -1.14747815873501e-28, -2.598344680660497e-30, + -2.3627992611907396e-32, -5.1591700948644842e-35, 1.8800399595171628e-36, 3.0702807946897477e-38, + 1.6580704713899132e-40, + }}, + {0.056567899442936193, 0.056314554141663635, 65.388646539543714, 19.206708146732904, 242.95405424178279, { + 0.056856850467489797, 0.00011701362795152814, 7.2195538470433857e-07, 7.5368293966123068e-09, + 8.5958249457622032e-11, 6.0461901619902558e-12, -6.7788550127481889e-13, 1.1797019734166367e-13, + -1.654566822457917e-14, 1.9376273689942378e-15, -1.5207246855059605e-16, -1.344451403036363e-18, + 3.0590290634332327e-18, -5.7122249647897576e-19, 4.2798559322600186e-20, 5.4772531490941244e-21, + -2.1191086318251584e-21, + }}, + {20464.544354342899, 0.28788971470417862, 8.483773080076314e-05, 2.2077060704325331, 10.470226987724185, { + 0.019988453087608882, 0.0092510921603180044, 0.0039603922779805224, 0.0015115511584015895, + 0.00047363080004910536, 9.0965305360069852e-05, -1.5698378949678407e-05, -2.3865388138452947e-05, + -8.2971443232641283e-06, 4.9778308612308817e-06, 1.0638613546689644e-05, 1.0515576531678736e-05, + 7.7142748821974088e-06, 4.5520686011493928e-06, 2.2001168037709656e-06, 9.5112540544150981e-07, + 6.096328896511779e-07, + }}, + {14264.894367095838, 7393.1223255459054, 0.57051672359541439, 109.4766319905295, 8958.7204893420912, { + 0.0093631492733329685, 5.2257179384444541e-07, 8.7496612010520811e-11, 2.4416603943222914e-14, + 9.5390981226556894e-18, 4.7915259330858964e-21, 2.9416469308674419e-24, 2.1343120452731696e-27, + 1.7867886779805449e-30, 1.6952983164172791e-33, 1.7977270331320374e-36, 2.1070122536296044e-39, + 2.7046988401487562e-42, 3.7738352861974704e-45, 5.6868362423267503e-48, 9.2043418043722778e-51, + 1.592496363040873e-53, + }}, + {6.2832816253877208, 6.2746637183390028, 23.19544022514869, 49.428421143426469, 3385.2273038420099, { + 2.5403305528202989e-37, 3.8506055121763648e-38, 5.8051138720282315e-39, 8.7034188768961305e-40, + 1.2975349122294484e-40, 1.9233078228632091e-41, 2.8341804562862254e-42, 4.1514632355117252e-43, + 6.0438465550928355e-44, 8.7438973764317199e-45, 1.2569384863948673e-45, 1.7950443824293602e-46, + 2.5463613214645787e-47, 3.5873716195048863e-48, 5.0184324329443701e-49, 6.969710446477893e-50, + 9.6079603582847104e-51, + }}, + {0.022502298038901723, 0.022502294945193031, 37.689471369312514, 9.1925393762735297, 31.964434422459412, { + 0.15675135617953068, 0.0024519935890897565, 0.00011504726213644223, 9.0085000669320812e-06, + 9.8074391057289521e-07, 1.4087670718021282e-07, 2.2990852520928505e-08, 5.1589720858072816e-09, + 1.0571904706824842e-09, 3.1483251139518879e-10, 9.2409196708229335e-11, 2.7360550576716899e-11, + 1.1217897142405691e-11, 4.2223546415711103e-12, 1.6564075091734988e-12, 8.2017381879223257e-13, + 3.9459832751624237e-13, + }}, + {227150.26233165391, 378.96032601152694, 0.10125496208353242, 1181.0466395148092, 1395935.5607489245, { + 0.00019653443889627715, 1.4626214393965e-07, 5.5903145793583194e-11, -3.0986162202577227e-14, + -5.1907052994476143e-17, 1.3395571385873596e-20, 7.952541603579054e-23, 6.5937823289476636e-27, + -1.6839821275103946e-28, -8.1493153116694417e-32, 4.5114736082692035e-34, 4.3726994469963788e-37, + -1.4468458312402332e-39, -2.2844170318209972e-42, 5.3325944168358248e-45, 1.271342958004725e-47, + -2.1793318839440247e-50, + }}, + {132964.6171080668, 2.4277187983974349, 0.00022326941760053237, 3.6865349788163781, 363.02143772002569, { + 2.3969866239861916e-106, 1.9405781710122116e-106, 1.57045694042297e-106, 1.2704255629108061e-106, + 1.0273055659782445e-106, 8.3037869975346123e-107, 6.7093091693203154e-107, 5.4188027835471606e-107, + 4.3747344837992165e-107, 3.5303822663862573e-107, 2.8478184394885467e-107, 2.2962661977107594e-107, + 1.8507611033745881e-107, 1.4910613110453502e-107, 1.2007606427941007e-107, 9.6656702784824818e-108, + 7.7771569774504321e-108, + }}, + {0.10941013600105379, 0.10824596704721813, 7.8438471792570619, 0.56073874851206451, 6.8039582273278318, { + 0.00056701837879366893, 0.00051655435153916758, 0.00046984009489084151, 0.00042664358118327088, + 0.00038674497231674416, 0.0003499360656388164, 0.00031601976197690585, 0.00028480955503275447, + 0.00025612904137491691, 0.00022981145029003704, 0.00020569919277755149, 0.00018364342899563252, + 0.00016350365348867974, 0.00014514729754849971, 0.00012844934808249763, 0.00011329198238276278, + 9.9564218209874768e-05, + }}, + {29251.213605338689, 125.22745684140034, 0.087560343688860134, 218.1786068094882, 52384.53977871473, { + 2.007159853989601e-54, 9.4243573745361942e-56, 4.4069074749213887e-57, 2.0521318503363983e-58, + 9.5157232086260172e-60, 4.393583840856443e-61, 2.0198133738920957e-62, 9.2446960381892122e-64, + 4.2124634783543349e-65, 1.9107945092132879e-66, 8.6277481361329698e-68, 3.8775459137133649e-69, + 1.7344502618744167e-70, 7.7211238813552843e-72, 3.4204262080958557e-73, 1.5077421571065584e-74, + 6.6128089240394495e-76, + }}, + {19692.773719360728, 794.75402519613738, 0.30358394013760454, 214.0491179323661, 44971.630287519234, { + 0.0041685418957512462, 1.5488197378042302e-07, -1.0102321783840648e-09, 8.2253325741242063e-12, + -5.4087578122731397e-14, 2.2891953396576516e-16, 2.7385460229807945e-19, -1.52059018982002e-20, + 1.222448348174968e-22, 5.123484409763934e-26, -1.1295323933252811e-26, 9.8640804278955845e-29, + 3.4573467963429578e-31, -1.5171936917507812e-32, 9.2672646245318472e-35, 1.337750187399676e-36, + -2.6532493596757421e-38, + }}, + {66.488421131237857, 53.66076332629526, 7.7353519434149689, 72.571167073037913, 4929.414633703901, { + 0.012617257643637558, 1.4718694365154773e-06, -6.0122637430204105e-09, 1.9408799722839762e-10, + -5.1440588748283212e-12, 1.1155948154858138e-13, -1.6078095769995741e-15, -3.1474653469675895e-18, + 1.0655169596315789e-18, -3.2309620453741294e-20, 1.9088294868511668e-22, 2.2887564584064119e-23, + -9.2117035248938169e-25, 4.1579387104958948e-27, 9.7104768646639104e-28, -3.595142127913283e-29, + -2.2978985666554542e-31, + }}, + {2301.849830239702, 56.77327427719753, 0.055783105214864522, 13.172727436374407, 290.41256645085883, { + 1.1618660146113516e-09, 2.7407010871702709e-10, 6.3180984819921244e-11, 1.4203819126475703e-11, + 3.1059965296552597e-12, 6.5856405234789278e-13, 1.348525787458144e-13, 2.6529124869867727e-14, + 4.9788012857963252e-15, 8.8244664964215532e-16, 1.454417745402778e-16, 2.170809719938638e-17, + 2.7803270356888798e-18, 2.6249588217635451e-19, 4.669672084768746e-21, -5.3277388746019483e-21, + -1.5520176745527182e-21, + }}, + {4.5602639107327949, 0.12457186365328925, 2.1847587441514875, 13.263680005085499, 796.83236269718986, { + 3.2842519561855695e-101, 1.74700146391831e-101, 9.2832867976997488e-102, 4.9278671735941994e-102, + 2.6131355389149198e-102, 1.3842277415899821e-102, 7.3247489092804143e-103, 3.8718131089739559e-103, + 2.0444161823940294e-103, 1.0783364857112619e-103, 5.6815403897344394e-104, 2.9902080311227877e-104, + 1.57201566066359e-104, 8.2552285005799589e-105, 4.3302667921513166e-105, 2.2688722642749717e-105, + 1.1874401661085348e-105, + }}, + {340.19249498105739, 0.19189639382894708, 0.017797980039490475, 1.4143900370321849, 191.03989645392898, { + 2.0240842005497599e-70, 1.828118931425708e-70, 1.6506442380840566e-70, 1.4899596660420454e-70, + 1.3445171846217697e-70, 1.21290786858034e-70, 1.093849714239213e-70, 9.8617649588831543e-71, + 8.8882757587165674e-71, 8.0083858878946071e-71, 7.2133292672387132e-71, 6.4951395835394865e-71, + 5.8465792031068372e-71, 5.2610742417084756e-71, 4.7326552713358189e-71, 4.2559031869711363e-71, + 3.8258997958381477e-71, + }}, + {132652.84489771849, 2326.2386357905875, 0.25976026055211737, 716.44622526610033, 510416.01759186009, { + 0.0012377110802578647, 9.7558364719290608e-09, -2.4033382437251204e-11, 5.9259358660425925e-14, + -1.1957634196772284e-16, 1.6209230625389578e-19, 1.2685640752991232e-23, -8.2852892653279426e-25, + 2.2502881653681156e-27, -6.3693200010331147e-31, -1.5825721262204702e-32, 4.9941869880380603e-35, + 1.4980918731859606e-38, -5.793280562822666e-40, 1.444499308194608e-42, 3.3292295705431789e-45, + -2.9084562103272123e-47, + }}, + {280107.12187589868, 114168.80819620333, 3.834254509826966, 3180.1502224263095, 10103287.439935764, { + 0.00027530807365976838, 2.0306855045983119e-09, -1.0044559386003003e-12, 4.0037900176394997e-16, + -9.9942411137195887e-20, -9.6784190365728727e-24, 2.45974621813941e-26, -9.8471290330565173e-30, + -2.4022190111018314e-33, 4.6056016145137725e-36, -1.3401918316584992e-39, -1.3842667576645505e-42, + 1.3514709811073143e-45, 8.1402037690483599e-50, -8.4286439762900615e-52, 3.6660479298962207e-55, + 4.0162864117136799e-58, + }}, + {2603.9259321548438, 0.19078497951246232, 0.0021122314123842192, 0.85253121742254101, 158.5608873915601, { + 2.3620978110944657e-63, 2.2175561750939901e-63, 2.0814534877813299e-63, 1.9533195062602209e-63, + 1.8327090070442291e-63, 1.7192005143840245e-63, 1.6123950903785564e-63, 1.5119151839989157e-63, + 1.417403536281387e-63, 1.3285221390684732e-63, 1.2449512447937311e-63, 1.1663884249183122e-63, + 1.0925476747343509e-63, 1.0231585623529894e-63, 9.5796541979306125e-64, 8.967265741804494e-64, + 8.392136171580861e-64, + }}, + {33.252648324796453, 33.142347339916888, 40.227829060957433, 231.23057203725838, 53991.128646198544, { + 0.00021022389566405865, 1.2952928930346969e-06, 6.3220414327655756e-09, 1.8925694638473113e-11, + -2.4404109079576897e-14, -6.4660327594132809e-16, -2.2827577609394652e-18, 1.8714825908304148e-20, + 2.1995857315985864e-22, -1.2178285992890911e-25, -1.6941301923492615e-26, -7.4579479117749622e-29, + 1.1984355107034525e-30, 1.3644332311605915e-32, -6.4663342720947051e-35, -1.9736494256166414e-36, + -1.6082682221499003e-39, + }}, + {132080.03446673267, 10.498605026636712, 0.00047371629517240144, 9.2655132536658655, 372.88885063263513, { + 1.9430637040942526e-47, 1.0183054632415227e-47, 5.324415200232397e-48, 2.7775000843399191e-48, + 1.4454719784499647e-48, 7.5045278723171201e-49, 3.8866854487936694e-49, 2.0079897068979355e-49, + 1.0347915905508301e-49, 5.3190779374656327e-50, 2.7270560429095535e-50, 1.3944650712312719e-50, + 7.1114491747855153e-51, 3.6168114596604578e-51, 1.8343856997245973e-51, 9.2775117100698335e-52, + 4.6787166705444951e-52, + }}, + {7308.7194358195011, 118.18342902215862, 0.038510721840007152, 39.927431805282161, 670.33101342880832, { + 0.034229478635430098, 2.5531773071592635e-05, 5.7132459695535648e-08, 2.1307554980671995e-10, + 1.1125315842226399e-12, 7.4685372281877532e-15, 6.1278613001836254e-17, 5.9420044206896562e-19, + 6.6482129369515432e-21, 8.4301350872958934e-23, 1.1947274066831828e-24, 1.8714094259202444e-26, + 3.2105344922055356e-28, 5.9868453567874982e-30, 1.2057089810482844e-31, 2.6080816603984951e-33, + 6.0306423134739816e-35, + }}, + {1595.3135686482347, 12.150468253796888, 0.45411425377557646, 214.18272170301557, 43194.591440743694, { + 0.0042641282219074107, 4.9359515620805642e-08, 1.714086670613734e-12, 9.920725105422446e-17, + 8.038630437584396e-21, 8.3747658797824146e-25, 1.0621989539055848e-28, 2.7855948959535712e-32, + -3.2813436473116358e-34, 9.1205001283523314e-36, -2.4697343367163952e-37, 6.5654942027004723e-39, + -1.7114174726653441e-40, 4.3687623027220494e-42, -1.0905110561178654e-43, 2.6571517377841299e-45, + -6.3069765911036073e-47, + }}, + {146.7729527309848, 119.84194205799446, 1.1077423150579067, 21.982003462596197, 220.57724901989246, { + 0.059671174030763406, 0.00013526139775499251, 9.1982331602200364e-07, 1.0425183468706814e-08, + 1.6542114974506375e-10, 3.3747595328186163e-12, 8.4148195305618387e-14, 2.4796903212687867e-15, + 8.4313670118128076e-17, 3.2490485724454353e-18, 1.3993266104922913e-19, 6.6611264494891714e-21, + 3.4728400978428508e-22, 1.9680412169555761e-23, 1.2044993658252161e-24, 7.9180400364218429e-26, + 5.5637624782757993e-27, + }}, + {0.026598917980722787, 0.026598900065803836, 65.411651232664568, 0.69237861758730901, 113.80844454540282, { + 1.4121940383802969e-46, 1.3333141443830061e-46, 1.2585511611883553e-46, 1.1877041273145916e-46, + 1.120581267238574e-46, 1.0569995976253541e-46, 9.9678454941333347e-47, 9.3976960515422476e-47, + 8.8579595102851062e-47, 8.3471214297782332e-47, 7.8637378641570602e-47, 7.4064322899757106e-47, + 6.973892659493729e-47, 6.5648685747257402e-47, 6.1781685776043243e-47, 5.8126575517749191e-47, + 5.4672542317043422e-47, + }}, + {17.674431603742239, 0.28466465102867705, 0.12402946062071656, 6.5619615595333025, 16.881346396534575, { + 0.21563940943194282, 0.0064229468193943359, 0.0005491967264379886, 9.2472924937904674e-05, + 1.4402014746789553e-05, 5.3528764087669943e-06, 1.5081852528518887e-06, 5.0467396058601932e-07, + 2.8025497350454272e-07, 1.3287317460864603e-07, 6.5393901209134507e-08, 4.2544805509125341e-08, + 2.8638928770568563e-08, 1.8898955905948081e-08, 1.359330982719048e-08, 1.0557859026305654e-08, + 8.3380958924725612e-09, + }}, + {15736.991120017041, 15176.750911137795, 17.105125116108749, 2189.3543049271202, 4774381.3346473034, { + 0.00040558920611620379, 4.2475993665341232e-11, 1.2521905809946976e-17, 1.5895781368392299e-21, + -2.9548174368750865e-24, 5.3412042280767766e-27, -9.315368294964647e-30, 1.5607664325690505e-32, + -2.4981574330490277e-35, 3.7910728572773312e-38, -5.3962460321608816e-41, 7.0867475367330937e-44, + -8.3480069864188012e-47, 8.3271314272928868e-50, -5.959687685716697e-53, 4.3965743612241916e-57, + 7.8614223041616298e-59, + }}, + {1.6989549191884179, 1.2250431289088841, 52.062863762052174, 77.535958819671976, 6386.5811980676244, { + 4.2286076548588953e-06, 1.3597045381270068e-07, 4.0716642617971884e-09, 1.1160569892206461e-10, + 2.7215154616626171e-12, 5.5864680103556312e-14, 8.346246035325784e-16, 3.3033231236430025e-18, + -2.9630162527834376e-19, -1.0947535361108319e-20, -1.49728732322839e-22, 3.1847689385376156e-24, + 2.1799970945678627e-25, 4.0253502340088612e-27, -8.1603823175994564e-29, -6.7244049561941123e-30, + -1.2164189383131703e-31, + }}, + {482.57891334402569, 0.81223981464729877, 0.0012595212772432921, 1.7694262548891022, 0.45484500024342361, { + 0.78755824771610594, 0.33531048623220305, 0.14428359499155444, 0.076839427074070321, + 0.060133525367842812, 0.058899027552088493, 0.058482105924879201, 0.054642290601476462, + 0.047709404514832715, 0.039461867410130913, 0.031638448820668158, 0.025377569941324788, + 0.021140459777503732, 0.018850853111822569, 0.018099657791352263, 0.018336590648848714, + 0.019015589646321718, + }}, + {65.685222889514648, 65.537799877303982, 0.68950875884693019, 8.3296166673661602, 31.298467466489974, { + 0.15840160548006818, 0.0025350042468166283, 0.00011929016185201481, 1.0497676089671391e-05, + 7.9966282161315515e-07, 2.3590200533071778e-07, 1.2067303045929503e-08, 5.9808085825804201e-09, + 2.4122584050725408e-09, 2.3323981844829462e-11, 1.0698535229379641e-10, 7.7396585048797632e-11, + 4.8851347905543226e-12, 2.1837946990112896e-12, 4.3538020843259685e-12, 1.2750921034597681e-12, + 1.879245190038893e-13, + }}, + {21212.422036015429, 3.4224319323779411, 0.010838866886208363, 131.88021413911818, 15445.908730198598, { + 0.0071307992239951282, 2.3083132719972516e-07, 2.2416744579270592e-11, 3.6282657386553078e-15, + 8.2215493481880063e-19, 2.3952603056960511e-22, 8.5290751818366974e-26, 3.5892345180399544e-29, + 1.7428080071483545e-32, 9.5910057436080357e-36, 5.8875670816777289e-39, 4.635926408850851e-42, + -3.1409500122580381e-44, 1.868667669114046e-45, -9.9940575108012585e-47, 5.2796236703634793e-48, + -2.7495825354554203e-49, + }}, + {15.203962940253065, 15.203587867889706, 76.833230739364126, 310.50501673452027, 89756.458098108444, { + 0.0029580947962933837, 1.6478451016081948e-08, 2.7538605074083055e-13, 7.6703687003730515e-18, + 2.9910149108113535e-22, 1.499565310825778e-26, 9.188875524173109e-31, 6.6544170940702399e-35, + 5.5603941223900167e-39, 5.2657325212915444e-43, 5.5733548328734576e-47, 6.519890265857603e-51, + 8.3535758480361161e-55, 1.1633669633700935e-58, 1.7497854013278153e-62, 2.8267479418075956e-66, + 4.8814975566589909e-70, + }}, +}; +// clang-format on diff --git a/include/libint2/boys.h b/include/libint2/boys.h index 9587502a9..69da55012 100644 --- a/include/libint2/boys.h +++ b/include/libint2/boys.h @@ -58,6 +58,7 @@ #include #include #include +#include #if HAVE_LAPACK // use F77-type interface for now, switch to LAPACKE later extern "C" void dgesv_(const int* n, const int* nrhs, double* A, const int* lda, @@ -349,7 +350,7 @@ class FmEval_Chebyshev7 { inline void eval(Real* Fm, Real x, int m_max) const { // large T => use upward recursion // cost = 1 div + 1 sqrt + (1 + 2*(m-1)) muls - if (x > T_crit) { + if (x >= T_crit) { const double one_over_x = 1 / x; Fm[0] = 0.88622692545275801365 * sqrt(one_over_x); // see Eq. (9.8.9) in Helgaker-Jorgensen-Olsen @@ -2109,6 +2110,13 @@ struct GenericGmEval : private GmEvalFunction { }; namespace detail { +/// the truncated Coulomb evaluator keeps its mutable work in Engine's scratch +template +struct CoreEvalScratch> + : os_core_ints::truncated_coulomb_gm_eval::scratch_type { + explicit CoreEvalScratch(int mmax) + : os_core_ints::truncated_coulomb_gm_eval::scratch_type(mmax) {} +}; /// Boys function scratch for core evaluators that need Boys function values /// in a separate buffer; lives on the caller's stack, so that evaluators /// can be invoked through a const reference from any number of threads diff --git a/include/libint2/boys_fwd.h b/include/libint2/boys_fwd.h index d03e440d1..a21bcb6c5 100644 --- a/include/libint2/boys_fwd.h +++ b/include/libint2/boys_fwd.h @@ -62,6 +62,8 @@ template struct erfx_coulomb_gm_eval; template struct q_gau_gm_eval; +template +class truncated_coulomb_gm_eval; } // namespace os_core_ints /* diff --git a/include/libint2/engine.h b/include/libint2/engine.h index badd60654..f2a1dabcc 100644 --- a/include/libint2/engine.h +++ b/include/libint2/engine.h @@ -273,6 +273,9 @@ enum class Operator { /// these components are not σ·σ-folded; the user multiplies them by -1 /// (see \ref rkb_operators). op_coulomb_op, + /// Spherically truncated Coulomb, theta(Rc-r12)/r12; Rc has length units. + /// Interpolation-free evaluation with binary64 order recurrences. + truncated_coulomb, // do not modify this invalid = -1, // !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!keep this @@ -280,7 +283,7 @@ enum class Operator { first_1body_oper = overlap, last_1body_oper = opemuop, first_2body_oper = delta, - last_2body_oper = op_coulomb_op, + last_2body_oper = truncated_coulomb, first_oper = first_1body_oper, last_oper = last_2body_oper }; @@ -731,6 +734,17 @@ struct operator_traits typedef const libint2::TennoGmEval core_eval_type; }; +template <> +struct operator_traits + : public detail::default_operator_traits { + /// Cutoff in the same length units as Shell::O. Zero gives zero integrals; + /// positive infinity gives ordinary Coulomb integrals. + typedef scalar_type oper_params_type; + static oper_params_type default_params() { return oper_params_type{0}; } + typedef const libint2::os_core_ints::truncated_coulomb_gm_eval + core_eval_type; +}; + /// the runtime version of \c operator_traits::default_params() libint2::any default_params(const Operator& oper); diff --git a/include/libint2/engine.impl.h b/include/libint2/engine.impl.h index e1a90f254..f678aee95 100644 --- a/include/libint2/engine.impl.h +++ b/include/libint2/engine.impl.h @@ -101,7 +101,8 @@ typename std::remove_all_extents::type* to_ptr1(T (&a)[N]) { (coulomb_op, /* coulomb_op */ \ (opop_coulomb_opop, /* opop_coulomb_opop */ \ (op_coulomb_op, /* op_coulomb_op */ \ - BOOST_PP_NIL)))))))))))))))))))))))))))))) + (eri, /* truncated_coulomb */ \ + BOOST_PP_NIL))))))))))))))))))))))))))))))) #define BOOST_PP_NBODY_OPERATOR_INDEX_TUPLE \ BOOST_PP_MAKE_TUPLE(BOOST_PP_LIST_SIZE(BOOST_PP_NBODY_OPERATOR_LIST)) @@ -980,7 +981,14 @@ __libint2_engine_inline void Engine::init_core_ints_params(const any& params) { "): geminal exponent " + std::to_string(p.first) + " is negative"); }; - if (oper_ == Operator::delcgtg2) { + if (oper_ == Operator::truncated_coulomb) { + const auto cutoff = any_cast< + operator_traits::oper_params_type>(params); + if (!(cutoff >= 0)) + throw std::invalid_argument( + "truncated Coulomb cutoff must be nonnegative"); + core_ints_params_ = params; + } else if (oper_ == Operator::delcgtg2) { // [g12,[- \Del^2, g12] = 2 (\Del g12) \cdot (\Del g12) // (\Del exp(-a r_12^2) \cdot (\Del exp(-b r_12^2) = 4 a b (r_{12}^2 exp(- // (a+b) r_{12}^2) ) @@ -1682,6 +1690,15 @@ __libint2_engine_inline const Engine::target_ptr_vec& Engine::compute2( .first(); core_eval_ptr->eval(gm_ptr, rho, T, mmax); } break; + case Operator::truncated_coulomb: { + auto& pack = any_cast&>(core_eval_pack_); + const auto cutoff = any_cast::oper_params_type>( + core_ints_params_); + pack.first()->eval(gm_ptr, rho, T, mmax, cutoff, + pack.second()); + } break; case Operator::erf_coulomb: { const auto& core_eval_ptr = any_cast +#include +#include +#include +#include +#include + +namespace libint2 { +namespace detail { +namespace truncated_coulomb { + +struct Workspace; // defined in truncated_coulomb.cc + +/// Mutable work storage of os_core_ints::truncated_coulomb_gm_eval: series +/// coefficients and recurrence buffers reused between calls. Not thread-safe; +/// Engine keeps one per Engine object. +class Scratch { + public: + explicit Scratch(int mmax); + ~Scratch(); + Scratch(const Scratch&); + Scratch& operator=(const Scratch&); + Scratch(Scratch&&) noexcept; + Scratch& operator=(Scratch&&) noexcept; + + private: + std::unique_ptr workspace_; + friend void evaluate(double*, double, double, int, double, Scratch&); +}; + +/// Computes \f$ G_m(R,T) \f$, m = 0..mmax, for the kernel +/// \f$ \theta(R - r)/r \f$ in units of \f$ \rho^{-1/2} \f$: +/// R = 0 gives zeros, R = infinity the Boys function \f$ F_m(T) \f$. +/// @param[out] Gm array of at least mmax+1 values +/// @param[in] precision target absolute error, >= double epsilon; the +/// F_m and zero limits are used only when exact to precision +/// relative to F_m +/// @throw std::invalid_argument for R < 0, NaN, or T < 0 or not finite +void evaluate(double* Gm, double R, double T, int mmax, double precision, + Scratch& scratch); +} // namespace truncated_coulomb +} // namespace detail + +namespace os_core_ints { +/// Interpolation-free core integral evaluator for the truncated Coulomb +/// kernel \f$ \theta(R_c - r_{12})/r_{12} \f$. All arithmetic is binary64. +/// The evaluator is immutable and can be shared between threads; mutable +/// storage is in a scratch_type object, one per thread. +template +class truncated_coulomb_gm_eval { + public: + typedef Real value_type; + typedef detail::truncated_coulomb::Scratch scratch_type; + + /// @param precision target absolute precision; values below double + /// epsilon are raised to it + truncated_coulomb_gm_eval(int mmax, Real precision) + : mmax_(mmax), + precision_(std::max(static_cast(precision), + std::numeric_limits::epsilon())) { + if (mmax < 0 || !std::isfinite(precision_)) + throw std::invalid_argument("invalid truncated Coulomb evaluator setup"); + } + + static std::shared_ptr instance( + int mmax, Real precision = std::numeric_limits::epsilon()) { + return std::make_shared(mmax, precision); + } + + /// @return the maximum m for which \f$ G_m \f$ can be computed + int max_m() const { return mmax_; } + /// @return the target absolute precision + Real precision() const { return static_cast(precision_); } + + /// @param[out] Gm array of at least mmax+1 values + /// @param[in] rho \f$ \rho = \zeta \eta / (\zeta + \eta) \f$ + /// @param[in] T \f$ \rho |PQ|^2 \f$ + /// @param[in] cutoff \f$ R_c \geq 0 \f$, in length units; infinity gives + /// Coulomb + void eval(Real* Gm, Real rho, Real T, int mmax, Real cutoff, + scratch_type& scratch) const { + if (mmax < 0 || mmax > mmax_ || !(rho > 0) || !std::isfinite(rho) || + !(cutoff >= 0)) + throw std::invalid_argument( + "invalid truncated Coulomb evaluator request"); + const double R = + static_cast(cutoff) * std::sqrt(static_cast(rho)); + evaluate(Gm, R, static_cast(T), mmax, scratch); + } + /// same as above, with temporary scratch (allocates on every call) + void eval(Real* Gm, Real rho, Real T, int mmax, Real cutoff) const { + scratch_type scratch(mmax); + eval(Gm, rho, T, mmax, cutoff, scratch); + } + + private: + int mmax_; + double precision_; + + void evaluate(double* Gm, double R, double T, int mmax, + scratch_type& scratch) const { + detail::truncated_coulomb::evaluate(Gm, R, T, mmax, precision_, scratch); + } + template + void evaluate(Output* Gm, double R, double T, int mmax, + scratch_type& scratch) const { + std::vector values(mmax + 1); + detail::truncated_coulomb::evaluate(values.data(), R, T, mmax, precision_, + scratch); + std::copy(values.begin(), values.end(), Gm); + } +}; +} // namespace os_core_ints +} // namespace libint2 +#endif diff --git a/src/lib/libint/CMakeLists.txt b/src/lib/libint/CMakeLists.txt index 8bd73f193..0ebaacabd 100644 --- a/src/lib/libint/CMakeLists.txt +++ b/src/lib/libint/CMakeLists.txt @@ -56,6 +56,7 @@ file(GLOB_RECURSE _populate_inputs CONFIGURE_DEPENDS "${PROJECT_SOURCE_DIR}/doc/progman/*" "${PROJECT_SOURCE_DIR}/cmake/modules/*.cmake" "${CMAKE_CURRENT_SOURCE_DIR}/*.h" + "${CMAKE_CURRENT_SOURCE_DIR}/truncated_coulomb*.cc" ) list(FILTER _populate_inputs EXCLUDE REGEX "/\\.[^/]*$") list(APPEND _populate_inputs diff --git a/src/lib/libint/populate.cmake b/src/lib/libint/populate.cmake index e1a4491f5..8fa5b25c2 100644 --- a/src/lib/libint/populate.cmake +++ b/src/lib/libint/populate.cmake @@ -198,6 +198,13 @@ file( # ==== /src =================================================================== +# Truncated-Coulomb kernel and the generator of its independent test data. +file(GLOB truncated_coulomb_sources + "${PROJECT_SOURCE_DIR}/src/lib/libint/truncated_coulomb*.cc") +file(INSTALL ${truncated_coulomb_sources} DESTINATION "${EXPORT_STAGE_DIR}/src") +file(INSTALL "${PROJECT_SOURCE_DIR}/tools/generate_truncated_coulomb_reference.py" + DESTINATION "${EXPORT_STAGE_DIR}/tools") + file( INSTALL "${PROJECT_SOURCE_DIR}/export/engine.cpp" diff --git a/src/lib/libint/truncated_coulomb.cc b/src/lib/libint/truncated_coulomb.cc new file mode 100644 index 000000000..110a73a6d --- /dev/null +++ b/src/lib/libint/truncated_coulomb.cc @@ -0,0 +1,1393 @@ +/* Binary64 truncated-Coulomb numerical kernel. + * SPDX-License-Identifier: LGPL-3.0-or-later + */ +#include +#include + +// The error bounds below rely on IEEE semantics: infinities as sentinels, +// compensated sums, exact products and directed outward rounding. The build +// compiles this file with value-safe options (see CMakeLists.txt.export); +// fail here rather than return wrong integrals if they were overridden. +#if defined(__FAST_MATH__) || \ + (defined(__FINITE_MATH_ONLY__) && __FINITE_MATH_ONLY__) +#error "truncated_coulomb.cc must be compiled without -ffast-math" +#endif + +#include +#include +#include +#include +#include +#include +#include + +namespace libint2 { +namespace detail { +namespace truncated_coulomb { +namespace { +constexpr double unit = std::numeric_limits::epsilon() / 2; +// Target precision and resource limits of a request; exhausting a limit +// throws std::runtime_error. +struct Options { + // C++11 aggregates have no default member initializers + explicit Options(double precision) : precision(precision) {} + double precision; + unsigned max_extra_orders = 4096; + unsigned max_taylor_terms = 4096; +}; +constexpr double infinity = std::numeric_limits::infinity(); +// Fused multiply-add is used only where the compiler emits it as an +// instruction. Elsewhere std::fma is a library call, emulated in software on +// CPUs without FMA (~150 ns on Ivy Bridge), so each use falls back to plain +// arithmetic with the extra rounding accounted for. Define +// LIBINT2_TRUNCATED_COULOMB_FAST_FMA=0 to test the fallback on FMA hardware. +#ifndef LIBINT2_TRUNCATED_COULOMB_FAST_FMA +#ifdef FP_FAST_FMA +#define LIBINT2_TRUNCATED_COULOMB_FAST_FMA 1 +#else +#define LIBINT2_TRUNCATED_COULOMB_FAST_FMA 0 +#endif +#endif +constexpr bool fast_fma = LIBINT2_TRUNCATED_COULOMB_FAST_FMA; +constexpr unsigned taylor_roundings = fast_fma ? 1 : 2; +// a*b+c, rounded once with fast FMA and twice otherwise; for arithmetic whose +// error budget does not depend on the number of roundings. +inline double mul_add(double a, double b, double c) { + return fast_fma ? std::fma(a, b, c) : a * b + c; +} +// Preserve independent vector lanes on portable targets as well. A scalar +// libm fma call would serialize them; the budget accounts for both roundings. +inline double taylor_add(double weight, double origin, double sum) { + return mul_add(weight, origin, sum); +} +// Numeric outward steps for nonnegative binary64 bounds. Exceptional and +// subnormal results retain nextafter's floating-point exception behavior. +inline double outward(double x, double direction) { + static_assert(sizeof(double) == sizeof(std::uint64_t) && + std::numeric_limits::is_iec559 && + std::numeric_limits::digits == 53, + "truncated Coulomb requires IEEE binary64 arithmetic"); + std::uint64_t bits; + std::memcpy(&bits, &x, sizeof(bits)); + // Integer classification also leaves quiet NaNs quiet until nextafter. + constexpr std::uint64_t min_normal = UINT64_C(0x0010000000000000); + constexpr std::uint64_t max_finite = UINT64_C(0x7fefffffffffffff); + if (bits >= min_normal && bits < max_finite && + (direction == infinity || bits > min_normal)) { + if (direction == infinity) + ++bits; + else + --bits; + std::memcpy(&x, &bits, sizeof(x)); + return x; + } + return std::nextafter(x, direction); +} +// Bound on a*b+c for a*b, c >= 0, rounded toward direction (0 or infinity): +// one outward step per rounding. +inline double outward_mul_add(double a, double b, double c, double direction) { + if (fast_fma) return outward(std::fma(a, b, c), direction); + return outward(outward(a * b, direction) + c, direction); +} +// T - R^2 with the product formed exactly. Near R = sqrt(T) the result is the +// small displacement sqrt(T) - R times sqrt(T) + R, so rounding R^2 first +// would cost a relative error of order T / (T - R^2). +inline double difference_of_square(double T, double R) { + if (fast_fma) return std::fma(-R, R, T); + // Dekker's exact product R*R = p + e via Veltkamp splitting of R. + const double p = R * R, c = 134217729. * R, hi = c - (c - R), lo = R - hi; + const double e = ((hi * hi - p) + 2 * hi * lo) + lo * lo; + return (T - p) - e; +} +// Double-word arithmetic: two_sum and two_product are error-free (Knuth, +// Dekker). Pairs are left unnormalized; callers bound their roundings. +struct Pair { + double hi, lo; +}; +inline Pair two_sum(double a, double b) { + const double s = a + b, t = s - a; + return {s, (a - (s - t)) + (b - t)}; +} +inline Pair fast_two_sum(double a, double b) { // |a| >= |b| + const double s = a + b; + return {s, b - (s - a)}; +} +// Veltkamp splitting, exact for |a| < 2^996. +inline Pair split(double a) { + const double c = 134217729. * a, hi = c - (c - a); + return {hi, a - hi}; +} +// a*b - p for p = fl(a*b), from split operands. +inline double product_error(Pair a, Pair b, double p) { + return ((a.hi * b.hi - p) + a.hi * b.lo + a.lo * b.hi) + a.lo * b.lo; +} +inline Pair two_product(double a, double b) { + const double p = a * b; + if (fast_fma) return {p, std::fma(a, b, -p)}; + return {p, product_error(split(a), split(b), p)}; +} +struct Sum { + double value = 0, correction = 0, magnitude = 0; + void add(double term) { + const double adjusted = term - correction; + const double next = value + adjusted; + correction = (next - value) - adjusted; + value = next; + magnitude += std::abs(term); + } +}; +void validate(double R, double T, int mmax, const Options& options) { + if (!(R >= 0) || !(T >= 0) || !std::isfinite(T) || mmax < 0 || + !(options.precision >= std::numeric_limits::epsilon()) || + !std::isfinite(options.precision)) + throw std::invalid_argument("invalid truncated Coulomb request"); +} +double root_difference(double R, double T, double s) { + if (R >= .5 * s && R <= 2 * s) return difference_of_square(T, R) / (s + R); + return s - R; +} +struct TaylorPlan { + int degree; + double remainder; +}; +struct Candidate { + double value; + double error; +}; +int zero_seed(double T, int n, double tolerance, unsigned limit) { + double product = 1; + const unsigned extra = std::min(limit, static_cast(INT_MAX - n)); + for (unsigned j = 1; j <= extra; ++j) { + const int M = n + static_cast(j); + product = outward(product * outward(T / (M - .5), infinity), infinity); + const double ratio = outward(T / (M + 1.5), infinity); + const double denominator = outward(1 - ratio, 0.); + const double magnitude = + outward(outward(2 / (2. * M + 1), infinity) / denominator, infinity); + const double error = outward(product * magnitude, infinity); + if (error <= tolerance) return M; + } + throw std::runtime_error("truncated Coulomb: zero-seed depth exhausted"); +} +} // namespace + +// Budgets contain the two row sums of the contractive propagation branches. +struct QTolerances { + int join = 0; + double upward = infinity, budget = 0, downward_h = 0; + double operator[](unsigned m) const { + if (m < static_cast(join)) return upward; + if (m == static_cast(join)) return infinity; + return budget * ((2. * m + 1) / downward_h); + } +}; + +// All mutable storage belongs to the caller and is reused between requests. +struct Workspace { + std::vector q, origin, hermite, weights, tails; + std::vector pair_terms; // origin_pair_taylor + std::vector inverse_integer, inverse_half; + QTolerances tolerance; + std::vector rhs, correction; + std::vector> matrix; + std::vector pivots; + std::shared_ptr> boys; + int boys_failed_order = INT_MAX; + double g0 = 0, seed = 0; + int upward_last = 0, seed_order = -1; + int origin_last = -1; + double origin_difference = 0, origin_envelope = 0; + double exp_minus_t = -1; + double sqrt_t = -1, displacement = 0; + bool displacement_ready = false; + bool gaussian_ready = false; + double gaussian_minus = 0, gaussian_plus = 0; + int limit_order = -1, limit_power = 0; + double limit_precision = 0, limit_gap = 0, limit_min_t = 0; + void reset_request() { + origin_last = -1; + exp_minus_t = -1; + sqrt_t = -1; + displacement_ready = false; + gaussian_ready = false; + weights.clear(); + tails.clear(); + } + void ensure_reciprocals(std::size_t size) { + if (inverse_integer.size() >= size) return; + static const ExpensiveNumbers numbers; + const auto begin = inverse_integer.size(); + inverse_integer.resize(size); + inverse_half.resize(size); + for (std::size_t k = begin; k < size; ++k) { + inverse_integer[k] = k ? 1. / k : 0; + inverse_half[k] = k < 128 ? 2 * numbers.twoi1[k] : 1 / (k + .5); + } + } + explicit Workspace(int n) { + const std::size_t capacity = static_cast(n) + 65; + for (auto* v : {&q, &origin, &hermite, &rhs, &correction, &weights, &tails}) + v->reserve(capacity); + ensure_reciprocals(capacity); + matrix.reserve(capacity); + pivots.reserve(capacity); + // The standard Engine order range uses exactly the ordinary Coulomb + // evaluator. Higher radial orders retain an uncapped binary64 limit. + boys = FmEval_Chebyshev7::instance(0); + try { + if (n > boys->max_m()) boys = FmEval_Chebyshev7::instance(n); + } catch (const std::invalid_argument&) { + boys_failed_order = n; + } + } +}; +Scratch::Scratch(int n) { + if (n < 0) + throw std::invalid_argument("negative truncated Coulomb scratch order"); + workspace_.reset(new Workspace(n)); +} +Scratch::~Scratch() = default; +Scratch::Scratch(const Scratch& other) + : workspace_(other.workspace_ ? new Workspace(*other.workspace_) + : nullptr) {} +Scratch& Scratch::operator=(const Scratch& other) { + if (this != &other) { + if (other.workspace_) { + if (workspace_) + *workspace_ = *other.workspace_; + else + workspace_.reset(new Workspace(*other.workspace_)); + } else + workspace_.reset(); + } + return *this; +} +Scratch::Scratch(Scratch&&) noexcept = default; +Scratch& Scratch::operator=(Scratch&&) noexcept = default; + +namespace { +enum class CutoffLimit { none, zero, boys }; + +// Upper bound on -log F_n(T), hence on -log F_m(T) for all m <= n. For +// T >= n+1/2 the regularized gamma P(n+1/2,T) >= 1/2 (a gamma median is +// below its mean), so F_n >= Gamma(n+1/2)/(4 T^(n+1/2)); otherwise +// F_n >= exp(-T)/(2n+1). The added 1 absorbs rounding of log and lgamma. +double boys_log_scale(double T, int n) { + const double scale = + T >= n + .5 ? (n + .5) * std::log(T) - std::lgamma(n + .5) + std::log(4.) + : T + std::log(2. * n + 1); + return std::max(0., scale) + 1; +} + +// The gap for an error bound exp(-gap^2) relative to F_n(T). +double relative_gap(double gap, double T, int n) { + const double square = + outward_mul_add(gap, gap, boys_log_scale(T, n), infinity); + return outward(std::sqrt(square), infinity); +} + +CutoffLimit finite_cutoff_limit(double R, double T, int n, + const Options& options, Workspace& w) { + // Adjacent binary64 bounds on log(2); no runtime transcendental bound is + // needed. frexp chooses an exact dyadic error budget <= precision / 32. + constexpr double log_two_lo = .6931471805599453; + constexpr double log_two_hi = .6931471805599454; + if (w.limit_order != n || w.limit_precision != options.precision) { + int exponent; + std::frexp(options.precision / 32, &exponent); + w.limit_power = std::max(0, 1 - exponent); + const double square = + outward_mul_add(double(w.limit_power), log_two_hi, 2. * n, infinity); + w.limit_gap = outward(std::sqrt(square), infinity); + w.limit_min_t = + outward_mul_add(double(w.limit_power), log_two_lo, double(n), 0.); + w.limit_order = n; + w.limit_precision = options.precision; + } + // For a=sqrt(T+n), Cauchy's derivative bound on the disk of radius n + // gives |F_m-G_m| <= exp(2n-(R-a)^2) when R>=a, simultaneously for + // m<=n. For n=0 the same inequality follows from the real Gaussian + // integral. Bound the distance downward and the required gap upward. + // Both limits are accepted only below precision * F_n(T), not precision + // alone: at large T the higher orders are far below precision, and an + // absolute test would replace them by values without a correct digit. + if (R >= w.limit_gap) { + const double center = outward(T + n, infinity); + const double a = outward(std::sqrt(center), infinity); + if (R >= a) { + const double distance = outward(R - a, 0.); + if (distance >= w.limit_gap && + distance >= relative_gap(w.limit_gap, T, n)) + return CutoffLimit::boys; + } + } + // In the other direction |G_m| <= R^2 exp(2n-(a-R)^2), R<=a. + // R^2 <= 1 for R<=1; otherwise frexp gives R^2 < 2^(2*exponent). + // T >= n + power*log(2) is necessary for this conservative test. + if (T >= w.limit_min_t) { + const double center = outward(T + n, 0.); + const double a = outward(std::sqrt(center), 0.); + if (a >= R) { + const double distance = outward(a - R, 0.); + if (distance >= w.limit_gap) { + double gap = w.limit_gap; + if (R > 1) { + int exponent; + std::frexp(R, &exponent); + const double square = + outward_mul_add(double(w.limit_power + 2 * exponent), log_two_hi, + 2. * n, infinity); + gap = outward(std::sqrt(square), infinity); + } + if (distance >= gap && distance >= relative_gap(gap, T, n)) + return CutoffLimit::zero; + } + } + } + return CutoffLimit::none; +} + +// Cache coefficients and their tail bounds within one request. G0 and Q +// retain their own first admissible degree, even when tolerances differ. +TaylorPlan taylor_plan(double T, double tolerance, unsigned max_terms, + Workspace& w) { + if (T >= max_terms) return {-1, infinity}; + if (w.weights.empty()) { + w.weights.push_back(1); + w.tails.push_back(infinity); + } + const double target = tolerance / 128; + const auto found = std::lower_bound(w.tails.begin() + 1, w.tails.end(), + target, std::greater()); + if (found != w.tails.end()) + return {static_cast(found - w.tails.begin()) - 1, *found}; + if (!std::isfinite(w.weights.back())) return {-1, infinity}; + for (unsigned k = static_cast(w.weights.size()); k <= max_terms; + ++k) { + const double weight = w.weights.back() * (-T / k); + w.weights.push_back(weight); + const double ratio = T / (k + 1.); + const double tail = ratio < 1 ? std::abs(weight) / (1 - ratio) : infinity; + w.tails.push_back(tail); + if (!std::isfinite(weight)) return {-1, infinity}; + if (tail <= target) return {static_cast(k) - 1, tail}; + } + return {-1, infinity}; +} +// The first omitted G_m coefficient is bounded by +// min(R^2, 2 / (2m + 2k + 1)). Reuse the Q Taylor weights and geometric +// tails, while choosing G's degree from its own, smaller remainder. +TaylorPlan taylor_g_plan(double x, double T, unsigned m, double tolerance, + unsigned max_terms, Workspace& w) { + if (T >= max_terms) return {-1, infinity}; + if (w.weights.empty()) { + w.weights.push_back(1); + w.tails.push_back(infinity); + } + const double target = tolerance / 128; + const auto tail_at = [&](std::size_t k) { + return w.tails[k] * std::min(x, 2 / (2. * m + 2. * k + 1)); + }; + // The cached tails and coefficient bound are both nonincreasing. Find + // the first admissible degree without rescanning G0's complete prefix. + std::size_t lo = 1, hi = std::min(w.tails.size(), + static_cast(max_terms) + 1); + const std::size_t end = hi; + while (lo < hi) { + const auto mid = lo + (hi - lo) / 2; + if (tail_at(mid) <= target) + hi = mid; + else + lo = mid + 1; + } + if (lo < end) return {static_cast(lo) - 1, tail_at(lo)}; + if (!std::isfinite(w.weights.back())) return {-1, infinity}; + for (unsigned k = static_cast(w.weights.size()); k <= max_terms; + ++k) { + const double weight = w.weights.back() * (-T / k); + w.weights.push_back(weight); + const double ratio = T / (k + 1.); + w.tails.push_back(ratio < 1 ? std::abs(weight) / (1 - ratio) : infinity); + if (!std::isfinite(weight)) return {-1, infinity}; + const double tail = tail_at(k); + if (tail <= target) return {static_cast(k) - 1, tail}; + } + return {-1, infinity}; +} +double exponential(double T, Workspace& w) { + if (w.exp_minus_t < 0) w.exp_minus_t = std::exp(-T); + return w.exp_minus_t; +} +double square_root(double T, Workspace& w) { + if (w.sqrt_t < 0) w.sqrt_t = std::sqrt(T); + return w.sqrt_t; +} +double displacement(double R, double T, Workspace& w) { + if (!w.displacement_ready) { + w.displacement = root_difference(R, T, square_root(T, w)); + w.displacement_ready = true; + } + return w.displacement; +} +void gaussians(double R, double T, Workspace& w) { + if (w.gaussian_ready) return; + const double z = displacement(R, T, w), plus = R + square_root(T, w); + w.gaussian_minus = std::exp(-z * z); + w.gaussian_plus = std::exp(-plus * plus); + w.gaussian_ready = true; +} +void boys_limit(double T, int n, Workspace& w, double* g) { + // A workspace reserved beyond the native table can later receive a lower + // request within it. Select by the actual order, independently of reserve + // capacity. Remember failed growth so unsupported orders do not keep + // throwing. + if (n > w.boys->max_m() && n < w.boys_failed_order) { + try { + w.boys = FmEval_Chebyshev7::instance(n); + } catch (const std::invalid_argument&) { + w.boys_failed_order = n; + } + } + if (w.boys && n <= w.boys->max_m()) { + w.boys->eval(g, T, n); + } else { + // Beyond the existing Boys table, join two contractive recurrences. + // In particular, an underflowed exp(-T) must not erase F0, and an + // upward recurrence must not grow a parasitic solution above its turn. + w.boys->eval(g, T, 0); + const double e = std::exp(-T); + int upward = 0; + while (upward < n && upward + .5 < T) { + g[upward + 1] = mul_add(2. * upward + 1, g[upward], -e) / (2 * T); + ++upward; + } + if (upward < n) { + double value = 0; + if (e != 0) { + double term = 1 / (2. * n + 1); + Sum sum; + sum.add(term); + for (int k = 1;; ++k) { + term *= T / (n + k + .5); + sum.add(term); + if (term <= unit * sum.value) break; + } + value = e * sum.value; + } + g[n] = value; + for (int m = n - 1; m > upward; --m) { + value = mul_add(2 * T, value, e) / (2. * m + 1); + g[m] = value; + } + } + } +} +void evaluate_limit(double T, int n, CutoffLimit limit, Workspace& w, + double* g) { + if (limit == CutoffLimit::boys) + boys_limit(T, n, w, g); + else + std::fill(g, g + static_cast(n) + 1, 0); +} +void origin_q(double R, unsigned last, Workspace& w) { + if (w.origin_last >= static_cast(last)) return; + w.origin.resize(static_cast(last) + 1); + const double x = R * R; + if (!std::isfinite(x)) { + std::fill(w.origin.begin(), w.origin.end(), 0); + w.origin_envelope = 0; + w.origin_last = static_cast(last); + return; + } + if (w.origin_last < 0) { + w.origin[0] = std::exp(-x); + w.origin_envelope = w.origin[0]; + w.origin_difference = 0; + w.origin_last = 0; + } + w.ensure_reciprocals(static_cast(last)); + for (unsigned k = static_cast(w.origin_last); k < last; ++k) { + w.origin_difference = + mul_add(k, w.origin_difference, -x * w.origin[k]) * w.inverse_half[k]; + w.origin[k + 1] = w.origin[k] + w.origin_difference; + w.origin_envelope = std::max(w.origin_envelope, std::abs(w.origin[k + 1])); + } + w.origin_last = static_cast(last); +} +Candidate r_series(double R, double T, unsigned m, bool g, double tolerance, + unsigned max_terms, Workspace& w) { + const double x = R * R, e = exponential(T, w); + if (!std::isfinite(x) || e == 0) return {0, infinity}; + double previous = 0, term = g ? x * e : e, majorant = std::abs(term); + const double alpha1 = g ? m + 1.5 : m + .5; + double ratio = x * (alpha1 + T) / (g ? 3. : .5); + Sum sum; + for (unsigned j = 0; j < max_terms; ++j) { + sum.add(term); + const double next_majorant = majorant * ratio; + const double next_denominator = + g ? (j + 3.) * (j + 2.5) : (j + 2.) * (j + 1.5); + const double next_ratio = x * (2 * j + 2 + alpha1 + T) / next_denominator; + if (!std::isfinite(next_majorant)) return {0, infinity}; + if (next_ratio < 1) { + const double tail = next_majorant / (1 - next_ratio); + if (tail <= tolerance / 128) + return {sum.value, 16 * unit * (j + 1) * sum.magnitude + tail}; + } + // No signed next term is needed after convergence or resource exhaustion. + if (j + 1 == max_terms) break; + const double denominator = g ? (j + 2.) * (j + 1.5) : (j + 1.) * (j + .5); + const double a = -x * (2 * j + alpha1 - T) / denominator; + const double b = j == 0 ? 0 + : g ? (m == 0 ? x * x * j / (denominator * (j + 1.)) + : x * x * j * (j + m + .5) / + (denominator * (j + 1.) * (j + .5))) + : x * x * (j + m - .5) / (denominator * (j - .5)); + const double next = mul_add(a, term, -b * previous); + if (!std::isfinite(next)) return {0, infinity}; + previous = term; + term = next; + majorant = next_majorant; + ratio = next_ratio; + } + return {0, infinity}; +} +Candidate taylor_g0(double R, const TaylorPlan& plan, Workspace& w) { + origin_q(R, static_cast(plan.degree), w); + Sum sum; + sum.add(-std::expm1(-R * R)); + for (int j = 1; j <= plan.degree; ++j) + sum.add(w.weights[j] * (1 - w.origin[j]) / (2. * j + 1)); + return {sum.value, 8 * unit * sum.magnitude + plan.remainder}; +} +Candidate direct_seed(double R, double T, unsigned m, const Options& options, + Workspace& w) { + const double x = R * R; + // The same conditioning bound as G0's radius series, at the seed order. + if (x * ((m + 1.5 + T) / 3) <= .25) { + const auto value = + r_series(R, T, m, true, options.precision, options.max_taylor_terms, w); + if (std::isfinite(value.error)) return value; + } + // The origin coefficients obey 0 <= G_k(R,0) <= 2/(2k+1). + // exp(T) <= 2 bounds the Taylor absolute-term sum, as for G0. + if (T <= std::log(2.)) { + const auto plan = + taylor_g_plan(x, T, m, options.precision, options.max_taylor_terms, w); + if (plan.degree >= 0) { + origin_q(R, m + static_cast(plan.degree), w); + Sum sum; + for (int j = 0; j <= plan.degree; ++j) + sum.add(w.weights[j] * (1 - w.origin[m + j]) / (2. * m + 2. * j + 1)); + if (std::isfinite(sum.value)) + return {sum.value, 8 * unit * sum.magnitude + plan.remainder}; + } + } + return {0, infinity}; +} +Candidate evaluate_g0(double R, double T, const Options& options, + Workspace& w) { + if (T == 0) { + const double value = -std::expm1(-R * R); + return {value, 4 * unit * std::abs(value)}; + } + const double x = R * R; + // r_0 <= 1/4 bounds the positive majorant / |sum| by two. + if (x * ((T + 1.5) / 3) <= .25) { + const auto value = + r_series(R, T, 0, true, options.precision, options.max_taylor_terms, w); + if (std::isfinite(value.error)) return value; + } + // For exp(T) <= 2 the origin Taylor series has a bounded absolute-term + // sum. It also avoids dividing rounded erf(sqrt(T)) by a tiny sqrt(T). + if (T <= std::log(2.)) { + const auto plan = + taylor_plan(T, options.precision, options.max_taylor_terms, w); + if (plan.degree >= 0) return taylor_g0(R, plan, w); + } + const double s = square_root(T, w), z = displacement(R, T, w); + const double factor = .44311346272637900682 / s; + double value, magnitude; + if (R >= s) { + double f; + w.boys->eval(&f, T, 0); + const double a = factor * std::erfc(-z), b = factor * std::erfc(R + s); + value = f - a + b; + magnitude = f + a + b; + // Near R = s this form cancels by up to 3. With F0 = 2 factor erf(s) + // the difference form below is the same function, and it is well + // conditioned once erfc(s) is small. + if (magnitude > 2 * std::abs(value)) { + const double c = std::erfc(z), d = std::erfc(s + R), e = 2 * std::erfc(s); + if (factor * (c + d + e) < magnitude) { + value = factor * (c + d - e); + magnitude = factor * (c + d + e); + } + } + } else { + const double a = std::erfc(z), b = std::erfc(s + R), c = 2 * std::erfc(s); + value = factor * (a + b - c); + magnitude = factor * (a + b + c); + } + if (magnitude <= 2 * std::abs(value) || magnitude <= options.precision / 16) + return {value, 16 * unit * magnitude}; + const auto radius = + r_series(R, T, 0, true, options.precision, options.max_taylor_terms, w); + if (std::isfinite(radius.error) && radius.error < 16 * unit * magnitude) + return radius; + const auto plan = + taylor_plan(T, options.precision, options.max_taylor_terms, w); + if (plan.degree >= 0) { + const auto taylor = taylor_g0(R, plan, w); + if (std::isfinite(taylor.error) && taylor.error < 16 * unit * magnitude) + return taylor; + } + return {value, 16 * unit * magnitude}; +} + +void gaussian_q0(double R, double T, Workspace& w) { + gaussians(R, T, w); + // This is also Miller's known boundary, so solving cannot improve it. + w.q[0] = .5 * (w.gaussian_minus + w.gaussian_plus); +} +void scaled_hermites(double R, double T, unsigned last, Workspace& w) { + auto& p = w.hermite; + p.resize(static_cast(last) + 1); + gaussians(R, T, w); + const bool minus_active = w.gaussian_minus != 0; + const bool plus_active = w.gaussian_plus != 0; + p[0] = w.gaussian_minus + w.gaussian_plus; + if (!minus_active && !plus_active) { + std::fill(p.begin(), p.end(), 0); + return; + } + if (!last) return; + const double s = square_root(T, w); + const double a_minus = minus_active ? displacement(R, T, w) / s : 0; + const double a_plus = plus_active ? (s + R) / s : 0; + double minus_previous = w.gaussian_minus, plus_previous = w.gaussian_plus; + double minus = a_minus * minus_previous, plus = a_plus * plus_previous; + p[1] = minus + plus; + const double inv2T = .5 / T; + for (unsigned k = 1; k < last; ++k) { + const double b = k * inv2T; + const double next_minus = + minus_active ? mul_add(a_minus, minus, -b * minus_previous) : 0; + const double next_plus = + plus_active ? mul_add(a_plus, plus, -b * plus_previous) : 0; + p[k + 1] = next_minus + next_plus; + minus_previous = minus; + plus_previous = plus; + minus = next_minus; + plus = next_plus; + } +} + +bool hermite_q(double R, double T, unsigned last, Workspace& w) { + scaled_hermites(R, T, last, w); + w.ensure_reciprocals(static_cast(last) + 1); + // Dividing .25 by T avoids overflowing a denominator 4*T at large T. + const double inv4T = .25 / T; + for (unsigned m = 1; m <= last; ++m) { + if (m == static_cast(w.upward_last)) continue; + Sum sum; + double weight = .5; + for (unsigned k = m; k > 0; --k) { + sum.add(weight * w.hermite[k]); + if (k > 1) + weight *= + (2. * m - k) * (k - 1.) * inv4T * w.inverse_integer[m - k + 1]; + } + if (!std::isfinite(sum.value) || std::abs(sum.value) > 1 || + unit * std::max(0., sum.magnitude - std::abs(sum.value)) > + w.tolerance[m]) + return false; + w.q[m] = sum.value; + } + return true; +} + +double taylor_absolute_weights(const TaylorPlan& plan, double gamma, + const Workspace& w) { + if (plan.degree == 0) return 1; + const double normal = std::numeric_limits::min(); + // |w_j| is unimodal. These endpoint tests ensure normal intermediates + // in the two-operation coefficient recurrence for every j <= degree. + if (std::abs(w.weights[plan.degree]) >= normal && + std::abs(w.weights[1]) >= normal * plan.degree && w.exp_minus_t > 0) { + if (taylor_roundings == 1) { + const double du = (2. * plan.degree) * unit; + gamma = outward(du / (1 - du), infinity); + } + const double denominator = outward(w.exp_minus_t, 0.); + if (denominator > 0) + return outward(outward(1 + gamma, infinity) / denominator, infinity); + } + // Subnormal coefficients need an absolute bound on the stored prefix. + double result = 1; + for (int j = 1; j <= plan.degree; ++j) + result = outward(result + std::abs(w.weights[j]), infinity); + return result; +} + +bool taylor_q(double R, unsigned last, const TaylorPlan& plan, Workspace& w) { + origin_q(R, last + static_cast(plan.degree), w); + const unsigned low_end = + std::min(last + 1, static_cast(w.upward_last)); + const unsigned high_begin = + std::min(last + 1, static_cast(w.upward_last) + 1); + unsigned fast_low_end = 1, fast_high_begin = last + 1; + const double operations = taylor_roundings * static_cast(plan.degree); + const double du = operations * unit; + const double gamma = du == 0 ? 0 : outward(du / (1 - du), infinity); + double largest_tolerance = low_end > 1 ? w.tolerance[1] : 0; + if (high_begin <= last) + largest_tolerance = std::max(largest_tolerance, w.tolerance[last]); + // A cheap exp(-T) predictor avoids constructing an accumulation bound + // when no requested order could use ordinary summation. The final test + // below uses the actual stored weights and origin envelope, rather than + // treating rounded coefficients as exact mathematical coefficients. + if (gamma * w.origin_envelope <= + (largest_tolerance - plan.remainder) * w.exp_minus_t) { + const double absolute_weights = taylor_absolute_weights(plan, gamma, w); + const double magnitude_bound = + outward(absolute_weights * w.origin_envelope, infinity); + const double underflow = outward( + operations * std::numeric_limits::denorm_min() / (1 - du), + infinity); + const double accumulation_error = outward( + outward(gamma * magnitude_bound, infinity) + underflow, infinity); + const double error = outward(accumulation_error + plan.remainder, infinity); + if (low_end > 1 && error <= w.tolerance[1]) fast_low_end = low_end; + // Downward tolerances increase with m, so the eligible orders form a + // single contiguous suffix. The unused join separates the two ranges. + unsigned lo = high_begin, hi = last + 1; + while (lo < hi) { + const unsigned mid = lo + (hi - lo) / 2; + if (error <= w.tolerance[mid]) + hi = mid; + else + lo = mid + 1; + } + fast_high_begin = lo; + } + std::copy(w.origin.begin() + 1, w.origin.begin() + fast_low_end, + w.q.begin() + 1); + std::copy(w.origin.begin() + fast_high_begin, w.origin.begin() + last + 1, + w.q.begin() + fast_high_begin); + // Independent m lanes permit vector multiply/add, or native vector FMA. + if (fast_low_end > 1 || fast_high_begin <= last) { + for (int j = 1; j <= plan.degree; ++j) { + const double weight = w.weights[j]; + for (unsigned m = 1; m < fast_low_end; ++m) + w.q[m] = taylor_add(weight, w.origin[m + j], w.q[m]); + for (unsigned m = fast_high_begin; m <= last; ++m) + w.q[m] = taylor_add(weight, w.origin[m + j], w.q[m]); + } + } + for (unsigned m = 1; m <= last; ++m) { + if (m == static_cast(w.upward_last)) continue; + if (m < fast_low_end || m >= fast_high_begin) { + if (!std::isfinite(w.q[m]) || std::abs(w.q[m]) > 1 + unit * 2) + return false; + continue; + } + Sum sum; + sum.add(w.origin[m]); + for (int j = 1; j <= plan.degree; ++j) { + sum.add(w.weights[j] * w.origin[m + j]); + } + if (!std::isfinite(sum.value) || std::abs(sum.value) > 1 + unit * 2) + return false; + w.q[m] = sum.value; + } + return true; +} +// Extending the artificial boundary reuses the outward recurrence prefix. +// At an accepted endpoint its product has already included that order, +// while the interval state still precedes the next recurrence step. +struct MillerBoundary { + double T, xlo, xhi; + unsigned last, order = 0; + double qlo = 0, qhi = 0, hlo = 0, hhi = 0, product = 1; + bool counted = false, failed = false; + MillerBoundary(double R, double t, unsigned l) : T(t), last(l) { + const double x = R * R; + xlo = outward(x, 0.); + xhi = outward(x, infinity); + } + double advance(unsigned n) { + if (failed) return infinity; + for (;;) { + if (!counted && order >= last) { + const double denominator = outward(qhi + 1, 0.); + const double ratio = std::min(1., outward(qhi / denominator, infinity)); + product = outward(product * ratio, infinity); + } + counted = true; + if (order == n) return product; + const double flo = outward_mul_add(double(order), hlo, + outward_mul_add(xlo, qlo, .5, 0.), 0.); + const double fhi = + outward_mul_add(double(order), hhi, + outward_mul_add(xhi, qhi, .5, infinity), infinity); + const double dlo = outward(T + flo, 0.), dhi = outward(T + fhi, infinity); + qlo = outward(outward(T * outward(qlo + 1, 0.), 0.) / dhi, 0.); + qhi = outward(outward(T * outward(qhi + 1, infinity), infinity) / dlo, + infinity); + // Retain the shared f in the monotone ratio f/(T+f). + hlo = outward(flo / outward(T + flo, infinity), 0.); + hhi = outward(fhi / outward(T + fhi, 0.), infinity); + if (!std::isfinite(qhi) || !std::isfinite(hhi)) { + failed = true; + return infinity; + } + ++order; + counted = false; + } + } +}; + +bool miller_q(double R, double T, unsigned last, double tolerance, + unsigned max_extra, Workspace& w) { + if (T == 0 || !std::isfinite(R * R)) return false; + gaussians(R, T, w); + const double q0 = .5 * (w.gaussian_minus + w.gaussian_plus); + if (q0 == 0) return false; + max_extra = std::min(max_extra, static_cast(INT_MAX) - last); + if (!max_extra) return false; + MillerBoundary boundary_state(R, T, last); + for (unsigned extra = std::min(16u, max_extra); extra > 0;) { + const unsigned n = last + extra; + const double boundary = boundary_state.advance(n); + if (!(boundary <= tolerance / 128)) { + if (extra == max_extra) break; + extra += std::min(extra, max_extra - extra); + continue; + } + // Column j of row i is stored at j-i+2. Fill-in after pivoting extends + // through column i+3; no dense matrix or external linear algebra is used. + auto& matrix = w.matrix; + matrix.resize(n); + auto& rhs = w.rhs; + rhs.resize(n); + auto& pivots = w.pivots; + pivots.resize(n); + for (unsigned m = 0; m < n; ++m) { + matrix[m].fill(0); + rhs[m] = 0; + const double coefficients[] = {-double(m), T + 2 * m + .5 - R * R, + -(2 * T + m + .5), T}; + for (unsigned k = 0; k < 4; ++k) { + const int order = static_cast(m) + static_cast(k) - 1; + if (order == 0) rhs[m] -= coefficients[k]; + if (order > 0 && order <= static_cast(n)) + matrix[m][k] = coefficients[k]; + } + } + bool finite = true; + for (unsigned k = 0; k < n; ++k) { + unsigned pivot = k; + for (unsigned i = k + 1; i < std::min(n, k + 3); ++i) + if (std::abs(matrix[i][k - i + 2]) > + std::abs(matrix[pivot][k - pivot + 2])) + pivot = i; + if (pivot != k) { + for (unsigned j = k; j < std::min(n, k + 4); ++j) + std::swap(matrix[k][j - k + 2], matrix[pivot][j - pivot + 2]); + std::swap(rhs[k], rhs[pivot]); + } + pivots[k] = pivot; + const double diagonal = matrix[k][2]; + if (diagonal == 0 || !std::isfinite(diagonal)) { + finite = false; + break; + } + for (unsigned i = k + 1; i < std::min(n, k + 3); ++i) { + const double factor = matrix[i][k - i + 2] / diagonal; + for (unsigned j = k + 1; j < std::min(n, k + 4); ++j) + matrix[i][j - i + 2] = + mul_add(-factor, matrix[k][j - k + 2], matrix[i][j - i + 2]); + rhs[i] = mul_add(-factor, rhs[k], rhs[i]); + matrix[i][k - i + 2] = factor; + } + } + if (!finite) return false; + for (unsigned kk = n; kk > 0; --kk) { + const unsigned k = kk - 1; + for (unsigned j = k + 1; j < std::min(n, k + 4); ++j) + rhs[k] = mul_add(-matrix[k][j - k + 2], rhs[j], rhs[k]); + rhs[k] /= matrix[k][2]; + } + for (auto& value : rhs) value *= q0; + auto& correction = w.correction; + correction.resize(n); + for (unsigned iteration = 0; iteration < 2; ++iteration) { + // Form the residual in differences, before collecting the diagonal + // coefficient. This retains R^2 when 2m is much larger and preserves + // a nearly constant hierarchy. Refinement uses the same double LU. + for (unsigned m = 0; m < n; ++m) { + const double current = m == 0 ? q0 : rhs[m - 1]; + const double previous = m == 0 ? 0 : m == 1 ? q0 : rhs[m - 2]; + const double next = rhs[m], after = m + 1 == n ? 0 : rhs[m + 1]; + const double difference = next - current; + const double residual = mul_add( + T, (after - next) - difference, + mul_add(double(m), current - previous, -(m + .5) * difference)); + correction[m] = mul_add(R * R, current, -residual); + } + for (unsigned k = 0; k < n; ++k) { + std::swap(correction[k], correction[pivots[k]]); + for (unsigned i = k + 1; i < std::min(n, k + 3); ++i) + correction[i] = + mul_add(-matrix[i][k - i + 2], correction[k], correction[i]); + } + for (unsigned kk = n; kk > 0; --kk) { + const unsigned k = kk - 1; + for (unsigned j = k + 1; j < std::min(n, k + 4); ++j) + correction[k] = + mul_add(-matrix[k][j - k + 2], correction[j], correction[k]); + correction[k] /= matrix[k][2]; + } + for (unsigned k = 0; k < n; ++k) rhs[k] += correction[k]; + } + double envelope = q0; + for (unsigned m = 1; m <= last; ++m) + envelope = std::max(envelope, std::abs(rhs[m - 1])); + if (std::isfinite(envelope)) { + w.q[0] = q0; + for (unsigned m = 1; m <= last; ++m) w.q[m] = rhs[m - 1]; + return true; + } + if (extra == max_extra) break; + extra += std::min(extra, max_extra - extra); + } + return false; +} + +// Uniform injection budgets, scaled by the exact recurrence coefficients. +// Along each contractive branch the coefficients increase in traversal +// order, so its positive row sums increase up to the join. H is the final +// row sum; no per-order tolerance array or maximum reduction is needed. +// Thus sum_j |d G_m / d Q_j| tolerance_j <= budget for every requested m, +// including a join arbitrarily close to an integer-plus-one-half. +double q_tolerances(double T, int last, double budget, Workspace& w) { + auto& tolerance = w.tolerance; + tolerance.join = w.upward_last; + tolerance.budget = budget; + double h = 0; + for (int m = 0; m < w.upward_last; ++m) { + h = mul_add((m + .5) / T, h, 1.); + } + tolerance.upward = h > 0 ? budget * (2 * T / h) : infinity; + h = 0; + for (int m = last; m > w.upward_last; --m) { + h = mul_add(T / (m + .5), h, 1.); + } + tolerance.downward_h = h; + double minimum = tolerance.join > 1 ? tolerance.upward : infinity; + if (last > tolerance.join) + minimum = std::min(minimum, tolerance[tolerance.join + 1]); + return minimum; +} + +void evaluate_q(double R, double T, unsigned last, double tolerance, + const Options& options, Workspace& w) { + w.q.resize(static_cast(last) + 1); + if (T == 0) { + origin_q(R, last, w); + std::copy(w.origin.begin(), w.origin.begin() + last + 1, w.q.begin()); + return; + } + // An unused Q0 needs no shifted Gaussians on either series path. + if (w.upward_last > 0) gaussian_q0(R, T, w); + if (last == 0) return; + const double x = R * R; + if (x * (2 * (last + .5 + T)) <= .25) { + bool finite = true; + for (unsigned m = 1; m <= last; ++m) { + if (m == static_cast(w.upward_last)) continue; + const double target = + std::isfinite(w.tolerance[m]) ? w.tolerance[m] : options.precision; + const auto value = + r_series(R, T, m, false, target, options.max_taylor_terms, w); + if (!std::isfinite(value.error)) { + finite = false; + break; + } + w.q[m] = value.value; + } + if (finite) return; + } + // |Q_m(R,0)| <= 1 bounds the Taylor absolute-term sum by exp(T). + // This predicts conditioning without first evaluating every candidate. + TaylorPlan plan{-1, infinity}; + bool taylor_tried = false; + if (unit <= tolerance * exponential(T, w)) { + taylor_tried = true; + plan = taylor_plan(T, tolerance, options.max_taylor_terms, w); + if (plan.degree >= 0 && taylor_q(R, last, plan, w)) return; + } + if (hermite_q(R, T, last, w)) return; + if (miller_q(R, T, last, tolerance, options.max_extra_orders, w)) return; + // Resource exhaustion or underflow in one representation is a fallback, + // rather than a reason to enumerate all representations on every call. + if (!taylor_tried) { + plan = taylor_plan(T, tolerance, options.max_taylor_terms, w); + if (plan.degree >= 0 && taylor_q(R, last, plan, w)) return; + } + throw std::runtime_error( + "truncated Coulomb: no finite bounded Q representation"); +} + +// Double-word tables of 1/((1/2)_k (2k+1)) and 1/k!. Each entry takes at +// most k double-word scalings by a double, each with a relative error of a +// few u^2 (Joldes, Muller, Popescu, ACM TOMS 44, 2017), and one Newton +// reciprocal. The bound 512 u^2 = 2^-97 used below exceeds the 9 u^2 +// measured against 60-digit values; FMA and Dekker products give the same +// entries. +constexpr int origin_pair_capacity = 160; +constexpr double origin_pair_max_t = 12; +struct OriginPairTables { + std::array coefficient, inverse_factorial; + // log (1/2)_k, and kappa_k = k! / (1/2)_{k+1} for the error estimate + std::array log_pochhammer, kappa; + static constexpr double error = 6.310887241768094e-30; // 2^-97 + OriginPairTables() { + Pair pochhammer{1, 0}, factorial{1, 0}; + for (int k = 0; k < origin_pair_capacity; ++k) { + coefficient[k] = reciprocal(scale(pochhammer, 2. * k + 1)); + inverse_factorial[k] = reciprocal(factorial); + log_pochhammer[k] = std::log(pochhammer.hi); + pochhammer = scale(pochhammer, k + .5); + kappa[k] = factorial.hi / pochhammer.hi; + factorial = scale(factorial, k + 1.); + } + } + static Pair scale(Pair a, double b) { + const Pair p = two_product(a.hi, b); + return fast_two_sum(p.hi, p.lo + a.lo * b); + } + static Pair reciprocal(Pair d) { + const double r = 1 / d.hi; + const Pair p = two_product(d.hi, r); + return fast_two_sum(r, r * (((1 - p.hi) - p.lo) - d.lo * r)); + } +}; + +// Taylor expansion about T = 0 for all orders at once, in double-word +// arithmetic: with x = R^2, +// G_m = F_m(T) - exp(-x) S_m, S_m = sum_j (-T)^j / j! p_{m+j}, +// p_k = P_k(x) / (2k+1) and P_k = 1F1(-k; 1/2; x). Q_k(R,0) = exp(-x) P_k, +// so |exp(-x) p_k| <= 1/(2k+1) (see direct_seed). The absolute terms of +// exp(-x) S_m sum to at most exp(T)/(2m+1), so binary64 summation would lose +// log10 exp(T) digits; the double-word error stays near u^2 exp(T). +// U_k = (1/2)_k P_k obeys the division-free recurrence +// U_{k+1} = c_k U_k + d_k U_{k-1}, c_k = 2k + 1/2 - x, d_k = -k(k - 1/2), +// and p_k = U_k / ((1/2)_k (2k+1)) uses a table. Returns false, leaving g +// unspecified, when the error estimate exceeds the budget or the request +// exceeds the tables. +bool origin_pair_taylor(double R, double T, int n, const Options& options, + Workspace& w, double* g) { + static const OriginPairTables tables; + constexpr int capacity = origin_pair_capacity; + // sum_{j>d} T^j / (j! (2j+1)) bounds the truncation of exp(-x) S_m. + const double target = options.precision / 128; + double term = 1, tail = infinity; + int degree = 0; + for (; n + degree + 1 < capacity; ++degree) { + term *= T / (degree + 1); + const double ratio = T / (degree + 2.); + if (ratio < 1) { + tail = term / ((1 - ratio) * (2. * degree + 3)); + if (tail <= target) break; + } + } + if (!(tail <= target)) return false; + const int last = n + degree; + const Pair x = two_product(R, R); + // |U_k| <= (1/2)_k exp(x) < 2^900 keeps every split exact, and the + // error estimate below finite. + if (!(x.hi <= 300) || tables.log_pochhammer[last] + x.hi > 600) return false; + const double ex = std::exp(-x.hi) * (1 - x.lo); + const double underflow = 16 * std::numeric_limits::denorm_min(); + + w.pair_terms.resize(14 * capacity); + double* const p_hi = w.pair_terms.data(); + double* const p_lo = p_hi + capacity; + double* const p_split_hi = p_lo + capacity; + double* const p_split_lo = p_split_hi + capacity; + double* const w_hi = p_split_lo + capacity; + double* const w_lo = w_hi + capacity; + double* const s_hi = w_lo + capacity; + double* const s_lo = s_hi + capacity; + double* const u_hi = s_lo + capacity; // U_k + double* const u_lo = u_hi + capacity; + double* const v = u_lo + capacity; // V_k + double* const alpha = v + capacity; + double* const P = alpha + capacity; // |exp(-x) P_k| + double* const Q = P + capacity; // |Q_k| + + // U_k, with c_k = c.hi + c_lo (the latter rounded once) and exact d_k. + // tau_{k+1} bounds the rounding in forming U_{k+1}, to first order: the + // products c.hi U_k.lo and c_lo U_k.hi and the sums after them (5 times), + // d U_{k-1}.lo (4 times), the error-free product and sum residuals (u + // times their leading parts, 3 and 1 times), the rounding of c_lo and the + // dropped c_lo U_k.lo. alpha_{k+1} = tau_{k+1} / k! is used below; it is + // independent of the recurrence and costs little here. + const double u = unit; + const Pair c0 = two_sum(.5, -x.hi); + u_hi[0] = 1; + u_lo[0] = 0; + u_hi[1] = c0.hi; + u_lo[1] = c0.lo - x.lo; + v[0] = 0; + v[1] = 1; + alpha[1] = u * std::abs(u_lo[1]) + underflow; + double magnitude_previous = 3 * u; // 5 |U_0.lo| + 3 u |U_0.hi| + for (int k = 1; k < last; ++k) { + const Pair c = two_sum(2. * k + .5, -x.hi); + const double c_lo = c.lo - x.lo, d = -k * (k - .5); + const Pair p1 = two_product(c.hi, u_hi[k]); + const Pair p2 = two_product(d, u_hi[k - 1]); + const double e1 = p1.lo + (c.hi * u_lo[k] + c_lo * u_hi[k]); + const double e2 = p2.lo + d * u_lo[k - 1]; + const Pair sum = two_sum(p1.hi, p2.hi); + u_hi[k + 1] = sum.hi; + u_lo[k + 1] = sum.lo + (e1 + e2); + v[k + 1] = mul_add(c.hi, v[k], d * v[k - 1]); + const double magnitude = 5 * std::abs(u_lo[k]) + 3 * u * std::abs(u_hi[k]); + const double tau = + u * (std::abs(c.hi) * magnitude - d * magnitude_previous + + u * std::abs(sum.hi)) + + std::abs(c_lo) * (std::abs(u_lo[k]) + 6 * u * std::abs(u_hi[k])) + + underflow; + alpha[k + 1] = tau * tables.inverse_factorial[k].hi; + magnitude_previous = magnitude; + } + // p_k = U_k c_k; P_k and Q_k = V_k / (1/2)_k for the estimate below. + p_hi[0] = 1; + p_lo[0] = 0; + P[0] = ex; + Q[0] = 0; + double max_p = 1, max_p_lo = 0; + for (int k = 1; k <= last; ++k) { + const Pair coefficient = tables.coefficient[k]; + const Pair product = two_product(u_hi[k], coefficient.hi); + p_hi[k] = product.hi; + p_lo[k] = + product.lo + (u_hi[k] * coefficient.lo + u_lo[k] * coefficient.hi); + const double order = 2. * k + 1; + P[k] = ex * order * std::abs(product.hi); + Q[k] = order * coefficient.hi * std::abs(v[k]); + max_p = std::max(max_p, std::abs(product.hi)); + max_p_lo = std::max(max_p_lo, std::abs(p_lo[k])); + } + // Rounding of p_k: with |c_k.lo| <= u c_k.hi the three roundings after + // the exact product are at most u (3 |p_k.lo| + 10 u |p_k.hi|), to first + // order; the table adds its relative error. + const double p_rounding = + ex * (u * (3 * max_p_lo + 10 * u * max_p) + tables.error * max_p) + + underflow; + + // An estimate of the error of exp(-x) P_k = (2k+1) exp(-x) p_k. The + // recurrence error follows its Green's function: a local error tau_i in + // U_i reaches U_k as + // tau_i (U_{i-1} V_k - V_{i-1} U_k) / C_{i-1}, + // with V a second solution (V_0 = 0, V_1 = 1) and the Casoratian + // C_i = U_i V_{i+1} - U_{i+1} V_i = i! (1/2)_i. With kappa_i = + // i! / (1/2)_{i+1} the error of exp(-x) P_k is at most + // sum_i alpha_i |g_ki|, where + // g_ki = exp(-x) (P_{i-1} Q_k - Q_{i-1} P_k) + // = exp(-x) Q_{i-1} Q_k sum_{l=i-1}^{k-1} kappa_l / (Q_l Q_{l+1}). + // Bounding the first form term by term suits oscillating U and V. Where + // both grow with the dominant solution (x >~ 4k) its terms cancel; the + // second form, as prefix sums of positive terms, then stays near the + // true error, and is evaluated when the first one misses the budget. + // This is a first-order estimate: V and the sums are rounded binary64, + // and the result enters doubled below. + const auto recurrence_error = [&](bool nested_form) -> double { + double sum_P = 0, sum_Q = 0, nested = 0, result = 0; + const double first = ex * alpha[1]; // the i = 1 term; Q_0 = 0 + for (int k = 1; k <= last; ++k) { + sum_P += P[k - 1] * alpha[k]; + sum_Q += Q[k - 1] * alpha[k]; + double bound = Q[k] * sum_P + P[k] * sum_Q; + if (nested_form) { + if (k > 1) nested += tables.kappa[k - 1] / (Q[k - 1] * Q[k]) * sum_Q; + bound = std::min(bound, Q[k] * (first + ex * nested)); + } + result = std::max(result, bound); + } + return result + p_rounding; + }; + + // (-T)^j / j!, with absolute error bounds. + Pair power{1, 0}; + double power_error = 0, weight_sum = 1, weight_lo_sum = 0, weight_error = 0; + w_hi[0] = 1; + w_lo[0] = 0; + for (int j = 1; j <= degree; ++j) { + const Pair p = two_product(power.hi, -T); + const double t = power.lo * -T, lo = p.lo + t; + power_error = T * power_error + u * (std::abs(t) + std::abs(lo)); + power = {p.hi, lo}; + const Pair f = tables.inverse_factorial[j]; + const Pair q = two_product(power.hi, f.hi); + const double ta = power.hi * f.lo, tb = power.lo * f.hi, tc = ta + tb; + w_hi[j] = q.hi; + w_lo[j] = q.lo + tc; + weight_error += + f.hi * power_error + + u * (std::abs(ta) + std::abs(tb) + std::abs(tc) + std::abs(w_lo[j])) + + tables.error * std::abs(q.hi) + underflow; + weight_sum += std::abs(q.hi); + weight_lo_sum += std::abs(w_lo[j]); + } + + // Rounding in the accumulation below, per order: |s_hi| <= A and + // |s_lo| <= AL; each step rounds the two lower-order products, their sum, + // e, the sum of e and the TwoSum error, and s_lo. + const double products = weight_sum * max_p_lo + weight_lo_sum * max_p; + const double exact_lo = u * weight_sum * max_p; + const double A = 1.01 * (weight_sum + weight_lo_sum) * (max_p + max_p_lo); + const double AL = 1.01 * (degree * u * A + products + exact_lo); + const double accumulation = + 1.02 * u * (6 * products + 3 * exact_lo + degree * (u * A + AL)); + const double fixed_error = + 2 * ex * (accumulation + weight_error * (max_p + max_p_lo)) + tail; + const double budget = options.precision / 4; + const double total = weight_sum + weight_lo_sum; + if (!(fixed_error + 2 * total * recurrence_error(false) <= budget) && + !(fixed_error + 2 * total * recurrence_error(true) <= budget)) + return false; + + // S_m = sum_j w_j p_{m+j}; independent m lanes vectorize. + std::copy(p_hi, p_hi + n + 1, s_hi); + std::copy(p_lo, p_lo + n + 1, s_lo); + if (fast_fma) { + for (int j = 1; j <= degree; ++j) { + const double wh = w_hi[j], wl = w_lo[j]; + const double *ph = p_hi + j, *pl = p_lo + j; + for (int m = 0; m <= n; ++m) { + const double p = wh * ph[m]; + const double e = std::fma(wh, ph[m], -p) + (wh * pl[m] + wl * ph[m]); + const double sum = s_hi[m] + p, t = sum - s_hi[m]; + s_lo[m] += ((s_hi[m] - (sum - t)) + (p - t)) + e; + s_hi[m] = sum; + } + } + } else { + // Split each operand once, outside the m loop. + for (int k = 0; k <= last; ++k) { + const Pair parts = split(p_hi[k]); + p_split_hi[k] = parts.hi; + p_split_lo[k] = parts.lo; + } + for (int j = 1; j <= degree; ++j) { + const double wh = w_hi[j], wl = w_lo[j]; + const Pair weight = split(wh); + const double *ph = p_hi + j, *pl = p_lo + j; + const double *qh = p_split_hi + j, *ql = p_split_lo + j; + for (int m = 0; m <= n; ++m) { + const double p = wh * ph[m]; + const double e = product_error(weight, {qh[m], ql[m]}, p) + + (wh * pl[m] + wl * ph[m]); + const double sum = s_hi[m] + p, t = sum - s_hi[m]; + s_lo[m] += ((s_hi[m] - (sum - t)) + (p - t)) + e; + s_hi[m] = sum; + } + } + } + boys_limit(T, n, w, g); + for (int m = 0; m <= n; ++m) g[m] = mul_add(-ex, s_hi[m] + s_lo[m], g[m]); + return true; +} + +void propagate_work(double T, int n, Workspace& w, double* g) { + g[0] = w.g0; + if (T == 0) { + for (int m = 1; m <= n; ++m) g[m] = (1 - w.q[m]) / (2. * m + 1); + return; + } + const double e = w.upward_last > 0 || w.seed_order > w.upward_last + 1 + ? exponential(T, w) + : 0; + for (int m = 0; m < w.upward_last; ++m) + g[m + 1] = mul_add(m + .5, g[m], .5 * (w.q[m] - e)) / T; + double next = w.seed; + if (w.seed_order == n) g[n] = next; + for (int m = w.seed_order - 1; m > w.upward_last; --m) { + next = mul_add(2 * T, next, e - w.q[m]) / (2. * m + 1); + if (m <= n) g[m] = next; + } +} + +void evaluate_work(double R, double T, int n, const Options& options, + Workspace& w, double* g) { + if (R == 0) { + evaluate_limit(T, n, CutoffLimit::zero, w, g); + return; + } + if (std::isinf(R)) { + evaluate_limit(T, n, CutoffLimit::boys, w, g); + return; + } + const auto limit = finite_cutoff_limit(R, T, n, options, w); + if (limit != CutoffLimit::none) { + evaluate_limit(T, n, limit, w, g); + return; + } + // Below the join at m ~ T - 1/2 the downward branch needs Q_m, which for + // T >~ 1 only Miller's boundary solve supplies, at several microseconds. + // The double-word origin expansion gives every order directly; its degree + // grows with T, which sets the upper end. Small R keeps the radius series. + if (n > 0 && T > 0 && T < n + .5 && T <= origin_pair_max_t && + !(R * R * ((n + 1.5 + T) / 3) <= .25) && + origin_pair_taylor(R, T, n, options, w, g)) + return; + w.reset_request(); + w.g0 = evaluate_g0(R, T, options, w).value; + if (n == 0) { + g[0] = w.g0; + return; + } + w.upward_last = 0; + w.seed_order = -1; + if (T > 0) + while (w.upward_last < n && w.upward_last + .5 < T) ++w.upward_last; + int last = T == 0 ? n : std::max(0, w.upward_last - 1); + if (T > 0 && w.upward_last < n) { + const auto seed = direct_seed(R, T, static_cast(n), options, w); + if (std::isfinite(seed.error)) { + w.seed = seed.value; + w.seed_order = n; + // The join Q_u is unused. If no downward step remains, only the + // upward inputs Q_0,...,Q_{u-1} are needed (possibly none). + last = n > w.upward_last + 1 ? n - 1 : w.upward_last - 1; + } else { + w.seed = 0; + w.seed_order = + zero_seed(T, n, options.precision / 32, options.max_extra_orders); + last = w.seed_order - 1; + } + } + const double tolerance = T != 0 && last >= 0 + ? q_tolerances(T, last, options.precision / 2, w) + : infinity; + if (last >= 0) evaluate_q(R, T, last, tolerance, options, w); + propagate_work(T, n, w, g); +} +} // namespace + +void evaluate(double* Gm, double R, double T, int mmax, double precision, + Scratch& scratch) { + const Options options{precision}; + validate(R, T, mmax, options); + if (!scratch.workspace_) // moved-from + scratch.workspace_.reset(new Workspace(mmax)); + evaluate_work(R, T, mmax, options, *scratch.workspace_, Gm); +} +} // namespace truncated_coulomb +} // namespace detail +} // namespace libint2 diff --git a/tools/generate_truncated_coulomb_reference.py b/tools/generate_truncated_coulomb_reference.py new file mode 100644 index 000000000..103dc17e8 --- /dev/null +++ b/tools/generate_truncated_coulomb_reference.py @@ -0,0 +1,168 @@ +#!/usr/bin/env python3 +"""Generate independent high-precision truncated-Coulomb test values. + +Requires mpmath. Differentiates the defining erf expression at fixed R; +T=0 uses the terminating confluent hypergeometric expression. Production +recurrences, CP2K code, tables and patch boundaries are not used. +""" +import argparse +import math +import random +from pathlib import Path +import mpmath as mp + + +def reference(r, t, mmax): + r, t = mp.mpf(r), mp.mpf(t) + if not t: + x = r * r + return [(1 - mp.exp(-x) * mp.hyp1f1(-m, mp.mpf('0.5'), x)) / + (2 * m + 1) for m in range(mmax + 1)] + + # Taylor-mode automatic differentiation of the defining expression. + # Guard against factorial scaling and cancellation near the sqrt branch + # point; this uses no G_m or boundary-derivative recurrence. + guard = int(mp.log10(mp.factorial(mmax))) + int( + mmax * max(0, -mp.log10(t))) + 20 + with mp.workdps(mp.mp.dps + guard): + n = mmax + root = [mp.sqrt(t)] + inverse_root = [1 / root[0]] + for k in range(1, n + 1): + root.append(root[-1] * (mp.mpf('1.5') - k) / (k * t)) + inverse_root.append(inverse_root[-1] * (mp.mpf('0.5') - k) / (k * t)) + + def erfc_jet(shift): + u = root.copy() + u[0] += shift + v = [-mp.fsum(u[j] * u[k-j] for j in range(k+1)) + for k in range(n+1)] + exponential = [mp.exp(v[0])] + for k in range(1, n+1): + exponential.append(mp.fsum(j*v[j]*exponential[k-j] + for j in range(1, k+1)) / k) + out = [mp.erfc(u[0])] + for k in range(1, n+1): + out.append(-2 / mp.sqrt(mp.pi) / k * mp.fsum( + exponential[j] * (k-j) * u[k-j] for j in range(k))) + return out + + minus, middle, plus = erfc_jet(-r), erfc_jet(0), erfc_jet(r) + h = [a - 2*b + c for a,b,c in zip(minus, middle, plus)] + return [(-1)**m * mp.factorial(m) * mp.sqrt(mp.pi) / 4 * + mp.fsum(h[j] * inverse_root[m-j] for j in range(m+1)) + for m in range(n+1)] + + + +def boys(t, mmax): + t = mp.mpf(t) + if not t: + return [mp.mpf(1) / (2 * m + 1) for m in range(mmax + 1)] + return [mp.gammainc(m + mp.mpf('0.5'), 0, t) / (2 * t**(m + mp.mpf('0.5'))) + for m in range(mmax + 1)] + + +def vrr_geometries(count, seed=20261025): + """Physical primitive geometries across the cutoff transition. + + Exponents zeta (bra) and eta (ket) are log-uniform in [1e-2, 1e6], the + cutoff Rc log-uniform in [1, 100], and sqrt(T) - R uniform in [-16, 16], + with R = Rc sqrt(rho), rho = zeta eta / (zeta + eta). W - P is the + Obara-Saika displacement for |PQ| = sqrt(T / rho). The first entries are + the worst cases of a 20000-geometry survey. + """ + def geometry(zeta, eta, rc, d): + rho = zeta * eta / (zeta + eta) + r = rc * math.sqrt(rho) + s = r + d + if s < 0 or r > 1e4: + return None + return zeta, rho, eta / (zeta + eta) * s / math.sqrt(rho), r, s * s + + out = [geometry(*g) for g in [(0.0112, 383, 8.65, 0.983), + (0.145, 1.26e5, 2.16, -0.0582), + (9.92e3, 8.12e4, 2.34, -4.32)]] + rng = random.Random(seed) + while len(out) < count: + g = geometry(10**rng.uniform(-2, 6), 10**rng.uniform(-2, 6), + 10**rng.uniform(0, 2), rng.uniform(-16, 16)) + if g: + out.append(g) + return out + + +def write_vrr(path, dps, mmax, count): + out = ['// Generated by tools/generate_truncated_coulomb_reference.py --vrr-output', + f'// mpmath {mp.__version__}; decimal precision {dps}; mmax {mmax}.', + '// G_m(R, T) for physical geometries across the cutoff transition;', + '// zeta, rho and W-P define the (L0|00) Obara-Saika recurrence.', + '// clang-format off', + 'struct TruncatedCoulombVRRReference {', + ' double zeta, rho, wp, r, t;', f' double g[{mmax + 1}];', '};', + 'static const TruncatedCoulombVRRReference truncated_coulomb_vrr_reference[] = {'] + for zeta, rho, wp, r, t in vrr_geometries(count): + values = reference(r, t, mmax) + with mp.workdps(dps + 40): + refined = reference(r, t, mmax) + f = boys(t, mmax) + # Stable to 1e-40 relative to the Coulomb scale max(|G_m|, F_m). + for a, b, c in zip(values, refined, f): + assert abs(a - b) <= mp.mpf('1e-40') * max(abs(b), c) + out.append(' {' + ', '.join(format(v, '.17g') for v in (zeta, rho, wp, r, t)) + ', {') + out.extend(' ' + ', '.join(format(float(v), '.17g') for v in values[k:k+4]) + ',' + for k in range(0, len(values), 4)) + out.append(' }},') + out.append('};') + out.append('// clang-format on') + path.write_text('\n'.join(out) + '\n') + + +def main(): + parser = argparse.ArgumentParser(description=__doc__) + unit_dir = Path('export/tests/unit') + if not unit_dir.is_dir(): + unit_dir = Path('tests/unit') + parser.add_argument('--output', type=Path, + default=unit_dir / 'truncated-coulomb-reference.h') + parser.add_argument('--dps', type=int, default=100) + parser.add_argument('--mmax', type=int, default=64) + parser.add_argument('--vrr-output', type=Path, + help='write the (L0|00) geometry set instead, mmax 16') + parser.add_argument('--vrr-count', type=int, default=96) + args = parser.parse_args() + mp.mp.dps = args.dps + if args.vrr_output: + write_vrr(args.vrr_output, args.dps, 16, args.vrr_count) + return + points = [('0.00000001', '0'), ('0.1', '0'), ('1', '0'), ('3', '0'), + ('0.1', '0.1'), ('1', '1'), ('3', '1.000001'), + ('1', '5'), ('3', '15'), ('5', '24'), ('5', '24.000001'), + ('5', '24.001'), ('8', '64'), ('8.000001', '64'), + ('5', '32'), ('5', '32.000001'), ('6', '36.000001'), + ('10', '40'), ('20', '400'), ('100', '10000')] + out = ['// Generated by tools/generate_truncated_coulomb_reference.py', + f'// mpmath {mp.__version__}; decimal precision {args.dps}; mmax {args.mmax}.', + '// Independent mathematical reference; rounded once to double.', + '// clang-format off', + 'struct TruncatedCoulombReference {', ' double r, t;', + f' double g[{args.mmax + 1}];', '};', + 'static const TruncatedCoulombReference truncated_coulomb_reference[] = {'] + for r, t in points: + values = reference(float(r), float(t), args.mmax) + # Verify stability under an additional 40 decimal working digits. + with mp.workdps(args.dps + 40): + refined = reference(float(r), float(t), args.mmax) + assert max(abs(a-b) for a, b in zip(values, refined)) < mp.mpf('1e-80') + out.append(' {' + r + ', ' + t + ', {') + out.extend(' ' + ', '.join(format(float(v), '.17g') for v in values[k:k+4]) + ',' + for k in range(0, len(values), 4)) + out.append(' }},') + print(f'generated R={r}, T={t}', flush=True) + out.append('};') + out.append('// clang-format on') + args.output.write_text('\n'.join(out) + '\n') + + +if __name__ == '__main__': + main()