Arcsine function for a real-time control program

Question

I have been developing control software for a three-phase induction motor. The main task of this C++ code is to control the torque of the motor. In one part of the algorithm (The control algorithm is executed with a period of 100 µs.), I need to calculate the arcsine function. Here is my implementation. I would appreciate your assessment.

Math.h

#include <cstdint>

class Math {

public:

/**
 * @brief Function calculates arcsine of a given argument.
 * @param x argument, \f$x\in\left<-1, 1\right>\f$
 * @return arcsine value of x
 */
static float arcsine(float x);

private:
static const uint16_t kArcSineLutSize = 101;
static const float arcsine_lut[kArcSineLutSize];

}

Math.cpp

float Math::arcsine(float x)
{
  bool neg_table_value = false;
  if (x < 0) {
    // arcsin(-x) = -arcsin(x)
    x = -x;
    neg_table_value = true;
  }
  uint16_t index = static_cast<uint16_t>(x * 100.0);
  float tmp =
      (arcsine_lut[index + 1] - arcsine_lut[index]) * (100.0 * x - index) +
      arcsine_lut[index];
  if (neg_table_value) {
    return -tmp;
  } else {
    return tmp;
  }
}

const float Math::arcsine_lut[Math::kArcSineLutSize] = {
    0.0000, 0.0100, 0.0200, 0.0300, 0.0400, 0.0500, 0.0600, 0.0701, 0.0801,
    0.0901, 0.1002, 0.1102, 0.1203, 0.1304, 0.1405, 0.1506, 0.1607, 0.1708,
    0.1810, 0.1912, 0.2014, 0.2116, 0.2218, 0.2321, 0.2424, 0.2527, 0.2630,
    0.2734, 0.2838, 0.2942, 0.3047, 0.3152, 0.3257, 0.3363, 0.3469, 0.3576,
    0.3683, 0.3790, 0.3898, 0.4006, 0.4115, 0.4225, 0.4334, 0.4445, 0.4556,
    0.4668, 0.4780, 0.4893, 0.5007, 0.5121, 0.5236, 0.5352, 0.5469, 0.5586,
    0.5704, 0.5824, 0.5944, 0.6065, 0.6187, 0.6311, 0.6435, 0.6561, 0.6687,
    0.6816, 0.6945, 0.7076, 0.7208, 0.7342, 0.7478, 0.7615, 0.7754, 0.7895,
    0.8038, 0.8183, 0.8331, 0.8481, 0.8633, 0.8788, 0.8947, 0.9108, 0.9273,
    0.9442, 0.9614, 0.9791, 0.9973, 1.0160, 1.0353, 1.0552, 1.0759, 1.0973,
    1.1198, 1.1433, 1.1681, 1.1944, 1.2226, 1.2532, 1.2870, 1.3252, 1.3705,
    1.4293, 1.5708};

Answer 1

Use a namespace instead of a class

There is no reason to use a class if you only have static member functions and variables. Consider creating a namespace instead:

namespace Math {
    static float arcsine(float x);
}

You don't even have to declare the private static members in the header file this way.

What are the requirements?

Why are there only 100 entries in the table? Why is every value only specified up to 4 decimals? This seems very arbitrary. The fewer points you have in your table and the less precise the points in the table are, the bigger the error is compared to the real arcsine. So you should ask yourself: what is the error budget?

A quick check of your functions versus std::asin() reveals that near 0, you have errors of about 0.06%, which is probably fine, but near ±1 you have an error of almost 2.4%. If that is not good enough for your purpose, you probably need to increase the number of points, use more precision for each point, and/or use something a bit more sophisticated than linear interpolation.

I also have to wonder why you are using a look-up table instead of just using std::asin(). If your CPU is not fast enough to do this calculation every 100 µs, then it's a good approach. If it's fast enough but it doesn't come with a standard library that contains an arcsine function, then you could consider implementing your own. I would then recommend using Chebyshev polynomials.

Calculate the table at compile-time

Instead of hardcoding the values in the look-up table, you can have the compiler generate the table for you. For example:

template<std::size_t N>
consteval auto generate_arcsine_lut() {
    std::array<float, N> values;
    for (std::size_t i = 0; i < N; ++i) {
        values[i] = std::asin(i / static_cast<float>(N));
    }
    return values;
}

static const auto arcsine_lut = generate_arcsine_lut<100>();

See this on godbolt.org. Note however that std::asin() is not a constexpr function, and while GCC is fine with this, Clang does not want to compile it. The way around this is to implement your own implementation of asin() (as mentioned above).

Do you really need the arcsine?

The main task of this C++ software is to control the torque of the motor. In a part of the control algorithm (the control algorithm is executed with 100μs period.) I need to calculate the arcsine function.

