boost spheric function

main.cpp

#include <boost/math/special_functions/spherical_harmonic.hpp>
#include <boost/math/special_functions/legendre.hpp>
#include <boost/math/constants/constants.hpp>

#include <cmath>
#include <iostream>
#include <type_traits>
#include <limits>

using namespace boost::math;

template<class T> inline T legendre_sp_00() {
	return static_cast<T>(1)/boost::math::constants::root_pi<T>()/2;
}
template<class T> inline T legendre_sp_next_mm(std::size_t m, T w, T sp_mm) {
	/* w \equiv \sqrt{1-x^2} */
	using CT = typename std::common_type<T, std::size_t, std::uintmax_t>::type;
	return -w * sp_mm * std::sqrt(static_cast<CT>(1) + static_cast<CT>(1)/static_cast<CT>(m)/static_cast<CT>(2));
}
template<class T> inline T legendre_sp_next_ml(std::size_t m, std::size_t l, T x, T sp_1, T sp_0) {
	using CT = typename std::common_type<T, std::size_t, std::uintmax_t>::type;
	const std::uintmax_t lpm = l+m;
	const std::size_t    lmm = l-m;
	const CT mp12 = static_cast<CT>(m)+static_cast<CT>(1)/static_cast<CT>(2);
	const CT lm32 = static_cast<CT>(l)-static_cast<CT>(3)/static_cast<CT>(2);
	const CT c1 = std::sqrt((static_cast<CT>(1)+mp12/lmm)*(static_cast<CT>(1)-mp12/lpm));
	const CT c2 = std::sqrt((static_cast<CT>(1)-static_cast<CT>(1)/lmm)*(static_cast<CT>(1)-static_cast<CT>(1)/lpm)*(static_cast<CT>(1)+static_cast<CT>(2)/lm32));
	return x * c1 * sp_1 * 2 - c2 * sp_0;
}

template<class T>
void print_sph(int max, T x) {
	static const T p_00 = legendre_sp_00<T>();
	T p_mm = p_00;
	T w = std::sqrt(1-x*x);

	std::cout << 0 << " " << 0 << " " << p_mm << std::endl;
	for(int m = 1; m <= max; ++m) {
		p_mm = legendre_sp_next_mm(m, w, p_mm);
		T p_0 = 0;
		T p_1 = p_mm;
		std::cout << m << " " << m << " " << p_1 << std::endl;
		for(int l = m + 1; l <= max; ++l) {
			const T p_l = legendre_sp_next_ml(m,l,x,p_1,p_0);
			p_0 = p_1;
			p_1 = p_l;
			std::cout << m << " " << l << " " << p_l << std::endl;

		}
	}
}

int main(int argc, char** argv) {

	print_sph(4, 0.75);

	return 0;
}