Locate Bessel function from libboost

If option `use_boost` is true, search the system for `Boost.Math` via
pkg-config. If not exist, download the Boost.Math header-only
library. Configure it in standalone mode.

Otherwise, if `use_boost=false`, check for `std::cyl_bessel_j` in the
C++17 standard.

.github/workflows/validate-build.yml

@@ -5,7 +5,7 @@ name: Compile and test C++
 
 on:
   push:
-    branches: [ master ]
+    branches: [ main ]
   pull_request:
     types: [edited, submitted]
 
@@ -17,16 +17,16 @@ jobs:
       matrix:
         #os: [ubuntu-20.04, ubuntu-latest, windows-latest, macos-latest]
         os: [ubuntu-20.04, ubuntu-latest]
-        python-version: ["3.10"]
+        use_boost: ["false", "true"]
 
     steps:
     - name: Fetch sources
       uses: actions/checkout@v4
 
-    - name: Set up Python ${{ matrix.python-version }}
+    - name: Set up Python 3.10
       uses: actions/setup-python@v5
       with:
-        python-version: ${{ matrix.python-version }}
+        python-version: "3.10"
         cache: pip
 
     - name: Install gcc toolchain
@@ -48,7 +48,10 @@ jobs:
         path: 3rdparty/packagecache/
 
     - name: Resolve C++ build dependencies (non-Windows)
-      run: meson setup build/
+      run: >-
+        meson setup
+        -Duse_boost=${{ matrix.use_boost }}
+        build/
       if: startsWith(matrix.os, 'window') == false
 
     #- name: Resolve C++ build dependencies (msvc toolchain)

3rdparty/armadillo-code.wrap

@@ -4,7 +4,7 @@ revision = 14.0.2
 depth = 1
 
 patch_directory = armadillo-code
-diff_files = armadillo-code/0001-Patch-armdillo-to-support-besselj.patch
+diff_files = 0001-Patch-armadillo-to-support-besselj.patch
 
 [provide]
-armadillo-code = armadillo_dep
+armadillo-code = armadillo_dep

3rdparty/boost-math.wrap

@@ -0,0 +1,10 @@
+[wrap-file]
+directory = math-boost-1.83.0
+source_url = https://github.com/boostorg/math/archive/boost-1.83.0.tar.gz
+source_filename = boost-1.83.0.tar.gz
+source_hash = 53e5f7539a66899fe0fca3080405cbd5f7959da5394ec13664746741aece1705
+
+patch_directory = boost-math
+
+[provide]
+boost_math = boost_math_dep

3rdparty/packagefiles/0001-Patch-armadillo-to-support-besselj.patch

@@ -0,0 +1,50 @@
+From de3cda8744e89b65b186b721853a611eb8d505cd Mon Sep 17 00:00:00 2001
+From: Antony Chan <[email protected]>
+Date: Tue, 26 Nov 2024 17:23:16 -0800
+Subject: [PATCH] Patch armadillo to support besselj
+
+---
+ include/armadillo_bits/eop_aux.hpp | 20 ++++++++++++++++++--
+ 1 file changed, 18 insertions(+), 2 deletions(-)
+
+diff --git a/include/armadillo_bits/eop_aux.hpp b/include/armadillo_bits/eop_aux.hpp
+index b63eba6..8da912b 100644
+--- a/include/armadillo_bits/eop_aux.hpp
++++ b/include/armadillo_bits/eop_aux.hpp
+@@ -20,6 +20,9 @@
+ //! @{
+
+
++#ifdef BOOST_HAS_BESSEL
++#include <boost/math/special_functions/bessel.hpp>
++#endif
+
+ //! use of the SFINAE approach to work around compiler limitations
+ //! http://en.wikipedia.org/wiki/SFINAE
+@@ -138,8 +141,21 @@ class eop_aux
+
+   template<typename T1, typename T2> arma_inline static typename   arma_integral_only<T1>::result pow (const T1 base, const T2 exponent) { return T1( std::pow( double(base), double(exponent) ) ); }
+   template<typename T1, typename T2> arma_inline static typename arma_real_or_cx_only<T1>::result pow (const T1 base, const T2 exponent) { return T1( std::pow(        base,         exponent  ) ); }
+-  
+-  
++
++// Wrap the scalar Bessel function of the first kind. Enable 32-bit and 64-bit
++// algorithm.
++template <uint8_t order, typename eT>
++arma_inline static typename arma_real_only<eT>::result
++besselj(const eT x) {
++#ifdef STD_HAS_BESSEL
++    using std::cyl_bessel_j;
++#elif defined(BOOST_HAS_BESSEL)
++    using boost::math::cyl_bessel_j;
++#endif
++    return cyl_bessel_j(double(order), x);
++}
++
++
+   template<typename eT>
+   arma_inline
+   static
+-- 
+2.25.1
+

3rdparty/packagefiles/boost-math/meson.build