You might have some mathematical formulas with an angle in them, but in the end you are not interested in that angle, just in how much you have to energize each coil of the motor. So probably you have some values coming in (from Hall sensors or some rotary encoder?), from which you create a vector of the current direction the axis is pointing in. If you just want to rotate that vector a little bit, then instead of calculating the angle, adding something to the angle, and converting that back into a vector, you can just multiply the first vector with a rotation matrix. The matrix is constant, so once precalculated you can just apply it to a vector using a few multiplication and additions.

If you want to vary the rotation you apply a lot, but if the rotation angle is always very small, then you can use the small-angle approximations of sin and cos to calculate that matrix very cheaply. Of course, keep your error budget in mind.

Answer 2

@G. Sliepen offers a good point: "Do you really need the arcsine?" as it is quite possible a simpler function will suffice.

Yet let us assume we want to get the best arcsine(x) via a small look-up table.

• Bug: Table too small. When x == 1.0, code attempts to access the table outside the range.

// 1.4293, 1.5708};
1.4293, 1.5708, /* something */ 1.4293 };

• Bug: Table does not well handle fabs(x) > 1.0

• To reduce spurious warnings, use float constants for a float table rather than double ones.

• Using a fixed 4-5 digit look-up table is slightly inferior to using a floating 9 digit one (for float). A 9 significant decimal digit table adds zero overhead and is more accurate. OP's 2.358% worst case error vs. potential 2.357%.

const float arcsine_lut[] = {
    0.0f, 0.0100001665f, 0.0200013332f, 0.0300045013f, 0.0400106721f, 0.0500208586f, 0.0600360557f,
    0.0700572953f, 0.0800855756f, 0.0901219472f, 0.100167423f, 0.110223047f, 0.120289877f, 0.130368978f,
    0.140461415f, 0.150568277f, 0.16069065f, 0.170829669f, 0.180986464f, 0.191162139f, 0.201357931f,
    0.211574957f, 0.221814469f, 0.232077688f, 0.242365852f, 0.252680242f, 0.263022184f, 0.273393035f,
    0.283794105f, 0.294226825f, 0.304692656f, 0.315193027f, 0.325729489f, 0.336303592f, 0.346916914f,
    0.357571095f, 0.368267894f, 0.379009038f, 0.389796287f, 0.400631577f, 0.411516845f, 0.422454059f,
    0.433445305f, 0.444492787f, 0.455598682f, 0.466765314f, 0.477995217f, 0.489290774f, 0.500654697f,
    0.512089789f, 0.52359879f, 0.535184741f, 0.54685092f, 0.558600545f, 0.570437133f, 0.582364261f,
    0.594385803f, 0.606505871f, 0.618728638f, 0.631058812f, 0.643501163f, 0.656060576f, 0.668742716f,
    0.681553185f, 0.694498241f, 0.707584381f, 0.720818818f, 0.734208822f, 0.74776262f, 0.761489034f,
    0.77539748f, 0.78949815f, 0.803802371f, 0.818322003f, 0.833070397f, 0.848062098f, 0.863313079f,
    0.878841102f, 0.894665778f, 0.91080904f, 0.927295208f, 0.944152117f, 0.961410999f, 0.979107618f,
    0.997283161f, 1.01598537f, 1.03526974f, 1.05520236f, 1.07586217f, 1.09734511f, 1.11976945f,
    1.14328408f, 1.16808057f, 1.19441283f, 1.22263026f, 1.25323582f, 1.28700209f, 1.32523096f,
    1.37046158f, 1.42925692f, 1.57079637f, 1.42925692f};

• Pedantic: use if (signbit(x)) instead of if (x < 0) to 1) return a -0.0 for arcsine(-0.0) and perhaps a tad faster.

• Using fixed points at x = 0.00, 0.01, 0.02, ... 1.00 is not the most optimal solution. TBD on that. (We likely can reduce the error by ~50%, but it takes some time/thought to form that table.)

• Avoid mixing double math in a float problem.

// uint16_t index = static_cast<uint16_t>(x * 100.0);
uint16_t index = static_cast<uint16_t>(x * 100.0f);

// float tmp = (arcsine_lut[index + 1] - arcsine_lut[index]) 
//       * (100.0 * x - index) + arcsine_lut[index];
float tmp = (arcsine_lut[index + 1] - arcsine_lut[index])
      * (100.0f * x - index) + arcsine_lut[index];

---

Follow up

• Rather than scale by 100.0, consider a power of 2, like 64 or 128. It makes for slightly more accurate answers and is sometimes easier to optimize.

• The biggest error comes from |x| > 0.99 as that is where the arcsine function becomes (nearly) vertical. I was able to improve the worst case by 8x by breaking the look-up-table into 2 groups: [0.0 ... 0.9885] and [0.9885... 1.0], each having the same size table. Other table distribution schemes that have more entries in higher x than low x (unlike the linear table OP is using) result in more accurate answers or smaller table size for the same benefit.