docs/blog/04-fixed-point.rst - pigweed/pigweed - Git at Google

 .. _blog-04-fixed-point:

 ========================================================================
 Pigweed Blog #4: Fixed-point arithmetic as a replacement for soft floats
 ========================================================================
 *Published on 2024-09-03 by Leonard Chan*

 Fixed-point arithmetic is an alternative to floating-point arithmetic for
 representing fractional data values (0.5, -1.125, 10.75, etc.). Fixed-point
 types can be represented as :ref:`scaled integers <blog-04-fixed-point-intro>`.
 The advantage here is many
 arithmetic operations (+, -, \*, /, etc.) can be implemented as normal integral
 instructions. This can be useful for embedded systems or other applications
 where floating-point operations can be too expensive or not available.
 Clang has support for fixed-point types according to `ISO TR 18037 <https://standards.iso.org/ittf/PubliclyAvailableStandards/c051126_ISO_IEC_TR_18037_2008.zip>`_.
 I work on a Google project that uses Pigweed.
 Recently, we replaced usage of floats in an embedded project running
 on Cortex M0+. This MCU doesn't provide any hardware support for floating-point
 operations, so floating-point operations are implemented in software. We’ve
 observed a **~2x speed improvement** in
 classification algorithms and a **small decrease in binary size** using fixed-point
 **without sacrificing correctness**.

 All Clang users can benefit from this by adding ``-ffixed-point`` as a
 compilation flag. All of the types and semantics described
 in ISO 18037 are supported.

 Problems with (soft) floats
 ===========================
 One of the ways to represent a decimal number is to use a significand, base,
 and exponent. This is how floating-point numbers according to IEEE 754 are
 represented. The base is two, but the significand and exponent are adjustable.

 When doing arithmetic with floats, many operations are needed to align the
 exponents of the operands and normalize the result. Most modern computers have
 Floating-Point Units (`FPUs <https://en.wikipedia.org/wiki/Floating-point_unit>`_) that do this arithmetic very quickly, but not all
 platforms (especially embedded ones) have this luxury. The alternative is to
 emulate floating-point arithmetic in software which can be costly.

 .. _blog-04-fixed-point-intro:

 ----------------------------------------------
 Intro to fixed-point types and scaled integers
 ----------------------------------------------
 `ISO TR 18037 <https://standards.iso.org/ittf/PubliclyAvailableStandards/c051126_ISO_IEC_TR_18037_2008.zip>`_
 describes 12 fixed-point types with varying scales and ranges:

 * ``unsigned short _Fract``
 * ``unsigned _Fract``
 * ``unsigned long _Fract``
 * ``signed short _Fract``
 * ``signed _Fract``
 * ``signed long _Fract``
 * ``unsigned short _Accum``
 * ``unsigned _Accum``
 * ``unsigned long _Accum``
 * ``signed short _Accum``
 * ``signed _Accum``
 * ``signed long _Accum``

 The scale represents the number of fractional bits in the type. Under the hood,
 Clang implements each of these types as integers. To get the decimal value of
 the type, we just divide the integer by the scale. For example, on x86_64
 Linux, an ``unsigned _Accum`` type is represented as an unsigned 32-bit integer
 with a scale of 16, so 16 bits are used to represent the fractional part and
 the remaining 16 are used to represent the integral part. An ``unsigned _Accum``
 type with a value of ``16.625`` would be represented as ``1089536`` because
 ``1089536 / 2**16 = 16.625``. Additionally, an ``unsigned _Accum`` type has a range
 of ``[0, 65535.99998474121]``, where the max value is represented as ``2^32-1``.
 The smallest possible value we can represent with a fixed-point type is
 ``1/2^scale``.

 The width and scale of each type are platform-dependent. In Clang, different
 targets are free to provide whatever widths
 and scales are best-suited to their needs and capabilities. Likewise,
 LLVM backends are free to lower the different fixed-point intrinsics to native
 fixed-point operations. The default configuration for all targets is the
 “typical desktop processor” implementation described in A.3 of ISO TR 18037.

 Each of these types has a saturating equivalent also (denoted by the ``_Sat``
 keyword). This means operations on a saturating type that would normally cause
 overflow instead clamps to the maximal or minimal value for that type.

 Fixed-point to the rescue
 =========================
 One of the main advantages of fixed-point types is their operations can be done
 with normal integer operations. Addition and subtraction between the same fixed-point
 types normally require just a regular integral ``add`` or ``sub`` instruction.
 Multiplication and division can normally be done with a regular ``mul`` and ``div``
 instruction (plus an extra shift to account for the scale). Platforms that
 don’t support hard floats normally need to make a call to a library function
 that implements the floating-point equivalents of these which can be large or
 CPU-intensive.

 Fixed-point types however are limited in their range and precision which are
 dependent on the number of integral and fractional bits. A ``signed _Accum``
 (which uses 15 fractional bits, 16 integral bits, and 1 sign bit on x86_64
 Linux) will have a range of ``[-65536, 65535.99996948242]`` and a precision of
 ``3.0517578125e-05``. Floats can effectively range from ``(-∞, +∞)`` and have
 precision as small as :math:`10^{-38}`. Maintaining correctness for your program
 largely depends on how much range and precision you intend your fractional
 numbers to have. For addressing ranges, there are sometimes a couple of clever
 math tricks you can do to get large values to fit within these ranges. Some of
 these methods will be discussed later.

 .. list-table:: Fixed-Point vs Floating-Point
    :widths: 10 45 45
    :header-rows: 1
    :stub-columns: 1

    * -
      - Fixed-Point
      - Floating-Point
    * - Range
      - Limited by integral and fractional bits
      - +/-∞
    * - Precision
      - Limited by the fractional bits (:math:`1/2^{scale}`)
      - :math:`10^{-38}`
    * - Binary Operations
      - | Typically fast
        | Usually just involves normal integral binary ops with some shifts
      - | *Hard floats* - Typically fast
        | *Soft floats* - Can be very slow/complex; depends on the implementation
    * - Correctness
      - | Always within 1 ULP of the true result
        | Will always produce the same result for a given op
      - | Always within 0.5 ULP of true result
        | May not always produce the same result for a given op due to rounding errors

 ---------------------
 Running on Cortex-M0+
 ---------------------
 The Arm Cortex-M0+ processor is a common choice for constrained embedded
 applications because it is very energy-efficient. Because it doesn’t have an
 FPU, many floating-point operations targeting this CPU depend on a helper
 `library <https://github.com/ARM-software/abi-aa/blob/7c2fbbd6d32c906c709804f873b67308d1ab46e2/rtabi32/rtabi32.rst#the-standard-compiler-helper-function-library>`_
 to define these as runtime functions. For our application, the code
 running on this processor runs different classification algorithms. One such
 classifier utilizes `softmax regression <http://ufldl.stanford.edu/tutorial/supervised/SoftmaxRegression/>`_.
 The algorithm is roughly:

 .. code-block:: c++

    // This is roughly a trimmed down version of the softmax regression
    // classifier.
    size_t Classify(std::span<const float> sample) const {
        // 1. Compute activations for each category.
        //    All activations are initially zero.
        std::array<float, NumCategories> activations{};
        for (size_t i = 0; i < NumCategories; i++) {
            for (size_t j = 0; j < sample.size(); j++) {
                activations[i] += coefficients_[i][j] * sample[j];
            }
            activations[i] += biases_[i];
        }
        float max_activation = *std::max_element(activations.begin(),
                                                 activations.end());

        // 2. Get e raised to each of these activations.
        std::array<T, NumCategories> exp_activations;
        for (size_t i = 0; i < NumCategories; i++) {
            exp_activations[i] = expf(activations[i]);
        }
        float sum_exp_activations = std::accumulate(exp_activations.begin(),
                                                    exp_activations.end(), 0);

        // 3. Compute each likelihood which us exp(x[i]) / sum(x).
        float max_likelihood = std::numeric_limits<float>::min();
        size_t prediction;
        for (size_t i = 0; i < NumCategories; i++) {
            // 0 <= result.likelihoods[i] < 1
            result.likelihoods[i] = exp_activations[i] / sum_exp_activations;
            if (result.likelihoods[i] > max_likelihood) {
                max_likelihood = result.likelihoods[i];
                prediction = i;  // The highest likelihood is the prediction.
            }
        }

        return prediction;  // Return the index of the highest likelihood.
    }

 Many of these operations involve normal floating-point comparison, addition,
 multiplication, and division. Each of these expands to a call to some
 ``__aeabi_*`` function. This particular function is on a hot path that
 activates, meaning soft float operations take up much computation and power.
 For the normal binary operations, fixed-point types might be a good fit because
 it’s very likely each of the binary operations will result in a value that
 can fit into one of the fixed-point types. (Each of
 the sample values is in the range ``[-256, +256]`` and there are at most 10
 categories. Likewise, each of the ``coefficients_`` and ``biases_`` are small values
 ranging from roughly ``(-3, +3)``.)

 If we wanted to completely replace floats for this function with fixed-point
 types, we’d run into two issues:

 #. There doesn’t exist an ``expf`` function that operates on fixed-point types.
 #. ``expf(x)`` can grow very large for even small for small values of x
    (``expf(23)`` would exceed the largest value possible for an
    ``unsigned long _Accum``).

 Fortunately, we can refactor the code as needed and we can make changes to
 llvm-libc, the libc implementation this device uses.

 llvm-libc and some math
 =======================
 The embedded device's codebase is dependent on a handful of functions that
 take floats: ``expf``, ``sqrtf``, and ``roundf``.
 ``printf`` is also used for printing floats, but support for the fixed-point format
 specifiers is needed. The project already uses llvm-libc, so we can implement
 versions of these functions that operate on fixed-point types.
 The llvm-libc team at Google has been very responsive
 and eager to implement these functions for us! Michael Jones (who implemented
 much of printf in llvm-libc) provided `printf support for each of the fixed
 point types <https://github.com/llvm/llvm-project/pull/82707>`_. Tue Ly
 provided implementations for `abs <https://github.com/llvm/llvm-project/pull/81823>`_,
 `roundf <https://github.com/llvm/llvm-project/pull/81994>`_,
 `sqrtf <https://github.com/llvm/llvm-project/pull/83042>`_,
 `integral sqrt with a fixed-point output <http://github.com/llvm/llvm-project/issues/83924>`_,
 and `expf <https://github.com/llvm/llvm-project/pull/84391>`_ for different
 fixed-point types. Now we have the necessary math functions which accept
 fixed-point types.

 With implementations provided, the next big step is avoiding overflow and
 getting results to fit within the new types. Let’s look at one case involving
 ``expf`` and one involving ``sqrtf``. (Tue Ly also provided these solutions.)

 The softmax algorithm above easily causes overflow with:

 .. code-block::

    exp_activations[i] = expf(activations[i]);

 But the larger goal is to get the likelihood that a sample matches each
 category. This is a value from [0, 1].

 .. math::

    likelihood(i) = \frac{e^{activations[i]}}{\sum_{k=0}^{NumCategories}{e^{activation[k]}}} = \frac{e^{activation[i]]}}{sum(e^{activation[...]}))}

 This can however be normalized by the max activation:

 .. math::

    likelihood(i) = \frac{e^{activation[i]]}}{sum(e^{activation[...]}))} = \frac{e^{activation[i]]}}{sum(e^{activation[...]}))} * \frac{max(e^{activation[...]})}{max(e^{activation[...]})} = \frac{e^{activation[i]-max(activation[...])]}}{sum(e^{activation[...]-max(activation[...])}))}

 This makes the exponent for each component at most zero and the result of
 ``exp(x)`` at most 1 which can definitely fit in the fixed-point types.
 Likewise, the sum of each of these is at most ``NumCategories`` (which happens
 to be 10 for us). For the above code, the only necessary change required is:

 .. code-block::

    exp_activations[i] = expf(activations[i] - max_activation);

 And lucky for us, the precision for a ``signed _Accum`` type is enough to get
 us the same results. One interesting thing is this change alone also improved
 the performance when using floats by 11%. The theory is the ``expf`` implementation
 has a quick switch to see if the exponents need to be adjusted for floats, and
 skips that scaling part when unnecessary.

 The code involving ``sqrtf`` can be adjusted similarly:

 .. code-block:: c++

    void TouchProcessor::CharacteriseRadialDeviation(Touch& touch) {
        // Compute center of touch.
        int32_t sum_x_w = 0, sum_y_w = 0, sum_w = 0;
        // touch.num_position_estimates is at most 100
        for (size_t i = 0; i < touch.num_position_estimates; i++) {
            sum_x_w += touch.position_estimates[i].position.x * 255;
            sum_y_w += touch.position_estimates[i].position.y * 255;
            sum_w += 255;
        }

        // Cast is safe, since average can't exceed any of the individual values.
        touch.center = Point{
            .x = static_cast<int16_t>(sum_x_w / sum_w),
            .y = static_cast<int16_t>(sum_y_w / sum_w),
        };

        // Compute radial deviation from center.
        float sum_d_squared_w = 0.0f;
        for (size_t i = 0; i < touch.num_position_estimates; i++) {
            int32_t dx = touch.position_estimates[i].position.x - touch.center.x;
            int32_t dy = touch.position_estimates[i].position.y - touch.center.y;
            sum_d_squared_w += static_cast<float>(dx * dx + dy * dy) * 255;
        }

        // deviation = sqrt(sum((dx^2 + dy^2) * w) / sum(w))
        touch.features[static_cast<size_t>(Touch::Feature::kRadialDeviation)]
            = sqrtf(sum_d_squared_w / sum_w);
    }

 We know beforehand the maximal values of ``dx`` and ``dy`` are +/-750, so
 ``sum_d_squared_w`` will require as much as 35 integral bits to represent the
 final sum. ``sum_w`` also needs 11 bits, so ``sqrtf`` would need to accept a value
 that can be held in ~24 bits. This can definitely fit into a
 ``signed long _Accum`` which has 32 integral bits, but ideally we’d like to
 use a ``_Sat signed _Accum`` for consistency. Similar to before, we can
 normalize the value by 255. That is, we can avoid multiplying by 255 when
 calculating ``sum_x_w``, ``sum_y_w``, and ``sum_w`` since this will result in
 the exact same ``touch.center`` results. Likewise, if we avoid the ``* 255`` when
 calculating ``sum_d_squared_w``, the final ``sqrtf`` result will be unchanged.
 This makes the code:

 .. code-block:: c++

    void TouchProcessor::CharacteriseRadialDeviation(Touch& touch) {
      // Compute center of touch.
      int32_t sum_x_w = 0, sum_y_w = 0, sum_w = touch.num_position_estimates;
      // touch.num_position_estimates is at most 100
      for (size_t i = 0; i < touch.num_position_estimates; i++) {
          sum_x_w += touch.position_estimates[i].position.x;
          sum_y_w += touch.position_estimates[i].position.y;
      }

      // Cast is safe, since average can't exceed any of the individual values.
      touch.center = Point{
          .x = static_cast<int16_t>(sum_x_w / sum_w),
          .y = static_cast<int16_t>(sum_y_w / sum_w),
      };

      // Compute radial deviation from center.
      float sum_d_squared_w = 0.0f;
      for (size_t i = 0; i < touch.num_position_estimates; i++) {
          int32_t dx = touch.position_estimates[i].position.x - touch.center.x;
          int32_t dy = touch.position_estimates[i].position.y - touch.center.y;
          sum_d_squared_w += static_cast<float>(dx * dx + dy * dy);
      }

      // deviation = sqrt(sum((dx^2 + dy^2) * w) / sum(w))
      touch.features[static_cast<size_t>(Touch::Feature::kRadialDeviation)]
          = sqrtf(sum_d_squared_w / sum_w);
    }

 This allows ``sum_d_squared_w`` to fit within 31 integral bits. We can go
 further and realize ``sum_d_squared_w`` will always have an integral value and
 replace its type with a ``uint32_t``. This will mean any fractional part from
 ``sum_d_squared_w / sum_w`` will be discarded, but for our use case, this is
 acceptable since the integral part of ``sum_d_squared_w / sum_w`` is so large
 that the fractional part is negligible. ``touch.features[i]`` also isn’t
 propagated significantly so we don’t need to worry about propagation of
 error. With this, we can replace the ``sqrtf`` with one that accepts an integral
 type and returns a fixed-point type:

 .. code-block:: c++

    touch.features[static_cast<size_t>(Touch::Feature::kRadialDeviation)]
        = isqrt(sum_d_squared_w / sum_w);

 -------
 Results
 -------

 Classification runtime performance
 ==================================
 Before any of this effort, a single run of classifying sensor input took on
 average **856.8 us**.

 Currently, with all floats substituted with fixed-point types, a single run is
 **412.845673 us** which is a **~2.1x speed improvement**.

 It’s worth noting that the optimizations we described above also improved
 performance when using floats, making a single run using floats
 **771.475464 us** which translates to **~1.9x speed improvement**.

 Size
 ====
 We’ve observed **~1.25 KiB** saved (out of a **~177 KiB** binary) when
 switching to fixed-point.

 .. code-block::

    $ bloaty app.elf.fixed -- app.elf.float
        FILE SIZE        VM SIZE
     --------------  --------------
      +0.5% +3.06Ki  [ = ]       0    .debug_line
      +0.1% +1.89Ki  [ = ]       0    .debug_str
      +0.2%    +496  [ = ]       0    .debug_abbrev
      +0.6%    +272  +0.6%    +272    .bss
      -0.3%     -16  -0.3%     -16    .data
      -0.3%    -232  [ = ]       0    .debug_frame
      -1.3%    -712  [ = ]       0    .debug_ranges
      -0.3% -1.11Ki  [ = ]       0    .debug_loc
      -1.2% -1.50Ki  -1.2% -1.50Ki    .text
      -1.6% -1.61Ki  [ = ]       0    .symtab
      -3.2% -3.06Ki  [ = ]       0    .pw_tokenizer.entries
      -1.2% -4.00Ki  [ = ]       0    .strtab
      -0.4% -19.6Ki  [ = ]       0    .debug_info
      -0.3% -26.1Ki  -0.7% -1.25Ki    TOTAL

 A large amount of this can be attributed to not needing many of the
 ``__aeabi_*`` float functions provided by libc. All instances of floats were
 removed so they don’t get linked into the final binary. Instead, we use the
 fixed-point equivalents provided by llvm-libc. Much of these fixed-point
 functions are also smaller and make fewer calls to other functions. For
 example, the ``sqrtf`` we used makes calls to ``__ieee754_sqrtf``,
 ``__aeabi_fcmpun``, ``__aeabi_fcmplt``, ``__errno``, and ``__aeabi_fdiv``
 whereas ``isqrtf_fast`` only makes calls to ``__clzsi2``, ``__aeabi_lmul``, and
 ``__aeabi_uidiv``. It’s worth noting that individual fixed-point binary
 operations are slightly larger since they involve a bunch of shifts or
 saturation checks. Floats tend to make a call to the same library function, but
 fixed-point types emit instructions at every callsite.

 .. code-block::

    float_add:
      push    {r7, lr}
      add     r7, sp, #0
      bl      __aeabi_fadd
      pop     {r7, pc}

    fixed_add:
      adds    r0, r0, r1
      bvc     .LBB1_2
      asrs    r1, r0, #31
      movs    r0, #1
      lsls    r0, r0, #31
      eors    r0, r1
    .LBB1_2:
      bx      lr

 Although the callsite is shorter, it’s worth noting that the soft float
 functions can’t be inlined because they tend to be very large. The
 ``__aeabi_fadd`` in particular is **444 bytes large**.

 Correctness
 ===========
 We observe the exact same results for our classification algorithms when
 using fixed-point types. We have **over 15000 classification tests** with none
 of them producing differing results. This effectively means we get all the
 mentioned benefits without any correctness cost.

 -----------------
 Using fixed-point
 -----------------
 Fixed-point arithmetic can be used out-of-the-box with Clang by adding
 ``-ffixed-point`` on the command line. All of the types and semantics described
 in ISO 18037 are supported. ``llvm-libc`` is the only ``libc`` that supports
 fixed-point printing and math functions at the moment. Not all libc functions
 described in ISO 18037 are supported at the moment, but they will be added
 eventually! Pigweed users can use these fixed-point functions by using
 the ``//pw_libc:stdfix`` target.

 -----------
 Future work
 -----------
 * An option we could take to reduce size (at the cost of even more performance)
   is to potentially teach LLVM to outline fixed-point additions. This could result
   in something similar to the float libcalls.
 * See if we can instead use a regular ``signed _Accum`` in the codebase instead
   of a ``_Sat signed _Accum`` to save some instructions (or even see if we can
   use a smaller type).
 * Investigate performance compared to hard floats. If we’re able to afford
   using an unsaturated type, then many native floating-point ops could be
   replaced with normal integral ops.
 * Full libc support for the fixed-point C functions.
	.. _blog-04-fixed-point:

	========================================================================
	Pigweed Blog #4: Fixed-point arithmetic as a replacement for soft floats
	========================================================================
	Published on 2024-09-03 by Leonard Chan

	Fixed-point arithmetic is an alternative to floating-point arithmetic for
	representing fractional data values (0.5, -1.125, 10.75, etc.). Fixed-point
	types can be represented as :ref:`scaled integers <blog-04-fixed-point-intro>`.
	The advantage here is many
	arithmetic operations (+, -, \*, /, etc.) can be implemented as normal integral
	instructions. This can be useful for embedded systems or other applications
	where floating-point operations can be too expensive or not available.
	Clang has support for fixed-point types according to `ISO TR 18037 <https://standards.iso.org/ittf/PubliclyAvailableStandards/c051126_ISO_IEC_TR_18037_2008.zip>`_.
	I work on a Google project that uses Pigweed.
	Recently, we replaced usage of floats in an embedded project running
	on Cortex M0+. This MCU doesn't provide any hardware support for floating-point
	operations, so floating-point operations are implemented in software. We’ve
	observed a ~2x speed improvement in
	classification algorithms and a small decrease in binary size using fixed-point
	without sacrificing correctness.

	All Clang users can benefit from this by adding ``-ffixed-point`` as a
	compilation flag. All of the types and semantics described
	in ISO 18037 are supported.

	Problems with (soft) floats
	===========================
	One of the ways to represent a decimal number is to use a significand, base,
	and exponent. This is how floating-point numbers according to IEEE 754 are
	represented. The base is two, but the significand and exponent are adjustable.

	When doing arithmetic with floats, many operations are needed to align the
	exponents of the operands and normalize the result. Most modern computers have
	Floating-Point Units (`FPUs <https://en.wikipedia.org/wiki/Floating-point_unit>`_) that do this arithmetic very quickly, but not all
	platforms (especially embedded ones) have this luxury. The alternative is to
	emulate floating-point arithmetic in software which can be costly.

	.. _blog-04-fixed-point-intro:

	----------------------------------------------
	Intro to fixed-point types and scaled integers
	----------------------------------------------
	`ISO TR 18037 <https://standards.iso.org/ittf/PubliclyAvailableStandards/c051126_ISO_IEC_TR_18037_2008.zip>`_
	describes 12 fixed-point types with varying scales and ranges:

	* ``unsigned short _Fract``
	* ``unsigned _Fract``
	* ``unsigned long _Fract``
	* ``signed short _Fract``
	* ``signed _Fract``
	* ``signed long _Fract``
	* ``unsigned short _Accum``
	* ``unsigned _Accum``
	* ``unsigned long _Accum``
	* ``signed short _Accum``
	* ``signed _Accum``
	* ``signed long _Accum``

	The scale represents the number of fractional bits in the type. Under the hood,
	Clang implements each of these types as integers. To get the decimal value of
	the type, we just divide the integer by the scale. For example, on x86_64
	Linux, an ``unsigned _Accum`` type is represented as an unsigned 32-bit integer
	with a scale of 16, so 16 bits are used to represent the fractional part and
	the remaining 16 are used to represent the integral part. An ``unsigned _Accum``
	type with a value of ``16.625`` would be represented as ``1089536`` because
	``1089536 / 2**16 = 16.625``. Additionally, an ``unsigned _Accum`` type has a range
	of ``[0, 65535.99998474121]``, where the max value is represented as ``2^32-1``.
	The smallest possible value we can represent with a fixed-point type is
	``1/2^scale``.

	The width and scale of each type are platform-dependent. In Clang, different
	targets are free to provide whatever widths
	and scales are best-suited to their needs and capabilities. Likewise,
	LLVM backends are free to lower the different fixed-point intrinsics to native
	fixed-point operations. The default configuration for all targets is the
	“typical desktop processor” implementation described in A.3 of ISO TR 18037.

	Each of these types has a saturating equivalent also (denoted by the ``_Sat``
	keyword). This means operations on a saturating type that would normally cause
	overflow instead clamps to the maximal or minimal value for that type.

	Fixed-point to the rescue
	=========================
	One of the main advantages of fixed-point types is their operations can be done
	with normal integer operations. Addition and subtraction between the same fixed-point
	types normally require just a regular integral ``add`` or ``sub`` instruction.
	Multiplication and division can normally be done with a regular ``mul`` and ``div``
	instruction (plus an extra shift to account for the scale). Platforms that
	don’t support hard floats normally need to make a call to a library function
	that implements the floating-point equivalents of these which can be large or
	CPU-intensive.

	Fixed-point types however are limited in their range and precision which are
	dependent on the number of integral and fractional bits. A ``signed _Accum``
	(which uses 15 fractional bits, 16 integral bits, and 1 sign bit on x86_64
	Linux) will have a range of ``[-65536, 65535.99996948242]`` and a precision of
	``3.0517578125e-05``. Floats can effectively range from ``(-∞, +∞)`` and have
	precision as small as :math:`10^{-38}`. Maintaining correctness for your program
	largely depends on how much range and precision you intend your fractional
	numbers to have. For addressing ranges, there are sometimes a couple of clever
	math tricks you can do to get large values to fit within these ranges. Some of
	these methods will be discussed later.

	.. list-table:: Fixed-Point vs Floating-Point
	:widths: 10 45 45
	:header-rows: 1
	:stub-columns: 1

	* -
	- Fixed-Point
	- Floating-Point
	* - Range
	- Limited by integral and fractional bits
	- +/-∞
	* - Precision
	- Limited by the fractional bits (:math:`1/2^{scale}`)
	- :math:`10^{-38}`
	* - Binary Operations
	- \| Typically fast
	\| Usually just involves normal integral binary ops with some shifts
	- \| Hard floats - Typically fast
	\| Soft floats - Can be very slow/complex; depends on the implementation
	* - Correctness
	- \| Always within 1 ULP of the true result
	\| Will always produce the same result for a given op
	- \| Always within 0.5 ULP of true result
	\| May not always produce the same result for a given op due to rounding errors

	---------------------
	Running on Cortex-M0+
	---------------------
	The Arm Cortex-M0+ processor is a common choice for constrained embedded
	applications because it is very energy-efficient. Because it doesn’t have an
	FPU, many floating-point operations targeting this CPU depend on a helper
	`library <https://github.com/ARM-software/abi-aa/blob/7c2fbbd6d32c906c709804f873b67308d1ab46e2/rtabi32/rtabi32.rst#the-standard-compiler-helper-function-library>`_
	to define these as runtime functions. For our application, the code
	running on this processor runs different classification algorithms. One such
	classifier utilizes `softmax regression <http://ufldl.stanford.edu/tutorial/supervised/SoftmaxRegression/>`_.
	The algorithm is roughly:

	.. code-block:: c++

	// This is roughly a trimmed down version of the softmax regression
	// classifier.
	size_t Classify(std::span<const float> sample) const {
	// 1. Compute activations for each category.
	// All activations are initially zero.
	std::array<float, NumCategories> activations{};
	for (size_t i = 0; i < NumCategories; i++) {
	for (size_t j = 0; j < sample.size(); j++) {
	activations[i] += coefficients_[i][j] * sample[j];
	}
	activations[i] += biases_[i];
	}
	float max_activation = *std::max_element(activations.begin(),
	activations.end());

	// 2. Get e raised to each of these activations.
	std::array<T, NumCategories> exp_activations;
	for (size_t i = 0; i < NumCategories; i++) {
	exp_activations[i] = expf(activations[i]);
	}
	float sum_exp_activations = std::accumulate(exp_activations.begin(),
	exp_activations.end(), 0);

	// 3. Compute each likelihood which us exp(x[i]) / sum(x).
	float max_likelihood = std::numeric_limits<float>::min();
	size_t prediction;
	for (size_t i = 0; i < NumCategories; i++) {
	// 0 <= result.likelihoods[i] < 1
	result.likelihoods[i] = exp_activations[i] / sum_exp_activations;
	if (result.likelihoods[i] > max_likelihood) {
	max_likelihood = result.likelihoods[i];
	prediction = i; // The highest likelihood is the prediction.
	}
	}

	return prediction; // Return the index of the highest likelihood.
	}

	Many of these operations involve normal floating-point comparison, addition,
	multiplication, and division. Each of these expands to a call to some
	``__aeabi_*`` function. This particular function is on a hot path that
	activates, meaning soft float operations take up much computation and power.
	For the normal binary operations, fixed-point types might be a good fit because
	it’s very likely each of the binary operations will result in a value that
	can fit into one of the fixed-point types. (Each of
	the sample values is in the range ``[-256, +256]`` and there are at most 10
	categories. Likewise, each of the ``coefficients_`` and ``biases_`` are small values
	ranging from roughly ``(-3, +3)``.)

	If we wanted to completely replace floats for this function with fixed-point
	types, we’d run into two issues:

	#. There doesn’t exist an ``expf`` function that operates on fixed-point types.
	#. ``expf(x)`` can grow very large for even small for small values of x
	(``expf(23)`` would exceed the largest value possible for an
	``unsigned long _Accum``).

	Fortunately, we can refactor the code as needed and we can make changes to
	llvm-libc, the libc implementation this device uses.

	llvm-libc and some math
	=======================
	The embedded device's codebase is dependent on a handful of functions that
	take floats: ``expf``, ``sqrtf``, and ``roundf``.
	``printf`` is also used for printing floats, but support for the fixed-point format
	specifiers is needed. The project already uses llvm-libc, so we can implement
	versions of these functions that operate on fixed-point types.
	The llvm-libc team at Google has been very responsive
	and eager to implement these functions for us! Michael Jones (who implemented
	much of printf in llvm-libc) provided `printf support for each of the fixed
	point types <https://github.com/llvm/llvm-project/pull/82707>`_. Tue Ly
	provided implementations for `abs <https://github.com/llvm/llvm-project/pull/81823>`_,
	`roundf <https://github.com/llvm/llvm-project/pull/81994>`_,
	`sqrtf <https://github.com/llvm/llvm-project/pull/83042>`_,
	`integral sqrt with a fixed-point output <http://github.com/llvm/llvm-project/issues/83924>`_,
	and `expf <https://github.com/llvm/llvm-project/pull/84391>`_ for different
	fixed-point types. Now we have the necessary math functions which accept
	fixed-point types.

	With implementations provided, the next big step is avoiding overflow and
	getting results to fit within the new types. Let’s look at one case involving
	``expf`` and one involving ``sqrtf``. (Tue Ly also provided these solutions.)

	The softmax algorithm above easily causes overflow with:

	.. code-block::

	exp_activations[i] = expf(activations[i]);

	But the larger goal is to get the likelihood that a sample matches each
	category. This is a value from [0, 1].

	.. math::

	likelihood(i) = \frac{e^{activations[i]}}{\sum_{k=0}^{NumCategories}{e^{activation[k]}}} = \frac{e^{activation[i]]}}{sum(e^{activation[...]}))}

	This can however be normalized by the max activation:

	.. math::

	likelihood(i) = \frac{e^{activation[i]]}}{sum(e^{activation[...]}))} = \frac{e^{activation[i]]}}{sum(e^{activation[...]}))} * \frac{max(e^{activation[...]})}{max(e^{activation[...]})} = \frac{e^{activation[i]-max(activation[...])]}}{sum(e^{activation[...]-max(activation[...])}))}

	This makes the exponent for each component at most zero and the result of
	``exp(x)`` at most 1 which can definitely fit in the fixed-point types.
	Likewise, the sum of each of these is at most ``NumCategories`` (which happens
	to be 10 for us). For the above code, the only necessary change required is:

	.. code-block::

	exp_activations[i] = expf(activations[i] - max_activation);

	And lucky for us, the precision for a ``signed _Accum`` type is enough to get
	us the same results. One interesting thing is this change alone also improved
	the performance when using floats by 11%. The theory is the ``expf`` implementation
	has a quick switch to see if the exponents need to be adjusted for floats, and
	skips that scaling part when unnecessary.

	The code involving ``sqrtf`` can be adjusted similarly:

	.. code-block:: c++

	void TouchProcessor::CharacteriseRadialDeviation(Touch& touch) {
	// Compute center of touch.
	int32_t sum_x_w = 0, sum_y_w = 0, sum_w = 0;
	// touch.num_position_estimates is at most 100
	for (size_t i = 0; i < touch.num_position_estimates; i++) {
	sum_x_w += touch.position_estimates[i].position.x * 255;
	sum_y_w += touch.position_estimates[i].position.y * 255;
	sum_w += 255;
	}

	// Cast is safe, since average can't exceed any of the individual values.
	touch.center = Point{
	.x = static_cast<int16_t>(sum_x_w / sum_w),
	.y = static_cast<int16_t>(sum_y_w / sum_w),
	};

	// Compute radial deviation from center.
	float sum_d_squared_w = 0.0f;
	for (size_t i = 0; i < touch.num_position_estimates; i++) {
	int32_t dx = touch.position_estimates[i].position.x - touch.center.x;
	int32_t dy = touch.position_estimates[i].position.y - touch.center.y;
	sum_d_squared_w += static_cast<float>(dx * dx + dy * dy) * 255;
	}

	// deviation = sqrt(sum((dx^2 + dy^2) * w) / sum(w))
	touch.features[static_cast<size_t>(Touch::Feature::kRadialDeviation)]
	= sqrtf(sum_d_squared_w / sum_w);
	}

	We know beforehand the maximal values of ``dx`` and ``dy`` are +/-750, so
	``sum_d_squared_w`` will require as much as 35 integral bits to represent the
	final sum. ``sum_w`` also needs 11 bits, so ``sqrtf`` would need to accept a value
	that can be held in ~24 bits. This can definitely fit into a
	``signed long _Accum`` which has 32 integral bits, but ideally we’d like to
	use a ``_Sat signed _Accum`` for consistency. Similar to before, we can
	normalize the value by 255. That is, we can avoid multiplying by 255 when
	calculating ``sum_x_w``, ``sum_y_w``, and ``sum_w`` since this will result in
	the exact same ``touch.center`` results. Likewise, if we avoid the ``* 255`` when
	calculating ``sum_d_squared_w``, the final ``sqrtf`` result will be unchanged.
	This makes the code:

	.. code-block:: c++

	void TouchProcessor::CharacteriseRadialDeviation(Touch& touch) {
	// Compute center of touch.
	int32_t sum_x_w = 0, sum_y_w = 0, sum_w = touch.num_position_estimates;
	// touch.num_position_estimates is at most 100
	for (size_t i = 0; i < touch.num_position_estimates; i++) {
	sum_x_w += touch.position_estimates[i].position.x;
	sum_y_w += touch.position_estimates[i].position.y;
	}

	// Cast is safe, since average can't exceed any of the individual values.
	touch.center = Point{
	.x = static_cast<int16_t>(sum_x_w / sum_w),
	.y = static_cast<int16_t>(sum_y_w / sum_w),
	};

	// Compute radial deviation from center.
	float sum_d_squared_w = 0.0f;
	for (size_t i = 0; i < touch.num_position_estimates; i++) {
	int32_t dx = touch.position_estimates[i].position.x - touch.center.x;
	int32_t dy = touch.position_estimates[i].position.y - touch.center.y;
	sum_d_squared_w += static_cast<float>(dx * dx + dy * dy);
	}

	// deviation = sqrt(sum((dx^2 + dy^2) * w) / sum(w))
	touch.features[static_cast<size_t>(Touch::Feature::kRadialDeviation)]
	= sqrtf(sum_d_squared_w / sum_w);
	}

	This allows ``sum_d_squared_w`` to fit within 31 integral bits. We can go
	further and realize ``sum_d_squared_w`` will always have an integral value and
	replace its type with a ``uint32_t``. This will mean any fractional part from
	``sum_d_squared_w / sum_w`` will be discarded, but for our use case, this is
	acceptable since the integral part of ``sum_d_squared_w / sum_w`` is so large
	that the fractional part is negligible. ``touch.features[i]`` also isn’t
	propagated significantly so we don’t need to worry about propagation of
	error. With this, we can replace the ``sqrtf`` with one that accepts an integral
	type and returns a fixed-point type:

	.. code-block:: c++

	touch.features[static_cast<size_t>(Touch::Feature::kRadialDeviation)]
	= isqrt(sum_d_squared_w / sum_w);

	-------
	Results
	-------

	Classification runtime performance
	==================================
	Before any of this effort, a single run of classifying sensor input took on
	average 856.8 us.

	Currently, with all floats substituted with fixed-point types, a single run is
	412.845673 us which is a ~2.1x speed improvement.

	It’s worth noting that the optimizations we described above also improved
	performance when using floats, making a single run using floats
	771.475464 us which translates to ~1.9x speed improvement.

	Size
	====
	We’ve observed ~1.25 KiB saved (out of a ~177 KiB binary) when
	switching to fixed-point.

	.. code-block::

	$ bloaty app.elf.fixed -- app.elf.float
	FILE SIZE VM SIZE
	-------------- --------------
	+0.5% +3.06Ki [ = ] 0 .debug_line
	+0.1% +1.89Ki [ = ] 0 .debug_str
	+0.2% +496 [ = ] 0 .debug_abbrev
	+0.6% +272 +0.6% +272 .bss
	-0.3% -16 -0.3% -16 .data
	-0.3% -232 [ = ] 0 .debug_frame
	-1.3% -712 [ = ] 0 .debug_ranges
	-0.3% -1.11Ki [ = ] 0 .debug_loc
	-1.2% -1.50Ki -1.2% -1.50Ki .text
	-1.6% -1.61Ki [ = ] 0 .symtab
	-3.2% -3.06Ki [ = ] 0 .pw_tokenizer.entries
	-1.2% -4.00Ki [ = ] 0 .strtab
	-0.4% -19.6Ki [ = ] 0 .debug_info
	-0.3% -26.1Ki -0.7% -1.25Ki TOTAL

	A large amount of this can be attributed to not needing many of the
	``__aeabi_*`` float functions provided by libc. All instances of floats were
	removed so they don’t get linked into the final binary. Instead, we use the
	fixed-point equivalents provided by llvm-libc. Much of these fixed-point
	functions are also smaller and make fewer calls to other functions. For
	example, the ``sqrtf`` we used makes calls to ``__ieee754_sqrtf``,
	``__aeabi_fcmpun``, ``__aeabi_fcmplt``, ``__errno``, and ``__aeabi_fdiv``
	whereas ``isqrtf_fast`` only makes calls to ``__clzsi2``, ``__aeabi_lmul``, and
	``__aeabi_uidiv``. It’s worth noting that individual fixed-point binary
	operations are slightly larger since they involve a bunch of shifts or
	saturation checks. Floats tend to make a call to the same library function, but
	fixed-point types emit instructions at every callsite.

	.. code-block::

	float_add:
	push {r7, lr}
	add r7, sp, #0
	bl __aeabi_fadd
	pop {r7, pc}

	fixed_add:
	adds r0, r0, r1
	bvc .LBB1_2
	asrs r1, r0, #31
	movs r0, #1
	lsls r0, r0, #31
	eors r0, r1
	.LBB1_2:
	bx lr

	Although the callsite is shorter, it’s worth noting that the soft float
	functions can’t be inlined because they tend to be very large. The
	``__aeabi_fadd`` in particular is 444 bytes large.

	Correctness
	===========
	We observe the exact same results for our classification algorithms when
	using fixed-point types. We have over 15000 classification tests with none
	of them producing differing results. This effectively means we get all the
	mentioned benefits without any correctness cost.

	-----------------
	Using fixed-point
	-----------------
	Fixed-point arithmetic can be used out-of-the-box with Clang by adding
	``-ffixed-point`` on the command line. All of the types and semantics described
	in ISO 18037 are supported. ``llvm-libc`` is the only ``libc`` that supports
	fixed-point printing and math functions at the moment. Not all libc functions
	described in ISO 18037 are supported at the moment, but they will be added
	eventually! Pigweed users can use these fixed-point functions by using
	the ``//pw_libc:stdfix`` target.

	-----------
	Future work
	-----------
	* An option we could take to reduce size (at the cost of even more performance)
	is to potentially teach LLVM to outline fixed-point additions. This could result
	in something similar to the float libcalls.
	* See if we can instead use a regular ``signed _Accum`` in the codebase instead
	of a ``_Sat signed _Accum`` to save some instructions (or even see if we can
	use a smaller type).
	* Investigate performance compared to hard floats. If we’re able to afford
	using an unsaturated type, then many native floating-point ops could be
	replaced with normal integral ops.
	* Full libc support for the fixed-point C functions.