Improve accuracy for atan2(y,x) with near-zero x

[grimaldini]

Hi,

I have the following setup:

#include <fpm/math.hpp>
using number = fpm::fixed<int64_t, __int128_t, 32>;
number x = number(50) * (number::pi() / number(100));
// x is 1.5708, really close to half pi but raw value is off by one (base type) from ::half_pi()
number cosX = cos(x);
// inside the cos->sin function we reduce x to a value of two but it is also off by one (base type) and hence returns the wrong result here.

Could this be caused by this:
static constexpr fixed pi() { return from_fixed_point<61>(7244019458077122842ll); }
I tried changing this to calculate fixed(355)/fixed(113) (approximation of pi) and use this instead and it seems to fix the issue but could the from_fixed_point function not be correctly converting between fixed point types?

Thanks,
Luis

[MikeLankamp]

Hi @grimaldini, since fixed-point math uses integers under the hood, off-by-one errors in the raw values are to be expected depending on how certain operations are performed. You must be mindful of loss of information.

pi() / 100 * 50 loses a few least significant bits. pi() / 2 does not. Let's run through this (all examples for your number type):

pi() has a raw value of 13493037705. Note that the "true" value should be π * 2^32 = 13493037704.5220...
half_pi() has a raw value of 6746518852.
pi() / 100 has a raw value of 134930377 (notice how we lost the 5 at the end)
pi() / 100 * 50 has a raw value of 6746518850. That's a difference of 2 with pi() / 2. That's exactly because of the 5 that we lost. It's also a difference of 0.0000000004 on 1.5708.

Note that the approximation of 355/133 is actually worse:
number(355) has a raw value of 1524713390080.
number(113) has a raw value of 485331304448.
number(355) / number(113) has a raw value of 13493038850 instead of 13493037704.
Half of that is 6746519425 instead of 6746518852. That's a difference of 573 instead of 2.

I'm not sure what you mean by "cos->sin reducing x to a value of two", but any off-by-one errors there are hardly suprising considering it's fixed-point math.

In conclusion, it seems everything's working as it should. half_pi() is correct. pi() / 100 * 50 is slightly off because of loss of information in intermediate steps. 355/133 is a bad approximation, as expected.

Please let me know if this answers your question.

[grimaldini]

that makes sense but these trigonometric functions don't seem to be able to handle such low results, for example:

    number x = number::from_raw_value(1); // simulating getting a loss of information case here
    number y = atan2(-number_one, x);
    cout << "x: " << x << ", y: " << y << endl; // x: 2.3283e-10, y: 0

y should be -pi or -1.5708 but instead it results in 0

[MikeLankamp]

y should be -pi or -1.5708 but instead it results in 0

Allright, that's a fair point. Apparently, atan2 doesn't work well if its inputs are nearly zero.

[MikeLankamp]

So after a quick check, atan2(-1, 0 + epsilon) has an acceptable error when using a 20.12 fixed point format (-1.5706 vs -1.5708), but not when using a 16.16 format (0 vs -1.5708) for this input.

I do think this is also a known consequence of using fixed-point math. atan2(y, x) is implemented as atan(y/x) for non-zero x. If x is very very small, then the division could easily overflow and result in 0. Assigning more bits to the integer part than the fraction part is a way to solve this, and asin, acos, atan and atan2 all rely on this.

It might be possible to have an implementation of atan2 that doesn't require this, but I'm seeing that as an improvement rather than a bug.

@grimaldini can you work around your problem by converting to a different fixed-point format before doing atan2?

[grimaldini]

I'm currently using 32.32 format though, so I do have more bits in the integer part. Do you mean that the integer part needs to be larger than the fractional part, if so why is that the case? is this a problem that might happen with other trig functions, does checking that all the trig functions of (0 + epsilon) can guarantee that a similar issue won't occur?

[MikeLankamp]

Do you mean that the integer part needs to be larger than the fractional part

Yes, indeed. At least, if you want atan2 to work with a very very small x. Because atan2 does atan(y/x) when x is non-zero, if x is the smallest it can be, then y/x is very, very large. Larger than what fits in an equal number of integer bits.

Example: when dividing 1 (the smallest non-zero integer) by 1*2^-N (the smallest non-zero fraction in a number with N fraction bits), the result is 2^N. This doesn't fit in an N-bit integer. So a quick workaround is having more integer bits than fraction bits.

Fixed-point math is not quite a drop-in replacement, since you have to care about the domain and range of functions.

That said, atan2 could try to detect this and approximate the answer.

[grimaldini]

Are you saying that for a number with N fraction bits, the integer part needs to be at least N+1?

Separately, I tried to use Q17.15 but now I am getting other errors such as

using number = fpm::fixed<std::int32_t, std::int64_t, 15>;
number z = number::from_raw_value(16384000);
number w = z * z + z * z;
cout << "w: " << w << std::endl; // w: -0.185303

How can fixed point arithmetics be reliably used in a complex system such as a physics simulation? where I cannot

[MikeLankamp]

Are you saying that for a number with N fraction bits, the integer part needs to be at least N+1?

That depends on what inputs you want to have work. If you want division by very small (near-zero) numbers to work, then yes, at least N+1. Likely more.

atan2 is special in that it makes an exception for x = 0 but not for x = 0.000001, hence the sudden discrepancy. I can look to mitigate this case for atan2, but that doesn't fix the other "problems" you've mentioned.

but now I am getting other errors such as

number::from_raw_value(16384000) is the same as number(500) for Q17.15. 2 * 500^2 is 500000 which requires 19 bits to store. So that doesn't fit in an integer of 17 bits.

How can fixed point arithmetics be reliably used in a complex system such as a physics simulation?

Fixed-point numbers are deterministic, which is their benefit. But this comes at the cost of having to care about underflow and overflow. You can't just pick a type and work with it as if It were a float. It is fundamental property of fixed point that you have to care about the range of operations. You may have to convert between different fixed-point types during a calculation.

[grimaldini]

Thanks for the fix! do you think this near-zero input issue might be something affecting other trigonometric functions?

[MikeLankamp]

This particular improvement also affected asin(x) (when x is near -1 or 1) and acos(x) (when x is near -1), so you might see some improvement there.