pfqn_grnmol.h File Reference

Normalizing constant by the closed-form Grundmann-Moeller rule. More...

#include <cstddef>
#include <vector>
#include "line/num/number.h"
#include "line/util/error.h"
#include "line/util/matrix.h"

Include dependency graph for pfqn_grnmol.h:

Go to the source code of this file.

Namespaces
namespace	line
namespace	line::pfqn

Functions
template<class T>
T	line::pfqn::pfqn_grnmol (const Matrix< T > &L, const std::vector< int > &N)
	Normalizing constant by the closed-form Grundmann-Moeller rule.

Detailed Description

Normalizing constant by the closed-form Grundmann-Moeller rule.

Templated port of matlab/src/api/pfqn/pfqn_grnmol.m, which applies the rule directly rather than through the successive-degree recursion of pfqn_cub:

G = [ (sum N + M - 1)! / prod_r N_r! ] * sum_{i=0}^{S} w_i H_i w_i = 2^{-2S} (-1)^i c_i^{2S+1} / (i! (i + c_i)!), c_i = 2(S-i) + M H_i = sum_{|b| = S-i} prod_r ( ((2b+1)/c_i)' L(:,r) )^{N_r}

with b ranging over the M-vectors of non-negative integers summing to S-i, which is what the reference's call to matlab/util/sprod.m enumerates. The point (2b+1)/c_i is barycentric by construction, since sum_m (2b_m+1) = c_i.

REFERENCE DEFECT, reproduced as a rejection rather than as a wrong number. pfqn_grnmol.m sets S = ceil(sum(N)-1)/2, which MATLAB parses as (ceil(sum(N)-1))/2 and NOT as ceil((sum(N)-1)/2), the value pfqn_cub.m uses. For an EVEN total population S is therefore a half-integer, c_i is a half-integer, and the weight calls factorial(i + c_i) on a non-integer, which MATLAB rejects outright ("N must be a matrix of non-negative integers"). pfqn_grnmol is thus callable only for an ODD total population, and the port throws InputError for an even one instead of silently choosing one of the two readings of the expression. Use pfqn_cub for even populations; it computes the same integral with the correct degree.

ARITHMETIC. Every ingredient – the barycentric points, the binomial-like weights, the integer powers – is rational, and the factorial prefactor is formed as an exact factorial rather than through gammaln. The routine is therefore EXACT in rational arithmetic and is deliberately left ungated: at the reference degree it returns the exact normalizing constant, and checking that against pfqn_ca is the sharpest test the rule admits.

Definition in file pfqn_grnmol.h.