std::cyl_neumann, std::cyl_neumannf, std::cyl_neumannl

#include <cassert>
#include <cmath>
#include <iostream>
#include <numbers>
 
const double π = std::numbers::pi; // or std::acos(-1) in pre C++20
 
// Para calcular la función cilíndrica de Neumann vía la función cilíndrica de Bessel de
// primera especie tenemos que implementar la J₋ᵥ, ya que la invocación directa de
// std::cyl_bessel_j(ν,x), de acuerdo con la fórmula anterior, genera 'std::domain_error'
// para una ν negativa.
// Mal argumento en __cyl_bessel_j.
 
double J₋ᵥ (double ν, double x) {
    return std::cos(-ν*π) * std::cyl_bessel_j(-ν,x)
          -std::sin(-ν*π) * std::cyl_neumann(-ν,x);
}
 
double J₊ᵥ (double ν, double x) { return std::cyl_bessel_j(ν,x); }
 
double Jᵥ (double ν, double x) { return ν < 0.0 ? J₋ᵥ(ν,x) : J₊ᵥ(ν,x); }
 
int main()
{
    std::cout << "verificaciones al azar para ν == 0.5\n" << std::fixed << std::showpos;
    double ν = 0.5;
    for (double x = 0.0; x <= 2.0; x += 0.333) {
        const double n = std::cyl_neumann(ν, x);
        const double j = (Jᵥ(ν, x)*std::cos(ν*π) - Jᵥ(-ν, x)) / std::sin(ν*π);
        std::cout << "N_.5(" << x << ") = " << n << ", calculado vía J = " << j << '\n';
        assert(n == j);
    }
}