From 612c5cbf31c0e816030fdb959ca2005ac1f50dfc Mon Sep 17 00:00:00 2001
From: Matt Borland <matt@mattborland.com>
Date: Mon, 4 Dec 2023 12:42:55 +0100
Subject: [PATCH] Add erfc for z < 3.5

---
 include/boost/decimal/detail/cmath/erf.hpp | 33 ++++++++++++++++++++++
 test/test_constants.cpp                    | 23 +++++++++++++++
 2 files changed, 56 insertions(+)

diff --git a/include/boost/decimal/detail/cmath/erf.hpp b/include/boost/decimal/detail/cmath/erf.hpp
index 7475c160c..fb2e1c786 100644
--- a/include/boost/decimal/detail/cmath/erf.hpp
+++ b/include/boost/decimal/detail/cmath/erf.hpp
@@ -544,6 +544,39 @@ constexpr auto erf_impl<decimal128>(decimal128 z, bool invert) noexcept -> decim
             constexpr decimal128 offset {UINT64_C(225), -2};
             result = Y + tools::evaluate_polynomial(P, z - offset) / tools::evaluate_polynomial(Q, z - offset);
         }
+        else if (z < decimal128{UINT64_C(35), -1})
+        {
+            // Maximum Deviation Found:                     8.126e-37
+            // Expected Error Term:                         -8.126e-37
+            // Maximum Relative Change in Control Points:   1.363e-04
+            // Max Error found at long double precision =   1.747062e-36
+            constexpr decimal128 Y {detail::uint128{UINT64_C(292937225141646), UINT64_C(6920050031251800064)}, -34};
+            constexpr std::array<decimal128, 9> P = {
+                decimal128{detail::uint128{UINT64_C(182706965924257), UINT64_C(1687510779571187718)}, -36, true},
+                decimal128{detail::uint128{UINT64_C(56892448168985), UINT64_C(572440462241151398)}, -35},
+                decimal128{detail::uint128{UINT64_C(80518338580783), UINT64_C(5160816315849708842)}, -35},
+                decimal128{detail::uint128{UINT64_C(442730178280838), UINT64_C(9281603077550627672)}, -36},
+                decimal128{detail::uint128{UINT64_C(135371629264938), UINT64_C(7268401433168016132)}, -36},
+                decimal128{detail::uint128{UINT64_C(252380094364866), UINT64_C(3735236004636993191)}, -37},
+                decimal128{detail::uint128{UINT64_C(287925910284089), UINT64_C(6157066008997322426)}, -38},
+                decimal128{detail::uint128{UINT64_C(186226232526489), UINT64_C(2794677292908361186)}, -39},
+                decimal128{detail::uint128{UINT64_C(526445427809093), UINT64_C(11759659595468142822)}, -41}
+            };
+            constexpr std::array<decimal128, 9> Q = {
+                decimal128{1},
+                decimal128{detail::uint128{UINT64_C(86688065670866), UINT64_C(11737797169918939734)}, -33},
+                decimal128{detail::uint128{UINT64_C(61583053693636), UINT64_C(8177778869190231158)}, -33},
+                decimal128{detail::uint128{UINT64_C(254010066013673), UINT64_C(14314255351052138662)}, -34},
+                decimal128{detail::uint128{UINT64_C(66581844722135), UINT64_C(10035464857808786462)}, -34},
+                decimal128{detail::uint128{UINT64_C(113662830747969), UINT64_C(10480615872240633506)}, -35},
+                decimal128{detail::uint128{UINT64_C(123515411355391), UINT64_C(5270626324694473614)}, -36},
+                decimal128{detail::uint128{UINT64_C(78194463948513), UINT64_C(4344969105995523842)}, -37},
+                decimal128{detail::uint128{UINT64_C(221048990718863), UINT64_C(13286283565256558792)}, -39}
+            };
+
+            constexpr decimal128 offset {UINT64_C(3), -1};
+            result = Y + tools::evaluate_polynomial(P, z - offset) / tools::evaluate_polynomial(Q, z - offset);
+        }
 
         decimal128 hi {};
         decimal128 lo {};
diff --git a/test/test_constants.cpp b/test/test_constants.cpp
index e1a02b9c8..f60938de6 100644
--- a/test/test_constants.cpp
+++ b/test/test_constants.cpp
@@ -236,6 +236,29 @@ int main()
     print_value("0.686843205749767250666787987163701209e-4"_DL, "Q8");
     print_value("0.192093541425429248675532015101904262e-5"_DL, "Q9");
 
+    std::cerr << "---------- z < 3.5 --------\n";
+    print_value("0.54037380218505859375"_DL, "Y");
+
+    print_value("-0.0033703486408887424921155540591370375"_DL, "P0");
+    print_value("0.0104948043110005245215286678898115811"_DL, "P1");
+    print_value("0.0148530118504000311502310457390417795"_DL, "P2");
+    print_value("0.00816693029245443090102738825536188916"_DL, "P3");
+    print_value("0.00249716579989140882491939681805594585"_DL, "P4");
+    print_value("0.0004655591010047353023978045800916647"_DL, "P5");
+    print_value("0.531129557920045295895085236636025323e-4"_DL, "P6");
+    print_value("0.343526765122727069515775194111741049e-5"_DL, "P7");
+    print_value("0.971120407556888763695313774578711839e-7"_DL, "P8");
+
+    print_value("1"_DL, "Q0");
+    print_value("1.59911256167540354915906501335919317"_DL, "Q1");
+    print_value("1.136006830764025173864831382946934"_DL, "Q2");
+    print_value("0.468565867990030871678574840738423023"_DL, "Q3");
+    print_value("0.122821824954470343413956476900662236"_DL, "Q4");
+    print_value("0.0209670914950115943338996513330141633"_DL, "Q5");
+    print_value("0.00227845718243186165620199012883547257"_DL, "Q6");
+    print_value("0.000144243326443913171313947613547085553"_DL, "Q7");
+    print_value("0.407763415954267700941230249989140046e-5"_DL, "Q8");
+
     return 1;
     #endif