Modernize legendre_convert Boost usage

ALPS3_MIGRATION_PLAN.md

@@ -387,19 +387,24 @@ unchanged. A text-input regression covers trailing blank lines.
 
 **2.3B: `legendre_convert`** (decided 2026-09-24), in this order:
 
-1. Add regression cases for `legendre_convert` first (it has no tests), with
-   references from the current Boost build, like the `kk` cases in 2.2.
-2. Replace Boost.Random by `<random>`. The `mt19937` engine gives identical
+1. **Done 2026-09-24:** add a deterministic regression case for
+   `legendre_convert`, with a reference from the current Boost build. It covers
+   the transform, tail enforcement, back-continuation, and Matsubara
+   convergence.
+2. **Done 2026-09-24:** replace Boost.Random by `<random>`. The `mt19937` engine gives identical
    numbers; the normal variates differ (different algorithm), which is
    harmless because the error estimate is seeded from the clock.
-3. Replace `boost::math::factorial` by a product (only small arguments occur).
-4. Replace `boost::math::legendre_p` by the standard three-term recurrence.
+3. **Done 2026-09-24:** replace `boost::math::factorial` by a product (only
+   small arguments occur).
+4. **Done 2026-09-24:** replace `boost::math::legendre_p` by the standard
+   three-term recurrence.
    `std::legendre` is not an option: libc++ (AppleClang) does not implement the
    C++17 special math functions.
 
 Keep `boost::math::sph_bessel` (also missing in libc++; our own version would
 need careful checking at large l and argument) and `program_options` (no
-standard equivalent).
+standard equivalent). Reconfirmed 2026-09-24: these are intentional retained
+Boost boundaries for now.
 
 **2.3C: command-line behavior** (updates the CLI references on purpose)
 

legendre_convert/gaussian_noise.hpp

@@ -0,0 +1,36 @@
+/*****************************************************************************
+*
+* ALPS Project Applications
+*
+* Copyright (C) 1998-2016 ALPS Collaboration
+*
+* ALPS Project: https://alps.comp-phys.org/
+* SPDX-License-Identifier: MIT
+*
+*****************************************************************************/
+
+#pragma once
+
+#include <random>
+#include <stdexcept>
+#include <vector>
+
+inline std::vector<double> generateGaussNoise(const std::vector<double>& data,
+                                              const std::vector<double>& error,
+                                              std::mt19937& rng) {
+    if (data.size() != error.size())
+        throw std::invalid_argument("data and error vectors must have equal size");
+
+    std::vector<double> noisy_data(data.size());
+    for (std::size_t i = 0; i < data.size(); ++i) {
+        if (error[i] < 0.0)
+            throw std::invalid_argument("Gaussian standard deviation must be nonnegative");
+        if (error[i] == 0.0)
+            noisy_data[i] = data[i];
+        else {
+            std::normal_distribution<> distribution(data[i], error[i]);
+            noisy_data[i] = distribution(rng);
+        }
+    }
+    return noisy_data;
+}

legendre_convert/legendre_convert.cpp

@@ -9,22 +9,40 @@
 *
 *****************************************************************************/
 
-#include <iostream>
+#include <cmath>
+#include <complex>
+#include <ctime>
 #include <fstream>
+#include <iomanip>
+#include <iostream>
+#include <random>
+#include <stdexcept>
+#include <string>
+#include <utility>
+#include <vector>
+
+#include "gaussian_noise.hpp"
 #include <boost/program_options.hpp>
-#include <boost/random/mersenne_twister.hpp>
-#include <boost/random/normal_distribution.hpp>
-#include <boost/random/variate_generator.hpp>
-#include <boost/random.hpp>
-#include <boost/math/special_functions/factorials.hpp>
-#include <boost/math/special_functions/legendre.hpp> //needed for Legendre transform
 #include <boost/math/special_functions/bessel.hpp> 
 
 typedef std::vector<double> vector_type;
 typedef std::pair<double,double> return_type;
 typedef std::complex<double> Complex ;
 namespace bmth = boost::math;
 
+double legendre_polynomial(int l, double x){
+    if(l == 0)
+        return 1.0;
+    double p_lm2 = 1.0;
+    double p_lm1 = x;
+    for(int n=2;n<=l;n++){
+        const double p_l = ((2*n-1)*x*p_lm1-(n-1)*p_lm2)/n;
+        p_lm2 = p_lm1;
+        p_lm1 = p_l;
+    }
+    return p_lm1;
+}
+
 
 ///E8 of Boehnke, et al
 ///high frequency expansion of Gl (see Tl)
