We might not even need numeric_limits_float128.h - the generic ones might work.

Debye integrals. matharr has some. Started a test.

We test double pretty hard. That is where most of the external data is.
We need to test using 32-bit float.
Find where it craps out. Make sure our limits-defined numbers work over a range of types.

Look at the results of test_bessel, test_sph_bessel, near the originas a function of.
You have to widen the near zero. Look at the near zero asymptotics and come up with a thing
as a function of underflow exponent and nu.
Neumann and K are just going to blow.
Should I go ahead and let J, I return denormal results? Why not.
And if the type doesn't do denorm just zero then and don't sweat.

Finish experfc = exp(x*x)erfc(x) (and maybe even sqrt(pi)xexp(xx)*erfc(x))
Call it erfc_scaled?

expint_en may be just sloppy. I think recursion may be bad?
Boost uses CF + series.

Implement polygamma.

We have both psi and digamma! Pick one!

Should I do the /_GLIBCXX_END_NAMESPACE_VERSION thing for __gnu_cxx?
Yes!

Things involving say std::cyl_bessel(__float128, __float128) won't compile.
We have __float128 in the promote thing.
Somehow, the promote thing is broken.
These were using the old school promote!

Lose __eps, __min, __max, __inf, __NaN from math_const in favor of numeric_limits.h

heuman_lambda - I've got this! - __heuman_lambda(_Tp __beta, _Tp __k)
template <class T1, class T2>
calculated-result-type heuman_lambda(T1 k, T2 phi);
and __jacobi_zeta(_Tp __k, _Tp __phi)
template <class T1, class T2, class Policy>
calculated-result-type jacobi_zeta(T1 k, T2 phi, const Policy&);
Done - __ellint_d(_Tp __k, _Tp __phi)
template <class T1, class T2>
calculated-result-type ellint_d(T1 k, T2 phi);

Implement inverse erf and erfc.
And for that matter, we want native implementations of erf and erfc as well. And for the complex plane too.

Test against spherical_harmonic(unsigned n, int m, T1 theta, T2 phi);
They (GSL) have spherical_harmonic_[r|i]
Note the signed m on all of them!

Do the thing.
Actually, it's done.

Why did check_lfactorial.cc get chopped after 171?
{ 711.71472580228999, 171, 0.0 },
- { 716.86222027910344, 172, 0.0 },

Add a state for Fock-airy functions and use it.

__riemann_zeta_alt should be tested with vanWijnGaarden.
These transforms work for other similar sums.

polynomial.h needs outline polynomial.tcc
Many algos are long.

Implement a static compile-time polynomial.
LambertW has a very good static polynomial I think - with static Horner.
They use macros... Sigh. Still, we could hack these into templates.

__log_gamma was a very sobering exercise in trig accuracy.
The error in the location of the zeros of, say,
std::sin(x) grow linearly with argument!
I had to test epsilon * |x| for zero to place poles!
The world really needs reperiodized trig functions!

sph_neumann is broken.

Tabulate riemann_zeta for integer arguments - various sums can use this.

Do the thing!

We will want lgamma_scaled - the Binet function.
Then tgamma_scaled - sometines called \Gamma^{*}.
Then Amparo, Gil, etc have a paper...

Implement inverse erf and erfc.
For that matter, I want native implementations of erf and erfc as well.

We test double pretty hard. That is where most of the external data is.
We need to test using long double and if different, __float128.

Owen's T function (a stat function)
template
calculated-result-type owens_t(T h, T a);

Expose Bernoulli numbers.

We should make a thing that returns the complement of ibeta.
Internally, we could return both in a struct.

Building and execution work but there are other general issues.

In mathar it turns out that
(2\pi)^k B_{2k}/(2k)! - 2 is small. It might be worth computing... carefully.
I think I just did this in __riemann_zeta_m_1_sum.
I might want a starting fit for small s?

I seem to have lost accuracy in cyl_bessel_i near the origin.
cyl_neumann also for < 2! All bessels in non-oscillatory region.

Gather 32-digit spot checks from:
http://functions.wolfram.com/webMathematica/FunctionEvaluation.jsp

__float128 should be a good foil.

I got derivs in the Bessel asymp!
Now I need a more principled max iteration.

W now ship bernoulli numbers, etc.

test_hankel build fails on
test_hankel.cpp:52:47: error: no matching function for call to '__airy(std::complex&, double)'
auto Airy = std::__detail::__airy(z, 1.0e-16);
It builds.
It's probably an API change.

Why did dilog lose the x=0 test?
It is because the answer dilog(0) == 0. So the fractional error blows up even when the answer is spot on.
But it doesn't look like such a point would be skipped though.
We need a way of trapping this...
We need a way of deciding when a point approximately crosses the real axis.
We don't want to trap when points asymptotically approach the axis.
So we need to store some previous numbers for a derivative approximation.
Maybe if |f| < eps and |f - fprev| > N eps switch to absolute diff (and mark it).

In test_riemann_zeta the negative argument differences between real zetac an zeta are much worse.
Reflection formula and reperiodized stuff?

Look at test_cyl_hankel - we got issues! vu=1, z=0.01 Y bombs on Steeds algorithm!
OK! Steeds algorithm fails for cyl_neumann for exatly nu=1 and no other. WTF!
nearbyint(nu + .5) -> nearbyint(nu). mu = nu - n was -1.

template <class T>
  calculated-result-type owens_t(T h, T a);

Boost has this.

Test asymptotic Hermite polys - adjust the crossover width!
GSL extras has these - there was struggle there too.

radpoly is nan ing out for odd m!

Implement inverse incomplete gammas:
[] We will want lgamma_scaled - the Binet function.
[] Then tgamma_scaled - sometines called \Gamma^{*}.
[] Then Amparo, Gil, etc have a paper...

beta for a=100 b>71 wrong sign. But dead on magnitude! We could crush it!

Boost has this.

Implement heuman_lambda - I've got this! - __heuman_lambda(_Tp __beta, _Tp __k)
template <class T1, class T2>
calculated-result-type heuman_lambda(T1 k, T2 phi);
and __jacobi_zeta(_Tp __k, _Tp __phi)
template <class T1, class T2, class Policy>
calculated-result-type jacobi_zeta(T1 k, T2 phi, const Policy&);
Done - __ellint_d(_Tp __k, _Tp __phi)
template <class T1, class T2>
calculated-result-type ellint_d(T1 k, T2 phi);

In testcase.tcc in the big string_view header I'd like to be able to replace
#include <specfun_testcase.h>
with
#include "../testcase.h"
for TR1 test generation...

Use dammit!
I'll wait until it shows up in std in libstdc++.

ibeta looks like it might have problems for large a, b near end points.
Check test_beta_inc.

bincoeff can take real upper thing - infinitely many k
Think (1 + x)^{1/2} as a generator.