@@ -0,0 +1,10 @@
+project('boost-math', 'cpp',
+    version: '1.83.0',
+)
+
+boost_math_inc = include_directories('include')
+
+boost_math_dep = declare_dependency(
+    include_directories: boost_math_inc,
+    compile_args: ['-DBOOST_MATH_STANDALONE'],
+)

besselj_armadillo_support/inc/armadillo_besselj_support.hpp

@@ -6,6 +6,10 @@
 #include <cmath>
 #include <cstdint>
 
+#ifdef BOOST_HAS_BESSEL
+#include <boost/math/special_functions/bessel.hpp>
+#endif
+
 #include <armadillo>
 
 namespace arma {

besselj_armadillo_support/meson.build

@@ -1,33 +1,69 @@
 armadillo_besselj_support_inc = include_directories('inc')
-armadillo_dep = dependency('armadillo-code')
+armadillo_dep = subproject('armadillo-code').get_variable('armadillo_dep')
 
 # TODO(Antony): use "cxx.has_function" when it supports std namespace.
-if not cxx.compiles(
-    '''#include <cmath>
+if get_option('use_boost')
+    warning('Fallback to Boost.Math for the Bessel function implementation.')
 
+    bessel_source = 'BOOST_HAS_BESSEL'
+    boost_math_dep = dependency('boost', modules: ['math'], fallback: [
+        'boost-math', 'boost_math_dep',
+    ])
+elif cxx.compiles(
+    '''#include <cmath>
 int main() {
     return int(std::cyl_bessel_j(0.0, 0.0));
 }
 ''')
+    bessel_source = 'STD_HAS_BESSEL'
+    boost_math_dep = []
+else
     error('Bessel function (of the first kind) does not exist in the C++ toolchain. Is this ISO C++17 compilant?')
 endif
 
 test_besselj_support_exe = executable('test_besselj_support',
     sources: [
         'tests/test-besselj.cpp',
+        #'tests/test-armadillo-support.cpp',
+    ],
+    cpp_args: [
+        '-D' + bessel_source,
+    ],
+    dependencies: [
+        catch2_dep,
+        boost_math_dep,
+    ],
+)
+
+test_besselj_wrapper_exe = executable('test_besselj_armadillo_wrapper',
+    sources: [
         'tests/test-armadillo-support.cpp',
     ],
-    include_directories: armadillo_besselj_support_inc,
+    cpp_args: [
+        '-D' + bessel_source,
+    ],
+    include_directories: [
+        armadillo_besselj_support_inc,
+    ],
     dependencies: [
         catch2_dep,
         armadillo_dep,
+        boost_math_dep,
     ],
 )
 
-test('Besselj wrapper',
+test('Besselj native support',
     test_besselj_support_exe,
     args: [
         '-r', 'tap',
     ],
     protocol: 'tap',
+)
+
+test('Besselj Armadillo wrapper',
+    test_besselj_wrapper_exe,
+    args: [
+        '-r', 'tap',
+    ],
+    protocol: 'tap',
 )

besselj_armadillo_support/tests/test-armadillo-support.cpp

@@ -1,5 +1,7 @@
-#include <armadillo>
 #include <armadillo_besselj_support.hpp>
+
+//
+#include <armadillo>
 #include <catch2/catch_test_macros.hpp>
 #include <cstdint>
 

besselj_armadillo_support/tests/test-besselj.cpp

@@ -3,15 +3,19 @@
 #include <catch2/catch_test_macros.hpp>
 #include <catch2/matchers/catch_matchers_floating_point.hpp>
 
+#ifdef STD_HAS_BESSEL
+using std::cyl_bessel_j;
+#elif defined(BOOST_HAS_BESSEL)
+#include <boost/math/special_functions/bessel.hpp>
+using boost::math::cyl_bessel_j;
+#endif
 TEST_CASE("Toolchain has std::cyl_bessel_j") {
-    using std::cyl_bessel_j;
     using std::isnan;
 
     REQUIRE(!isnan(cyl_bessel_j(0.0, 0.0)));
 }
 
 TEST_CASE("Besselj spot check") {
-    using std::cyl_bessel_j;
     using Catch::Matchers::WithinAbs;
     using Catch::Matchers::WithinRel;
 

examples/generate-psf.cpp

@@ -23,7 +23,7 @@ int main() {
 
     const auto psf =
         makePSF(params, {0.1_um, 0.25_um}, {256, 128}, 0.610_um, precision);
-#ifdef  ARMA_USE_HDF5
+#ifdef ARMA_USE_HDF5
     std::cout << "Saving volume to HDF5...\n";
     psf.save(hdf5_name("psf.h5", "psf", hdf5_opts::trans));
     std::cout << R"(Done.
@@ -52,4 +52,4 @@ Import to Python:
         std::cout << "PSF XZ plane saved to 'psf_xz.pgm'.\n";
     }
     return 0;
-}
+}

meson.options

@@ -1 +1,2 @@
-option('install_examples', type: 'boolean', value: false)
+option('install_examples', type: 'boolean', value: false)
+option('use_boost', type: 'boolean', value: false)