@@ -38,7 +56,10 @@ double tl(int l, int p){
         double qsum =1;
         for(int q=-p+2;q<p;q++)
             qsum*=l+q;
-        return std::pow(-1.0,p)*2*std::sqrt(2*l+1)*qsum/bmth::factorial<double>(p-1);
+        double factorial = 1.0;
+        for(int factor=2;factor<p;factor++)
+            factorial*=factor;
+        return std::pow(-1.0,p)*2*std::sqrt(2*l+1)*qsum/factorial;
     }
 }
 ///G(i\omega)=\sumT_{n\ell}G_\ell
@@ -83,36 +104,36 @@ double generateGl(const vector_type &gtau,const vector_type &tau_points, int l,d
 
     	--ndat_;
     	double dtau = tau_points[1]-tau_points[0];
-    	gsum+=gtau[0]*bmth::legendre_p(l,2*tau_points[0]/beta-1.0);
+        gsum+=gtau[0]*legendre_polynomial(l,2*tau_points[0]/beta-1.0);
 	    for(int i=1;i<ndat_/2;i++)
-	    	gsum+=2*gtau[2*i]*bmth::legendre_p(l,2*tau_points[2*i]/beta-1);
+        gsum+=2*gtau[2*i]*legendre_polynomial(l,2*tau_points[2*i]/beta-1);
 	    for(int i=1;i<ndat_/2+1;i++)
-	    	gsum+=4*gtau[2*i-1]*bmth::legendre_p(l,2*tau_points[2*i-1]/beta-1);
+        gsum+=4*gtau[2*i-1]*legendre_polynomial(l,2*tau_points[2*i-1]/beta-1);
 
-	    gsum+=gtau[ndat_]*bmth::legendre_p(l,2*tau_points[ndat_]/beta-1.0);
+	    gsum+=gtau[ndat_]*legendre_polynomial(l,2*tau_points[ndat_]/beta-1.0);
 	    gsum*= dtau/3;
     }
     else{
 
 		//Riemann sums
 	    /*for(int i=0;i<ndat_-1;i++){
 	        double dtau = tau_points[i+1]-tau_points[i];
-	        gsum+=dtau*(bmth::legendre_p(l, 2*tau_points[i]/beta-1)*gtau[i]);
+	        gsum+=dtau*(legendre_polynomial(l, 2*tau_points[i]/beta-1)*gtau[i]);
 	    }
 	    //endpoint:
 	    double dtau=tau_points[ndat_-1]-tau_points[ndat_-2];
-	    gsum+=dtau*(bmth::legendre_p(l, 2*tau_points[ndat_-1]/beta-1)*gtau[ndat_-1]);*/
+	    gsum+=dtau*(legendre_polynomial(l, 2*tau_points[ndat_-1]/beta-1)*gtau[ndat_-1]);*/
 
 	    //trapezoidal rule
     /*    double dtau = tau_points[1]-tau_points[0];
-	    gsum+=dtau*(bmth::legendre_p(l, 2*tau_points[0]/beta-1)*gtau[0]);
+	    gsum+=dtau*(legendre_polynomial(l, 2*tau_points[0]/beta-1)*gtau[0]);
 
         for(int i=1;i<ndat_-1;i++){
         	 double dtau = tau_points[i+1]-tau_points[i];
-	        gsum+=2*dtau*(bmth::legendre_p(l, 2*tau_points[i]/beta-1)*gtau[i]);
+	        gsum+=2*dtau*(legendre_polynomial(l, 2*tau_points[i]/beta-1)*gtau[i]);
         }
      	dtau = tau_points[ndat_-1]-tau_points[ndat_-2];
-        gsum+=dtau*(bmth::legendre_p(l, 2*tau_points[ndat_-1]/beta-1)*gtau[ndat_-1]);
+        gsum+=dtau*(legendre_polynomial(l, 2*tau_points[ndat_-1]/beta-1)*gtau[ndat_-1]);
         gsum/=2; */
         //try Simpson's rule within ndat-1 region, then trapezoidal for endpoint
 
@@ -122,25 +143,25 @@ double generateGl(const vector_type &gtau,const vector_type &tau_points, int l,d
         double gsum1=0.0,gsum2=0.0;
         ndat_-=2;
         double dtau = tau_points[1]-tau_points[0];
-    	gsum1+=gtau[0]*bmth::legendre_p(l,2*tau_points[0]/beta-1.0);
+        gsum1+=gtau[0]*legendre_polynomial(l,2*tau_points[0]/beta-1.0);
 	    for(int i=1;i<ndat_/2;i++)
-	    	gsum1+=2*gtau[2*i]*bmth::legendre_p(l,2*tau_points[2*i]/beta-1);
+        gsum1+=2*gtau[2*i]*legendre_polynomial(l,2*tau_points[2*i]/beta-1);
 	    for(int i=1;i<ndat_/2+1;i++)
-	    	gsum1+=4*gtau[2*i-1]*bmth::legendre_p(l,2*tau_points[2*i-1]/beta-1);
+        gsum1+=4*gtau[2*i-1]*legendre_polynomial(l,2*tau_points[2*i-1]/beta-1);
 
-	    gsum1+=gtau[ndat_]*bmth::legendre_p(l,2*tau_points[ndat_]/beta-1.0);
+	    gsum1+=gtau[ndat_]*legendre_polynomial(l,2*tau_points[ndat_]/beta-1.0);
 	    gsum1*= dtau/3;
 
       //-----------------------
       ndat_++;
       dtau = tau_points[1]-tau_points[0];
-      gsum2+=gtau[0]*bmth::legendre_p(l,2*tau_points[0]/beta-1.0);
+      gsum2+=gtau[0]*legendre_polynomial(l,2*tau_points[0]/beta-1.0);
       for(int i=2;i<ndat_/2;i++)
-        gsum2+=2*gtau[2*i]*bmth::legendre_p(l,2*tau_points[2*i]/beta-1);
+        gsum2+=2*gtau[2*i]*legendre_polynomial(l,2*tau_points[2*i]/beta-1);
       for(int i=2;i<ndat_/2+1;i++)
-        gsum2+=4*gtau[2*i-1]*bmth::legendre_p(l,2*tau_points[2*i-1]/beta-1);
+        gsum2+=4*gtau[2*i-1]*legendre_polynomial(l,2*tau_points[2*i-1]/beta-1);
 
-      gsum2+=gtau[ndat_]*bmth::legendre_p(l,2*tau_points[ndat_]/beta-1.0);
+      gsum2+=gtau[ndat_]*legendre_polynomial(l,2*tau_points[ndat_]/beta-1.0);
       gsum1*= dtau/3;
       //----------------------
       gsum = (gsum1+gsum2)/2;
@@ -155,7 +176,7 @@ double Gt(const vector_type &gl_in, double tau, double beta ){
     double gsum = 0.0;
     int N = gl_in.size();
     for(int l=0;l<N;l++){
-    	gsum+=std::sqrt(2*l+1)/beta*bmth::legendre_p(l,2.0*tau/beta-1.0)*gl_in[l];
+        gsum+=std::sqrt(2*l+1)/beta*legendre_polynomial(l,2.0*tau/beta-1.0)*gl_in[l];
     }
     return gsum;
 }
@@ -241,24 +262,6 @@ struct argstruct{
         double beta;int l;double tau; 
         std::pair<double,double> tails; vector_type *tau_points;
  };
-///genereate a normally distrubted noisy vector. 
-//i.e. output[i] = normally dist number with mean=data[i] and stddev=err[i]
-vector_type generateGaussNoise(vector_type data, vector_type err,boost::mt19937 &rng){
-
-    typedef boost::variate_generator<boost::mt19937&,boost::normal_distribution<> > ran_gen;
-    //notice the & in the first template argument and function rng argument.
-    //If we omit this, it will compile and run
-    //however, the numbers will be less(/not) random b/c it will copy the generator
-    //each time, outputting the mean with some noise, rather than truly random
-
-    const int N = data.size();
-    vector_type data_noise(N);
-    for(int i=0;i<N;i++){
-        boost::normal_distribution<> s(data[i],err[i]);
-        data_noise[i] = ran_gen(rng,s)();
-    }
-    return data_noise;
-}
 ///Generalized bootstrap routine. Requires the non-linear function to be
 // f(vector_type v,void *arg) where *arg is most easily a struct
 return_type bootstrap(double (*f)(vector_type,void*),
@@ -268,7 +271,7 @@ return_type bootstrap(double (*f)(vector_type,void*),
     //and determinging the variation on the output
     std::vector<double> newData(maxit);
     std::cout << std::setprecision(14);
-    boost::mt19937 rng;
+    std::mt19937 rng;
     rng.seed(static_cast<unsigned int>(std::time(0)));
     for(int i=0;i<maxit;i++){
          vector_type temp_data= generateGaussNoise(data, err,rng);

test/CMakeLists.txt

@@ -28,6 +28,7 @@ set(unit_tests
   backcontTest
   default_modelTest
   gridTest
+  legendreConvertTest
   paramFailureTest
   paramsTest
   simulationTest)

test/legendreConvertTest.cpp

@@ -0,0 +1,38 @@
+/*****************************************************************************
+*
+* ALPS Project Applications
+*
+* Copyright (C) 1998-2016 ALPS Collaboration
+*
+* ALPS Project: https://alps.comp-phys.org/
+* SPDX-License-Identifier: MIT
+*
+*****************************************************************************/
+
+#include "../legendre_convert/gaussian_noise.hpp"
+
+#include <gtest/gtest.h>
+
+#include <cmath>
+#include <random>
+#include <vector>
+
+TEST(LegendreConvert, ZeroDeviationDoesNotAdvanceGenerator) {
+  std::mt19937 rng(1234);
+  const std::mt19937 untouched = rng;
+  const std::vector<double> data{1.0, -2.0};
+
+  EXPECT_EQ(generateGaussNoise(data, {0.0, 0.0}, rng), data);
+  EXPECT_EQ(rng, untouched);
+}
+
+TEST(LegendreConvert, PositiveDeviationUsesDistribution) {
+  std::mt19937 rng(1234);
+  const std::mt19937 untouched = rng;
+
+  const auto noisy = generateGaussNoise({1.0}, {0.5}, rng);
+
+  ASSERT_EQ(noisy.size(), 1U);
+  EXPECT_TRUE(std::isfinite(noisy[0]));
+  EXPECT_NE(rng, untouched);
+}

test/regression/CMakeLists.txt

@@ -53,6 +53,20 @@ if(MAXENT_BUILD_UTILITIES)
             ${CMAKE_CURRENT_SOURCE_DIR}/reference ${results} --sets kk)
   set_tests_properties(regression_kk_generate PROPERTIES FIXTURES_SETUP regression_kk LABELS regression-fast)
   set_tests_properties(regression_kk_compare PROPERTIES FIXTURES_REQUIRED regression_kk LABELS regression-fast)
+
+  set(results ${CMAKE_CURRENT_BINARY_DIR}/results-legendre-convert)
+  add_test(NAME regression_legendre_convert_generate
+    COMMAND ${Python3_EXECUTABLE} ${CMAKE_CURRENT_SOURCE_DIR}/generate.py
+            --maxent $<TARGET_FILE:maxent>
+            --legendre-convert $<TARGET_FILE:legendre_convert>
+            --out ${results} --sets legendre_convert)
+  add_test(NAME regression_legendre_convert_compare
+    COMMAND ${Python3_EXECUTABLE} ${CMAKE_CURRENT_SOURCE_DIR}/compare.py
+            ${CMAKE_CURRENT_SOURCE_DIR}/reference ${results} --sets legendre_convert)
+  set_tests_properties(regression_legendre_convert_generate PROPERTIES
+    FIXTURES_SETUP regression_legendre_convert LABELS regression-fast)
+  set_tests_properties(regression_legendre_convert_compare PROPERTIES
+    FIXTURES_REQUIRED regression_legendre_convert LABELS regression-fast)
 endif()
 if(MAXENT_REGRESSION_FULL)
   maxent_regression_test(full full regression-full)

test/regression/MANIFEST.md

@@ -85,6 +85,12 @@ Generated by `manifest.py` from `reference/`, `cases.py` and `tolerances.json`.
 | `kk_imag_to_real_self` | kk: Im Sigma(omega) of self_u1 (average spectrum) (`--input_file=input.dat --output_file=output.dat --imag_to_real`) | `test/regression/inputs/kk_imag_to_real_self` | 6 | 1.1 s |  |
 | `kk_real_to_imag` | kk: Re G(omega) of a unit Gaussian spectrum (Dawson function) (`--input_file=input.dat --output_file=output.dat --real_to_imag`) | `test/regression/inputs/kk_real_to_imag` | 6 | 1.0 s |  |
 
+## legendre_convert utility
+
+| Case | Covers | Inputs | Datasets | Runtime | Flags |
+|---|---|---|---|---|---|
+| `legendre_convert_transform` | legendre_convert: transform, tails, back-continuation, and Matsubara convergence (`--beta=5 --input_gtau_file=input.dat --output_gl_file=Gl.dat --tail1=1 --tail2=0 --lmax=12 --maxit=10 --backcontinue --plot_convergence=12 --maxn=3 --noerr`) | `test/regression/inputs/legendre_convert_transform` | 9 | 0.3 s |  |
+
 ## Components
 
 | Case | Covers | Inputs | Datasets | Runtime | Flags |
@@ -97,6 +103,7 @@ Rules are matched in order against '<case>:<dataset>'; the first match applies.
 
 | Pattern (case:dataset) | atol | rtol | Why |
 |---|---|---|---|
+| `^legendre_convert_transform:` | 1e-14 | 1e-12 | direct utility regression; allows only rounding-level variation |
 | `:cli/` | 0 | 0 | command-line output and exit code must not change |
 | `:log/scalars/(n_singular|max_it_warnings)$` | 0 | 0 | counts; identical in all variants |
 | `:log/scalars/minimal_chi2$` | 1e-12 | 1e-06 | values below ~1e-12 are rounding noise (observed relative differences up to 1 at 1e-17) |
@@ -112,7 +119,7 @@ Rules are matched in order against '<case>:<dataset>'; the first match applies.
 | `\.chi2\.dat$` | 0 | 1e-07 | observed 1.7e-9 |
 | `:files/` | 1e-15 | 2e-06 | spectra, back-continuations, grids, all other outputs; observed 1.5e-7 (GCC 15, u2_time_fast) |
 
-Total compared datasets: 1183 in 59 cases.
+Total compared datasets: 1192 in 60 cases.
 
 ## Known limitations
 

test/regression/PROVENANCE.md

@@ -45,6 +45,15 @@ form of the Legendre kernel reproduces the GSL kernels in `components.h5` to
 3.3e-16 relative, all Legendre cases pass, and the new `kk` spline reproduces
 the `kk` references bit for bit.
 
+## legendre_convert reference (added 2026-09-24)
+
+`legendre_convert_transform.h5` was generated from `modernize/step2` merge
+commit `a57f150` using AppleClang 21, C++17, Boost 1.88, and ALPSCore
+`560ae112`. The utility source was unchanged from the Boost implementation.
+Its synthetic Gaussian-spectrum input has zero error bars, so the clock-seeded
+bootstrap is deterministic; three consecutive generate-and-compare runs
+passed.
+
 ## Determinism
 
 Two independent generations with the reference build were compared with

test/regression/README.md

@@ -33,6 +33,7 @@ result, and is flagged in [`MANIFEST.md`](MANIFEST.md).
 | `cli` | `--help`, `--help.models`, `--help.grids`, missing `BETA` | instant | command-line behavior |
 | `full` | the 11 example runs as shipped | about 1 min (Legendre dominates) | opt-in |
 | `kk` | the `kk` utility on 3 inputs (Im G and Im Σ from the references, an analytic Re G) | seconds | registered when `MAXENT_BUILD_UTILITIES=ON` |
+| `legendre_convert` | deterministic transform, tails, back-continuation and Matsubara convergence | instant | protects the utility before Boost replacements |
 | components | `components.h5` from `dump_components` | instant | building blocks |
 
 ## Running
@@ -52,9 +53,10 @@ them against `reference/`. By hand:
 ```bash
 python3 test/regression/generate.py --maxent build/maxent \
     --components build/test/regression/dump_components --out /tmp/maxent-results \
-    [--kk build/kk/kk] [--sets fast,targeted,cli,full,kk] [--cases REGEX]
+    [--kk build/kk/kk] [--legendre-convert build/legendre_convert/legendre_convert] \
+    [--sets fast,targeted,cli,full,kk,legendre_convert] [--cases REGEX]
 python3 test/regression/compare.py test/regression/reference /tmp/maxent-results \
-    [--sets fast,targeted,cli,kk,components] [--cases REGEX]
+    [--sets fast,targeted,cli,kk,legendre_convert,components] [--cases REGEX]
 ```
 
 `compare.py` exits with status 1 if any case fails. `--report` prints the