Merge pull request #1 from boostorg/develop

Automatic merge of math/develop

Author: ryanelandt

.drone.star

@@ -16,13 +16,15 @@ windowsglobalimage="cppalliance/dronevs2019"
 def main(ctx):
 
   things_to_test = [ "special_fun", "distribution_tests", "mp", "misc", "interpolators", "quadrature", "autodiff", "long-running-tests", "float128_tests" ]
+  gcc13_things_to_test = [ "special_fun", "distribution_tests", "mp", "misc", "interpolators", "quadrature", "autodiff", "long-running-tests", "float128_tests", "new_floats" ]
   sanitizer_test = [ "special_fun", "distribution_tests", "misc", "interpolators", "quadrature", "float128_tests" ]
   gnu_5_stds = [ "gnu++14", "c++14" ]
   gnu_6_stds = [ "gnu++14", "c++14", "gnu++17", "c++17" ]
   clang_6_stds = [ "c++14", "c++17" ]
   gnu_9_stds = [ "gnu++14", "c++14", "gnu++17", "c++17", "gnu++2a", "c++2a" ]
   clang_10_stds = [ "c++14", "c++17", "c++2a" ]
   gnu_non_native = [ "gnu++17" ]
+  gcc13_stds = [ "c++23" ]
 
   result = []
 
@@ -59,6 +61,9 @@ def main(ctx):
       result.append(linux_cxx("Ubuntu g++ ARM64" + cxx + " " + suite, "g++", packages="g++", buildtype="boost", image="cppalliance/droneubuntu2204:multiarch", arch="arm64", environment={'TOOLSET': 'gcc', 'COMPILER': 'g++', 'CXXSTD': cxx, 'TEST_SUITE': suite, }, globalenv=globalenv))
     for cxx in gnu_non_native:
       result.append(osx_cxx("M1 Clang " + cxx + " " + suite, "clang++", buildscript="drone", buildtype="boost", xcode_version="14.1", environment={'TOOLSET': 'clang', 'CXXSTD': cxx, 'TEST_SUITE': suite, 'DEFINE': 'BOOST_MATH_NO_REAL_CONCEPT_TESTS,BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS,BOOST_MATH_MULTI_ARCH_CI_RUN', }, globalenv=globalenv))
+  for suite in gcc13_things_to_test:
+    for cxx in gcc13_stds:
+      result.append(linux_cxx("Ubuntu g++-13 " + cxx + " " + suite, "g++-13", packages="g++-13", buildtype="boost", image="cppalliance/droneubuntu2304:1", environment={'TOOLSET': 'gcc', 'COMPILER': 'g++-13', 'CXXSTD': cxx, 'TEST_SUITE': suite, }, globalenv=globalenv))
 
   return result
 

.github/workflows/ci.yml

@@ -9,6 +9,7 @@ on:
     branches:
       - master
       - develop
