C++ file that defines various crude approximations of floating point functions.

fast_float_funcs.hpp

/**
 * MIT License
 * 
 * Copyright (c) 2021 Sam Belliveau
 * 
 * Permission is hereby granted, free of charge, to any person obtaining a copy
 * of this software and associated documentation files (the "Software"), to deal
 * in the Software without restriction, including without limitation the rights
 * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
 * copies of the Software, and to permit persons to whom the Software is
 * furnished to do so, subject to the following conditions:
 * 
 * The above copyright notice and this permission notice shall be included in all
 * copies or substantial portions of the Software.
 * 
 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
 * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
 * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
 * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
 * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
 * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
 * SOFTWARE.
 */

#ifndef SAM_BELLIVEAU_FAST_FLOATING_POINT_FUNCS
#define SAM_BELLIVEAU_FAST_FLOATING_POINT_FUNCS 1

#include <bit>
#include <limits>
#include <cstdint>
#include <cmath>

/**
 * Context:
 * 
 * This code is heavily inspired by the famous Q_rsqrt(float) 
 * function from quake. It offers some other definitions of 
 * functions that do similar things in both 32 bit and 64 bit
 * versions.
 * 
 * Floating point numbers are really like, an exponential version
 * of integers when you really get into it. Many of the different
 * equations are various log identities, such as:
 * 
 *      2 * log(x) = x ^ 2
 *      -log(x) = 1 / x
 *      0.5 * log(x) = sqrt(x)
 *      -0.5 * log(x) = 1 / sqrt(x)
 * 
 * By using the fact that casting a float to an int is similarish
 * to taking a logarithm, and casting it back is sort of like reversing
 * that, we can take advantage of these identities to make crude
 * approximations.
 */

/**
 * Resources:
 * 
 *  - Wikipedia:
 *      - https://en.wikipedia.org/wiki/Fast_inverse_square_root
 * 
 *  - Desmos (created by yours truly):
 *      - https://www.desmos.com/calculator/2yhvcg8nxn
 *      - This desmos really helps visualize whats happening
 */

namespace sb 
{

    /********************/
    /*** FAST SQUARES ***/
    /********************/

    [[nodiscard]] inline constexpr 
    float fast_square32(const float x) noexcept
    {
        static_assert(
            std::numeric_limits<float>::is_iec559, 
            "32-bit Floats on this architecture are not IEC559 compliant!"
        );

        constexpr std::uint32_t MAGIC_ADJ = (std::uint32_t(1) << 30) - (std::uint32_t(1) << 23);
        
        return std::bit_cast<float>(
            (std::bit_cast<std::uint32_t>(x) << 1) - MAGIC_ADJ
        );
    }

    [[nodiscard]] inline constexpr 
    double fast_square64(const double x) noexcept
    {
        static_assert(
            std::numeric_limits<double>::is_iec559, 
            "64-bit Floats on this architecture are not IEC559 compliant!"
        );

        constexpr std::uint64_t MAGIC_ADJ = (std::uint64_t(1) << 62) - (std::uint64_t(1) << 52);
        
        return std::bit_cast<double>(
            (std::bit_cast<std::uint64_t>(x) << 1) - MAGIC_ADJ
        );
    }


    /*************************/
    /*** FAST SQUARE ROOTS ***/
    /*************************/

    [[nodiscard]] inline constexpr 
    float fast_sqrt32(const float x) noexcept 
    {
        static_assert(
            std::numeric_limits<float>::is_iec559, 
            "32-bit Floats on this architecture are not IEC559 compliant!"
        );

        constexpr std::uint32_t MAGIC_ADJ = (std::uint32_t(1) << 29) - (std::uint32_t(1) << 22);

        return std::bit_cast<float>(
            MAGIC_ADJ + (std::bit_cast<std::uint32_t>(x) >> 1)
        );
    }

    [[nodiscard]] inline constexpr 
    double fast_sqrt64(const double x) noexcept 
    {
        static_assert(
            std::numeric_limits<double>::is_iec559, 
            "64-bit Floats on this architecture are not IEC559 compliant!"
        );

        constexpr std::uint64_t MAGIC_ADJ = (std::uint64_t(1) << 61) - (std::uint64_t(1) << 51);

        return std::bit_cast<double>(
            MAGIC_ADJ + (std::bit_cast<std::uint64_t>(x) >> 1)
        );
    }


    /***********************/
    /*** FAST RECIPROCAL ***/
    /***********************/

    [[nodiscard]] inline constexpr 
    float fast_rcpl32(const float x) noexcept 
    {
        static_assert(
            std::numeric_limits<float>::is_iec559, 
            "32-bit Floats on this architecture are not IEC559 compliant!"
        );

        constexpr std::uint32_t MAGIC_ADJ = (std::uint32_t(1) << 31) - (std::uint32_t(1) << 24);

        return std::bit_cast<float>(
            MAGIC_ADJ - std::bit_cast<std::uint32_t>(x)
        );
    }

    [[nodiscard]] inline constexpr 
    double fast_rcpl64(const double x) noexcept 
    {
        static_assert(
            std::numeric_limits<double>::is_iec559, 
            "64-bit Floats on this architecture are not IEC559 compliant!"
        );

        constexpr std::uint64_t MAGIC_ADJ = (std::uint64_t(1) << 63) - (std::uint64_t(1) << 53);

        return std::bit_cast<double>(
            MAGIC_ADJ - std::bit_cast<std::uint64_t>(x)
        );
    }


    /*********************************/
    /*** FAST INVERSE SQUARE ROOTS ***/
    /*********************************/

    [[nodiscard]] inline constexpr 
    float fast_rcpl_sqrt32(const float x) noexcept 
    {
        static_assert(
            std::numeric_limits<float>::is_iec559, 
            "32-bit Floats on this architecture are not IEC559 compliant!"
        );

        constexpr std::uint32_t MAGIC_ADJ = (std::uint32_t(1) << 31) - (std::uint32_t(1) << 29) - (std::uint32_t(1) << 23) - (std::uint32_t(1) << 22);

        return std::bit_cast<float>(
            MAGIC_ADJ - (std::bit_cast<std::uint32_t>(x) >> 1)
        );
    }

    [[nodiscard]] inline constexpr 
    double fast_rcpl_sqrt64(const double x) noexcept 
    {
        static_assert(
            std::numeric_limits<double>::is_iec559, 
            "64-bit Floats on this architecture are not IEC559 compliant!"
        );

        constexpr std::uint64_t MAGIC_ADJ = (std::uint64_t(1) << 63) - (std::uint64_t(1) << 61) - (std::uint64_t(1) << 52) - (std::uint64_t(1) << 51);

        return std::bit_cast<double>(
            MAGIC_ADJ - (std::bit_cast<std::uint64_t>(x) >> 1)
        );
    }
}

#endif