libcruft-util/test/maths/fast.cpp

#include <cruft/util/maths/fast.hpp>
#include <cruft/util/tap.hpp>

#include <iostream>
#include <iomanip>


int
main ()
{
    struct {
        unsigned sin = 0;
        unsigned exp = 0;
    } successes;

    // 64+1 evenly spaced evaluations of sine in the interval [0,pi/4]
    // calculated with higher than float64 precision using python's mpmath.
    //
    // >>> import mpmath
    // >>> steps=64
    // [mpmath.nstr(mpmath.sin(mpmath.pi/4*(x/steps)), 8) for x in range(steps+1)]
    static constexpr float sins[] = {
        0x0.000000p+00f, 0x1.921d1fp-07f, 0x1.92155fp-06f, 0x1.2d8657p-05f,
        0x1.91f65fp-05f, 0x1.f656e7p-05f, 0x1.2d5209p-04f, 0x1.5f6d00p-04f,
        0x1.917a6bp-04f, 0x1.c3785cp-04f, 0x1.f564e5p-04f, 0x1.139f0cp-03f,
        0x1.2c8106p-03f, 0x1.45576bp-03f, 0x1.5e2144p-03f, 0x1.76dd9dp-03f,
        0x1.8f8b83p-03f, 0x1.a82a02p-03f, 0x1.c0b826p-03f, 0x1.d934fep-03f,
        0x1.f19f97p-03f, 0x1.04fb80p-02f, 0x1.111d26p-02f, 0x1.1d3443p-02f,
        0x1.294062p-02f, 0x1.35410cp-02f, 0x1.4135c9p-02f, 0x1.4d1e24p-02f,
        0x1.58f9a7p-02f, 0x1.64c7ddp-02f, 0x1.708853p-02f, 0x1.7c3a93p-02f,
        0x1.87de2ap-02f, 0x1.9372a6p-02f, 0x1.9ef794p-02f, 0x1.aa6c82p-02f,
        0x1.b5d100p-02f, 0x1.c1249dp-02f, 0x1.cc66e9p-02f, 0x1.d79775p-02f,
        0x1.e2b5d3p-02f, 0x1.edc195p-02f, 0x1.f8ba4dp-02f, 0x1.01cfc8p-01f,
        0x1.073879p-01f, 0x1.0c9704p-01f, 0x1.11eb35p-01f, 0x1.1734d6p-01f,
        0x1.1c73b3p-01f, 0x1.21a799p-01f, 0x1.26d054p-01f, 0x1.2bedb2p-01f,
        0x1.30ff7fp-01f, 0x1.36058bp-01f, 0x1.3affa2p-01f, 0x1.3fed95p-01f,
        0x1.44cf32p-01f, 0x1.49a449p-01f, 0x1.4e6cabp-01f, 0x1.532829p-01f,
        0x1.57d693p-01f, 0x1.5c77bbp-01f, 0x1.610b75p-01f, 0x1.659192p-01f,
        0x1.6a09e6p-01f
    };

    // [hex(mpmath.exp(mpmath.log(2)/2/steps*x), 8) for x in range(steps+1)]
    static constexpr float exps[] = {
        0x1.000000p+00f, 0x1.0163dap+00f, 0x1.02c9a3p+00f, 0x1.04315ep+00f,
        0x1.059b0dp+00f, 0x1.0706b2p+00f, 0x1.087451p+00f, 0x1.09e3ecp+00f,
        0x1.0b5586p+00f, 0x1.0cc922p+00f, 0x1.0e3ec3p+00f, 0x1.0fb66ap+00f,
        0x1.11301dp+00f, 0x1.12abdcp+00f, 0x1.1429aap+00f, 0x1.15a98cp+00f,
        0x1.172b83p+00f, 0x1.18af93p+00f, 0x1.1a35bep+00f, 0x1.1bbe08p+00f,
        0x1.1d4873p+00f, 0x1.1ed502p+00f, 0x1.2063b8p+00f, 0x1.21f499p+00f,
        0x1.2387a6p+00f, 0x1.251ce4p+00f, 0x1.26b456p+00f, 0x1.284dfep+00f,
        0x1.29e9dfp+00f, 0x1.2b87fdp+00f, 0x1.2d285ap+00f, 0x1.2ecafap+00f,
        0x1.306fe0p+00f, 0x1.32170fp+00f, 0x1.33c08bp+00f, 0x1.356c55p+00f,
        0x1.371a73p+00f, 0x1.38cae6p+00f, 0x1.3a7db3p+00f, 0x1.3c32dcp+00f,
        0x1.3dea64p+00f, 0x1.3fa450p+00f, 0x1.4160a2p+00f, 0x1.431f5dp+00f,
        0x1.44e086p+00f, 0x1.46a41ep+00f, 0x1.486a2bp+00f, 0x1.4a32afp+00f,
        0x1.4bfdadp+00f, 0x1.4dcb29p+00f, 0x1.4f9b27p+00f, 0x1.516daap+00f,
        0x1.5342b5p+00f, 0x1.551a4cp+00f, 0x1.56f473p+00f, 0x1.58d12dp+00f,
        0x1.5ab07dp+00f, 0x1.5c9268p+00f, 0x1.5e76f1p+00f, 0x1.605e1bp+00f,
        0x1.6247ebp+00f, 0x1.643463p+00f, 0x1.662388p+00f, 0x1.68155dp+00f,
        0x1.6a09e6p+00f,
    };
    (void)exps;

    cruft::TAP::logger tap;

    for (std::size_t i = 0; i < std::size (sins); ++i) {
        const auto theta = cruft::pi<float> / 4 * static_cast<float> (i) / (std::size (sins) - 1);
        const auto value = cruft::maths::fast::sin (theta);

        const auto abserr = std::abs (sins[i] - value);

        successes.sin += (abserr < 1e-32f || cruft::relatively_equal (sins[i], value, 1e-6f)) ? 1 : 0;
    }

    for (std::size_t i = 0; i < std::size (exps); ++i) {
        static constexpr float half_log2 = 0.34657359027997265470861606072908828403775006718012f;
        const auto exponent = half_log2 * static_cast<float> (i) / (std::size (exps)-1);
        const auto value = cruft::maths::fast::exp (exponent);

        const auto abserr = std::abs (exps[i] - value);
        successes.exp += (abserr < 1e-32f || cruft::relatively_equal (exps[i], value, 1e-4f)) ? 1 : 0;
    }

    //for (int i = 0; i < 64; ++i) {
    //    const float expected = expf (i);
    //    const float observed = cruft::maths::fast::exp (i);
    //
    //    const float relerr = 1 - observed / expected;
    //    std::cout << observed << '\t' << expected << '\t' << relerr << '\n';
    //    tap.expect_lt (relerr, 1e-4f, "relative exp({}) error", i);
    //}

    tap.expect_eq (successes.sin, std::size (sins), "sin evaluation, {}/{}", successes.sin, std::size (sins));
    tap.expect_eq (successes.exp, std::size (exps), "exp evaluation, {}/{}", successes.exp, std::size (exps));

    return tap.status ();
}