+      - feature/**
   pull_request:
   release:
     types: [published, created, edited]
@@ -47,7 +48,7 @@ jobs:
         if: steps.retry1.outcome=='failure'
         run: sudo apt-add-repository -y "ppa:ubuntu-toolchain-r/test"
       - name: Install packages
-        run: sudo apt install g++-12 clang-14 libgmp-dev libmpfr-dev libfftw3-dev
+        run: sudo apt-get install -y g++-12 clang-14 libgmp-dev libmpfr-dev libfftw3-dev
       - name: Checkout main boost
         run: git clone -b develop --depth 1 https://github.com/boostorg/boost.git ../boost-root
       - name: Update tools/boostdep
@@ -105,7 +106,7 @@ jobs:
         if: steps.retry1.outcome=='failure'
         run: sudo apt-add-repository -y "ppa:ubuntu-toolchain-r/test"
       - name: Install packages
-        run: sudo apt install g++-9 g++-11 clang-9 clang-10 libgmp-dev libmpfr-dev libfftw3-dev
+        run: sudo apt-get install -y g++-9 g++-11 clang-9 clang-10 libgmp-dev libmpfr-dev libfftw3-dev
       - name: Checkout main boost
         run: git clone -b develop --depth 1 https://github.com/boostorg/boost.git ../boost-root
       - name: Update tools/boostdep
@@ -364,7 +365,7 @@ jobs:
         if: steps.retry1.outcome=='failure'
         run: sudo apt-add-repository -y "ppa:ubuntu-toolchain-r/test"
       - name: Install packages
-        run: sudo apt install g++-10 libgmp-dev libmpfr-dev libfftw3-dev
+        run: sudo apt-get install -y g++-10 libgmp-dev libmpfr-dev libfftw3-dev
       - name: Checkout main boost
         run: git clone -b develop --depth 1 https://github.com/boostorg/boost.git ../boost-root
       - name: Update tools/boostdep
@@ -404,7 +405,7 @@ jobs:
         if: steps.retry1.outcome=='failure'
         run: sudo apt-add-repository -y "ppa:ubuntu-toolchain-r/test"
       - name: Install packages
-        run: sudo apt install clang-10 libgmp-dev libmpfr-dev libfftw3-dev
+        run: sudo apt-get install -y clang-10 libgmp-dev libmpfr-dev libfftw3-dev libtbb-dev
       - name: Checkout main boost
         run: git clone -b develop --depth 1 https://github.com/boostorg/boost.git ../boost-root
       - name: Update tools/boostdep
@@ -448,7 +449,7 @@ jobs:
         if: steps.retry1.outcome=='failure'
         run: sudo apt-add-repository -y "ppa:ubuntu-toolchain-r/test"
       - name: Install packages
-        run: sudo apt install g++-10 libgmp-dev libmpfr-dev libfftw3-dev
+        run: sudo apt-get install -y g++-10 libgmp-dev libmpfr-dev libfftw3-dev
       - name: Checkout main boost
         run: git clone -b develop --depth 1 https://github.com/boostorg/boost.git ../boost-root
       - name: Update tools/boostdep

doc/graphs/fourier_transform_daubechies.png

@@ -127,6 +127,48 @@ The 2 vanishing moment scaling function.
 [$../graphs/daubechies_8_scaling.svg]
 The 8 vanishing moment scaling function.
 
+Boost.Math also provides numerical evaluation of the Fourier transform of these functions.
+This is useful in sparse recovery problems where the measurements are taken in the Fourier basis.
+The usage is exhibited below:
+
+    #include <boost/math/special_functions/fourier_transform_daubechies_scaling.hpp>
+    using boost::math::fourier_transform_daubechies_scaling;
+    // Evaluate the Fourier transform of the 4-vanishing moment Daubechies scaling function at ω=1.8:
+    std::complex<float> hat_phi = fourier_transform_daubechies_scaling<float, 4>(1.8f);
+
+The Fourier transform convention is unitary with the sign of the imaginary unit being given in Daubechies Ten Lectures.
+In particular, this means that `fourier_transform_daubechies_scaling<float, p>(0.0)` returns 1/sqrt(2π).
+
+The implementation computes an infinite product of trigonometric polynomials as can be found from recursive application of the identity 𝓕[φ](ω) = m(ω/2)𝓕[φ](ω/2).
+This is neither particularly fast nor accurate, but there appears to be no literature on this extremely useful topic, and hence the naive method must suffice.
+
+[$../graphs/fourier_transform_daubechies.png]
+
+A benchmark can be found in `reporting/performance/fourier_transform_daubechies_performance.cpp`; the results on a ~2021 M1 Macbook pro are presented below:
+
+
+Run on (10 X 24.1212 MHz CPU s)
+CPU Caches:
+  L1 Data 64 KiB (x10)
+  L1 Instruction 128 KiB (x10)
+  L2 Unified 4096 KiB (x5)
+Load Average: 1.33, 1.52, 1.62
+-----------------------------------------------------------
+Benchmark                                              Time 
+-----------------------------------------------------------
+FourierTransformDaubechiesScaling<double, 1>        70.3 ns
+FourierTransformDaubechiesScaling<double, 2>         330 ns
+FourierTransformDaubechiesScaling<double, 3>         335 ns
+FourierTransformDaubechiesScaling<double, 4>         364 ns
+FourierTransformDaubechiesScaling<double, 5>         386 ns
+FourierTransformDaubechiesScaling<double, 6>         436 ns
+FourierTransformDaubechiesScaling<double, 7>         447 ns
+FourierTransformDaubechiesScaling<double, 8>         473 ns
+FourierTransformDaubechiesScaling<double, 9>         503 ns
+FourierTransformDaubechiesScaling<double, 10>        554 ns
+
+Due to the low accuracy of this method, `float` precision is arg-promoted to `double`, and hence takes just as long as `double` precision to execute.
+
 [heading References]
 
 * Daubechies, Ingrid. ['Ten Lectures on Wavelets.] Vol. 61. Siam, 1992.

doc/sf/daubechies.qbk

@@ -0,0 +1,43 @@
+//  (C) Copyright Nick Thompson 2023.
+//  Use, modification and distribution are subject to the
+//  Boost Software License, Version 1.0. (See accompanying file
+//  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
+
+#include <utility>
+#include <boost/math/filters/daubechies.hpp>
+#include <boost/math/tools/polynomial.hpp>
+#include <boost/multiprecision/cpp_bin_float.hpp>
+#include <boost/math/constants/constants.hpp>
+
+using std::pow;
+using boost::multiprecision::cpp_bin_float_100;
+using boost::math::filters::daubechies_scaling_filter;
+using boost::math::tools::polynomial;
+using boost::math::constants::half;
+using boost::math::constants::root_two;
+
+template<typename Real, size_t N>
+std::vector<Real> get_constants() {
+   auto h = daubechies_scaling_filter<cpp_bin_float_100, N>();
+   auto p = polynomial<cpp_bin_float_100>(h.begin(), h.end());
+
+   auto q = polynomial({half<cpp_bin_float_100>(), half<cpp_bin_float_100>()});
+   q = pow(q, N);
+   auto l = p/q;
+   return l.data(); 
+}
+
+template<typename Real>
+void print_constants(std::vector<Real> const & l) {
+   std::cout << std::setprecision(std::numeric_limits<Real>::digits10 -10);
+   std::cout << "return std::array<Real, " << l.size() << ">{";
+   for (size_t i = 0; i < l.size() - 1; ++i) {
+       std::cout << "BOOST_MATH_BIG_CONSTANT(Real, std::numeric_limits<Real>::digits, " << l[i]/root_two<Real>() << "), ";
+   }
+   std::cout << "BOOST_MATH_BIG_CONSTANT(Real, std::numeric_limits<Real>::digits, " << l.back()/root_two<Real>() << ")};\n";
+}
+
+int main() {
+   auto constants = get_constants<cpp_bin_float_100, 1>();
+   print_constants(constants);
+}

example/calculate_fourier_transform_daubechies_constants.cpp

@@ -0,0 +1,60 @@
+// boost-no-inspect
+//  (C) Copyright Nick Thompson 2023.
+//  Use, modification and distribution are subject to the
+//  Boost Software License, Version 1.0. (See accompanying file
+//  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
+
+#include <boost/math/special_functions/fourier_transform_daubechies.hpp>
+#include <boost/math/tools/ulps_plot.hpp>
+
+using boost::math::fourier_transform_daubechies_scaling;
+using boost::math::tools::ulps_plot;
+
+template<int p>
+void real_part() {
+   auto phi_real_hi_acc = [](double omega) {
+       auto z = fourier_transform_daubechies_scaling<double, p>(omega);
+       return z.real();
+   };
+
+   auto phi_real_lo_acc = [](float omega) {
+      auto z = fourier_transform_daubechies_scaling<float, p>(omega);
+      return z.real();
+   };
+   auto plot = ulps_plot<decltype(phi_real_hi_acc), double, float>(phi_real_hi_acc, float(0.0), float(100.0), 20000);
+   plot.ulp_envelope(false);
+   plot.add_fn(phi_real_lo_acc);
+   plot.clip(100);
+   plot.title("Accuracy of 𝔑(𝓕[𝜙](ω)) with " + std::to_string(p) + " vanishing moments.");
+   plot.write("real_ft_daub_scaling_"  + std::to_string(p) + ".svg");
+ 
+}
+
+template<int p>
+void imaginary_part() {
+   auto phi_imag_hi_acc = [](double omega) {
+       auto z = fourier_transform_daubechies_scaling<double, p>(omega);
+       return z.imag();
+   };
+
+   auto phi_imag_lo_acc = [](float omega) {
+      auto z = fourier_transform_daubechies_scaling<float, p>(omega);
+      return z.imag();
+   };
+   auto plot = ulps_plot<decltype(phi_imag_hi_acc), double, float>(phi_imag_hi_acc, float(0.0), float(100.0), 20000);
+   plot.ulp_envelope(false);
+   plot.add_fn(phi_imag_lo_acc);
+   plot.clip(100);
+   plot.title("Accuracy of 𝕴(𝓕[𝜙](ω)) with " + std::to_string(p) + " vanishing moments.");
+   plot.write("imag_ft_daub_scaling_"  + std::to_string(p) + ".svg");
+ 
+}
+
+
+int main() {
+   real_part<3>();
+   imaginary_part<3>();
+   real_part<6>();
+   imaginary_part<6>();
+   return 0;
+}

example/fourier_transform_daubechies_ulp_plot.cpp

@@ -45,6 +45,10 @@
 #  include <cstdio>
 #endif
 
+#if __has_include(<stdfloat>)
+#  include <stdfloat>
+#endif
+
 namespace boost{ namespace math{
 
 namespace concepts
@@ -79,6 +83,12 @@ class real_concept
 #ifdef BOOST_MATH_USE_FLOAT128
    real_concept(BOOST_MATH_FLOAT128_TYPE c) : m_value(c){}
 #endif
+#ifdef __STDCPP_FLOAT32_T__
+   real_concept(std::float32_t c) : m_value(static_cast<real_concept_base_type>(c)){}
+#endif
+#ifdef __STDCPP_FLOAT64_T__
+   real_concept(std::float64_t c) : m_value(static_cast<real_concept_base_type>(c)){}
+#endif
 
    // Assignment:
    real_concept& operator=(char c) { m_value = c; return *this; }
@@ -96,6 +106,12 @@ class real_concept
    real_concept& operator=(float c) { m_value = c; return *this; }
    real_concept& operator=(double c) { m_value = c; return *this; }
    real_concept& operator=(long double c) { m_value = c; return *this; }
+   #ifdef __STDCPP_FLOAT32_T__
+   real_concept& operator=(std::float32_t c) { m_value = c; return *this; }
+   #endif
+   #ifdef __STDCPP_FLOAT64_T__
+   real_concept& operator=(std::float64_t c) { m_value = c; return *this; }
+   #endif
 
    // Access:
    real_concept_base_type value()const{ return m_value; }

include/boost/math/concepts/real_concept.hpp

@@ -24,7 +24,7 @@
 #pragma GCC system_header
 #endif
 
-  #if defined(BOOST_CSTDFLOAT_HAS_INTERNAL_FLOAT128_T) && defined(BOOST_MATH_USE_FLOAT128) && !defined(BOOST_CSTDFLOAT_NO_LIBQUADMATH_SUPPORT)
+  #if defined(BOOST_CSTDFLOAT_HAS_INTERNAL_FLOAT128_T) && defined(BOOST_MATH_USE_FLOAT128) && !defined(BOOST_CSTDFLOAT_NO_LIBQUADMATH_SUPPORT) && (!defined(__GNUC__) || (defined(__GNUC__) && __GNUC__ < 14))
 
     #include <limits>
     #include <boost/math/tools/nothrow.hpp>

include/boost/math/cstdfloat/cstdfloat_limits.hpp

@@ -1491,7 +1491,7 @@ fvar<RealType, Order> sqrt(fvar<RealType, Order> const& cr) {
   BOOST_IF_CONSTEXPR (order == 0)
     return fvar<RealType, Order>(*derivatives);
   else {
-    root_type numerator = 0.5;
+    root_type numerator = root_type(0.5);
     root_type powers = 1;
 #ifndef BOOST_NO_CXX17_IF_CONSTEXPR
     derivatives[1] = numerator / *derivatives;

include/boost/math/differentiation/autodiff.hpp

@@ -153,7 +153,7 @@ namespace detail {
 
       const Real eps = (numeric_limits<Real>::epsilon)();
       // Error bound ~eps^4/5
-      Real h = pow(11.25*eps, static_cast<Real>(1) / static_cast<Real>(5));
+      Real h = pow(Real(11.25)*eps, static_cast<Real>(1) / static_cast<Real>(5));
       h = detail::make_xph_representable(x, h);
       Real ymth = f(x - 2 * h);
       Real yth = f(x + 2 * h);
@@ -222,7 +222,7 @@ namespace detail {
       // Mathematica code to get the error:
       // Series[(f[x+h]-f[x-h])*(4/5) + (1/5)*(f[x-2*h] - f[x+2*h]) + (4/105)*(f[x+3*h] - f[x-3*h]) + (1/280)*(f[x-4*h] - f[x+4*h]), {h, 0, 9}]
       // If we used Kahan summation, we could get the max error down to h^8|f^(9)(x)|/630 + |f(x)|eps/h.
-      Real h = pow(551.25*eps, static_cast<Real>(1) / static_cast<Real>(9));
+      Real h = pow(Real(551.25)*eps, static_cast<Real>(1) / static_cast<Real>(9));
       h = detail::make_xph_representable(x, h);
 
       Real yh = f(x + h);