Incorrect asin

[libbooze]

I did not have a lot of time to investigate, so here is brief recap of what I figured out for now.

template<typename T>
void test_asin()
{
    using nl = std::numeric_limits<T>;
    std::print("{}\n", boost::typeindex::type_id<T>().pretty_name());
    auto tiny0 = nl::epsilon()*999/1000;
    auto tiny1 = nl::epsilon();
    auto tiny2 = nl::epsilon()*1000/999;
    assert(tiny0!=tiny1);
    assert(tiny1!=tiny2);
    std::println("{:.36e} -> {:.36e}", tiny0, asin(tiny0));
    std::println("{:.36e} -> {:.36e}", tiny1, asin(tiny1));
    std::println("{:.36e} -> {:.36e}", tiny2, asin(tiny2));
}

prints

float
1.190900817960027779918164014816284180e-07 -> 1.190900817960027779918164014816284180e-07
1.192092895507812500000000000000000000e-07 -> 1.192092895507812500000000000000000000e-07
1.193286180978248012252151966094970703e-07 -> 1.193286180978248012252151966094970703e-07
double
2.218225603201062765794263809792929471e-16 -> 2.218225603201062765794263809792929471e-16
2.220446049250313080847263336181640625e-16 -> 2.220446049250313080847263336181640625e-16
2.222668717968281558635048558466794015e-16 -> 2.222668717968281558635048558466794015e-16
long double
1.083117970313018929550875856533283850e-19 -> 1.083117970313018929550875856533283850e-19
1.084202172485504434007452800869941711e-19 -> 1.084202172485504434007452800869941711e-19
1.085287459945449883881930776840203615e-19 -> 1.085287459945449883881930776840203615e-19
boost::decimal::decimal32_t
9.990000000000000000000000000000000000e-07 -> 9.990000000000000000000000000000000000e-07
1.000000000000000000000000000000000000e-06 -> 1.000000000000000000000000000000000000e-06
1.001001000000000000000000000000000000e-06 -> 1.001075000000000000000000000000000000e-06
boost::decimal::decimal64_t
9.990000000000000000000000000000000000e-16 -> 9.990000000000000000000000000000000000e-16
1.000000000000000000000000000000000000e-15 -> 1.000000000000000000000000000000000000e-15
1.001001001001001000000000000000000000e-15 -> 1.001609728980325000000000000000000000e-15
boost::decimal::decimal128_t
9.990000000000000000000000000000000000e-34 -> 9.990000000000000000000000000000000000e-34
1.000000000000000000000000000000000000e-33 -> 1.000000000000000000000000000000000000e-33
1.001001001001001001001001001001001000e-33 -> 1.109451248382801572892973052948674000e-33

I think that depending on precision you guarantee you could claim all results are fine, but I do not know what you guarantee so I opened an issue.

Issues I can see:
https://www.wolframalpha.com/input?i=asin%281.001001001001001000000000000000000000e-15+%29
https://www.wolframalpha.com/input?i=asin%281.001001001001001001001001001001001000e-33+%29

In other words tiny2 for decimal64 and decimal128 seems to have very large error.
Beside that function does not seem "smooth". What I mean by this is that code jumps from identity(tiny0, tiny1) to values that are quite higher from identity for small increment over epsilon since different code branch is used to compute tiny0, tiny1 versus branch for tiny2.

[mborland]

@ckormanyos when you have time could you take a look at this one? Right now the control logic is:

if (absx <= std::numeric_limits<T>::epsilon())
{
    result = absx;
}
else if (absx <= T{5, -1})
{
    result = asin_series(absx);
}
else if (absx <= T{1, 0})
{
    constexpr T half_pi {numbers::pi_v<T> / 2};
    result = half_pi - 2 * asin_series(sqrt((1 - absx) / 2));
}

I used a remez polynomial from 0, 0.5 that should have enough terms for low error, but maybe there needs to be another expansion at the extremes?

[ckormanyos]

Thx @libbooze for pointing this out. I think there is room for improvement here.

@ckormanyos when you have time could you take a look at this one?

Yes @mborland, I am happy to go through this one.

Sometimes for small values, it can be wise to include several terms of a Taylor-expansion.

Give me a couple days to find the right approximation(s) and their crossovers. Should be straightforward. Usually it is just a matter of bringing these pesky edge-cases into focus and then dealing with them.

[ckormanyos]

The intended repair is cycling in #1145. We might need similar treatment of args near epsilon for other funcs as the library evolves.

I put in a check for args less than cube root of epsilon.

Cc: @libbooze and @mborland

[ckormanyos]

See also:

Series[ArcSin[x], {x, 0, 3}]

$$ \sin^{-1}{x} {\approx} x + \dfrac{x^3}{6} $$

for $x$ small.