Just dropping this on here for documentation/searchability. I do not think we should merge this functionality. In fact, I do not think the paper that introduces this solver strategy should have been published at all.
The TL;DR is that the scalars used to determine step sizes are literally unrepresentable even for fairly loose tolerances.
Pasting in my commentary from the implementation as detailed rationale for trashing the whole solver strategy.
//! Due to the use of the term `1 / (1 + err^2)^0.5`, the solver update here will
//! crash any time the error is less than around sqrt(<T>::EPSILON), where epsilon
//! is around 1e-16 for f64 and around 1e-7 for f32. Empirically, this sets a practical
//! lower bound on tolerances of about 1e-6 for f64 and about 0.1 for f32 after which
//! the solver will produce NaN values regardless of the quality of the inputs.
I also took a look at using Decimal or BigRational types to circumvent this limitation, since it's only necessary for a small number of scalar calculations. Unfortunately, because the resulting rational is still unrepresentable as a float, it makes no meaningful difference.
Just dropping this on here for documentation/searchability. I do not think we should merge this functionality. In fact, I do not think the paper that introduces this solver strategy should have been published at all.
The TL;DR is that the scalars used to determine step sizes are literally unrepresentable even for fairly loose tolerances.
Pasting in my commentary from the implementation as detailed rationale for trashing the whole solver strategy.
I also took a look at using Decimal or BigRational types to circumvent this limitation, since it's only necessary for a small number of scalar calculations. Unfortunately, because the resulting rational is still unrepresentable as a float, it makes no meaningful difference.