Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Add safe acos and safe asin #76876

Closed
wants to merge 1 commit into from
Closed
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
3 changes: 1 addition & 2 deletions core/math/basis.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -807,8 +807,7 @@ void Basis::get_axis_angle(Vector3 &r_axis, real_t &r_angle) const {
z = (rows[1][0] - rows[0][1]) / s;

r_axis = Vector3(x, y, z);
// CLAMP to avoid NaN if the value passed to acos is not in [0,1].
r_angle = Math::acos(CLAMP((rows[0][0] + rows[1][1] + rows[2][2] - 1) / 2, (real_t)0.0, (real_t)1.0));
r_angle = Math::safe_acos((rows[0][0] + rows[1][1] + rows[2][2] - 1) / 2);
}

void Basis::set_quaternion(const Quaternion &p_quaternion) {
Expand Down
24 changes: 20 additions & 4 deletions core/math/math_funcs.h
Original file line number Diff line number Diff line change
Expand Up @@ -74,11 +74,27 @@ class Math {
static _ALWAYS_INLINE_ double tanh(double p_x) { return ::tanh(p_x); }
static _ALWAYS_INLINE_ float tanh(float p_x) { return ::tanhf(p_x); }

static _ALWAYS_INLINE_ double asin(double p_x) { return ::asin(p_x); }
static _ALWAYS_INLINE_ float asin(float p_x) { return ::asinf(p_x); }
// GUIDELINES on use of asin / acos.
// Use safe version if you are SURE a clamp is needed.
// Use unsafe version if you are SURE a clamp is not needed, OR Nan is dealt with (and likely to be a bottleneck).
// Use generic version if you are not sure.
static _ALWAYS_INLINE_ double safe_asin(double p_x) { return p_x < -1 ? (-Math_PI / 2) : (p_x > 1 ? (Math_PI / 2) : ::asin(p_x)); }
static _ALWAYS_INLINE_ float safe_asin(float p_x) { return p_x < -1 ? (-Math_PI / 2) : (p_x > 1 ? (Math_PI / 2) : ::asinf(p_x)); }

static _ALWAYS_INLINE_ double acos(double p_x) { return ::acos(p_x); }
static _ALWAYS_INLINE_ float acos(float p_x) { return ::acosf(p_x); }
static _ALWAYS_INLINE_ double unsafe_asin(double p_x) { return ::asin(p_x); }
static _ALWAYS_INLINE_ float unsafe_asin(float p_x) { return ::asinf(p_x); }

static _ALWAYS_INLINE_ double asin(double p_x) { return safe_asin(p_x); }
static _ALWAYS_INLINE_ float asin(float p_x) { return safe_asin(p_x); }
Comment on lines +87 to +88
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I understand the intent to have both an explicit safe_asin and the main asin which is safe by default, but having two synonymous methods just means that we'll introduce inconsistency across the codebase where some contributors use asin and others use safe_asin.

Copy link
Member Author

@lawnjelly lawnjelly May 10, 2023

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This is a tricky one. I added the explicit safe_acos to replace parts of code that were already clamping. This signifies areas where the author / bugfix has shown clamping to be absolutely necessary.

In general there are three cases rather than two:

  1. You absolutely know you must clamp (safe_acos)
  2. You aren't sure, and not a bottleneck, so might as well (acos)
  3. You have determined the Nan condition can't happen and / or this is a bottleneck and Nans are dealt with (unsafe_acos)

We could alternatively use a comment to specify that a clamp is necessary instead of safe_acos (1), but having all three is clear and easily searchable.

I can adjust the PR to use comments to distinguish (1) and (2) but imo I would caution against it. What happens when you later come to optimize a function that uses acos? You now no longer know whether it must be safe or not (unless you re-analyse the code, which may be complex for math) - has someone missed a comment etc. There is essentially a danger of e.g. 10 different programmers examine the same bit of code in the following years and waste time trying to work out the intent again and again, when the code just could have been explicit from day 1.

I can add some comment explanation to the math_funcs.h file to make it clearer when to use each version, as a guide for contributors. 🤔

The other alternative which I may have mentioned earlier, is to disallow the generic acos and force all (core) use to be explicit. This is a compromise but you do lose a little information, but might be viable. We can't easily remove the generic to prevent compat breaking in modules etc.

Anyway I am happy to modify the PR if broad opinion favours one of the other options, but just want to caution that "simplifying" potentially has a cost in terms of future maintenance / bug prevention. 🙂

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think this might benefit from some deeper analysis of the uses of the trigonometric functions, from the case in #8111 there might be reasons to use atan2 instead in situations where it is possible, i.e. where a well defined x and y can be found, I think fixing the immediate problem of NaNs might be sufficient to merge this with a placeholder function and then investigate individual cases to follow up

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

(Not sure about real code) But at least in micro benchmarks the unsafe version is only about 0.5% faster then the branched version (if all inputs are actually in the -1 - 1 range, otherwise the branched version quickly wins), so I am not sure if the unsafe version of these functions is actually necessary at all.

----------------------------------------------------------
Benchmark                Time             CPU   Iterations
----------------------------------------------------------
BenchAsinClamp       30815 ns        30815 ns        22743
BenchAsinBranch      27037 ns        27037 ns        25850
BenchAsinUnsafe      26934 ns        26934 ns        26773

Copy link
Member Author

@lawnjelly lawnjelly May 10, 2023

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@RedworkDE that's actually a good point. The unsafe version by definition shouldn't be hitting the -1 / +1 limit. I guess we could do this. The other option is that the unsafe version is so rare (not even used as yet) that we could directly call ::acos at the call site with a comment. 🤔

(This is difficult to benchmark because the benchmark with that data will presumably be using perfect branch prediction, and in the real world it is more likely to be mixed instructions.)

The other option is to just use safe version everywhere, and delay this discussion until we actually find a case where profiling shows would benefit from the unsafe(!). 😄

I'll create another simpler PR with just the safe version and will let you guys decide which to go for.

Copy link
Member Author

@lawnjelly lawnjelly May 10, 2023

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@RedworkDE comparing unsafe to the branching, I'm getting unsafe being around:
4% faster (when some are outside range)
10% faster (when all inside range)

This probably depends very much on the benchmark / CPU / compiler / compiler options etc. But in general the clamped version shouldn't be that much slower (and in practice, often cache is the bottleneck rather than these instructions).

Also, as ever, benchmarking is notoriously easy to get wrong 😄 , yours may be more reliable if you are using google benchmark, I'm doing this in a test app.


static _ALWAYS_INLINE_ double safe_acos(double p_x) { return p_x < -1 ? Math_PI : (p_x > 1 ? 0 : ::acos(p_x)); }
static _ALWAYS_INLINE_ float safe_acos(float p_x) { return p_x < -1 ? Math_PI : (p_x > 1 ? 0 : ::acosf(p_x)); }

static _ALWAYS_INLINE_ double unsafe_acos(double p_x) { return ::acos(p_x); }
static _ALWAYS_INLINE_ float unsafe_acos(float p_x) { return ::acosf(p_x); }

static _ALWAYS_INLINE_ double acos(double p_x) { return safe_acos(p_x); }
static _ALWAYS_INLINE_ float acos(float p_x) { return safe_acos(p_x); }

static _ALWAYS_INLINE_ double atan(double p_x) { return ::atan(p_x); }
static _ALWAYS_INLINE_ float atan(float p_x) { return ::atanf(p_x); }
Expand Down
2 changes: 1 addition & 1 deletion core/math/quaternion.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -35,7 +35,7 @@

real_t Quaternion::angle_to(const Quaternion &p_to) const {
real_t d = dot(p_to);
return Math::acos(CLAMP(d * d * 2 - 1, -1, 1));
return Math::safe_acos(d * d * 2 - 1);
}

Vector3 Quaternion::get_euler(EulerOrder p_order) const {
Expand Down
4 changes: 2 additions & 2 deletions doc/classes/@GlobalScope.xml
Original file line number Diff line number Diff line change
Expand Up @@ -65,7 +65,7 @@
<return type="float" />
<param index="0" name="x" type="float" />
<description>
Returns the arc cosine of [param x] in radians. Use to get the angle of cosine [param x]. [param x] must be between [code]-1.0[/code] and [code]1.0[/code] (inclusive), otherwise, [method acos] will return [constant @GDScript.NAN].
Returns the arc cosine of [param x] in radians. Use to get the angle of cosine [param x]. [param x] will be clamped between [code]-1.0[/code] and [code]1.0[/code] (inclusive), in order to prevent [method acos] from returning [constant @GDScript.NAN].
[codeblock]
# c is 0.523599 or 30 degrees if converted with rad_to_deg(c)
var c = acos(0.866025)
Expand All @@ -76,7 +76,7 @@
<return type="float" />
<param index="0" name="x" type="float" />
<description>
Returns the arc sine of [param x] in radians. Use to get the angle of sine [param x]. [param x] must be between [code]-1.0[/code] and [code]1.0[/code] (inclusive), otherwise, [method asin] will return [constant @GDScript.NAN].
Returns the arc sine of [param x] in radians. Use to get the angle of sine [param x]. [param x] will be clamped between [code]-1.0[/code] and [code]1.0[/code] (inclusive), in order to prevent [method asin] from returning [constant @GDScript.NAN].
[codeblock]
# s is 0.523599 or 30 degrees if converted with rad_to_deg(s)
var s = asin(0.5)
Expand Down
10 changes: 6 additions & 4 deletions tests/core/math/test_math_funcs.h
Original file line number Diff line number Diff line change
Expand Up @@ -157,15 +157,17 @@ TEST_CASE_TEMPLATE("[Math] asin/acos/atan", T, float, double) {
CHECK(Math::asin((T)0.1) == doctest::Approx((T)0.1001674212));
CHECK(Math::asin((T)0.5) == doctest::Approx((T)0.5235987756));
CHECK(Math::asin((T)1.0) == doctest::Approx((T)1.5707963268));
CHECK(Math::is_nan(Math::asin((T)1.5)));
CHECK(Math::is_nan(Math::asin((T)450.0)));
CHECK(Math::asin((T)2.0) == doctest::Approx((T)1.5707963268));
CHECK(Math::is_nan(Math::unsafe_asin((T)1.5)));
CHECK(Math::is_nan(Math::unsafe_asin((T)450.0)));

CHECK(Math::acos((T)-0.1) == doctest::Approx((T)1.670963748));
CHECK(Math::acos((T)0.1) == doctest::Approx((T)1.4706289056));
CHECK(Math::acos((T)0.5) == doctest::Approx((T)1.0471975512));
CHECK(Math::acos((T)1.0) == doctest::Approx((T)0.0));
CHECK(Math::is_nan(Math::acos((T)1.5)));
CHECK(Math::is_nan(Math::acos((T)450.0)));
CHECK(Math::acos((T)2.0) == doctest::Approx((T)0.0));
CHECK(Math::is_nan(Math::unsafe_acos((T)1.5)));
CHECK(Math::is_nan(Math::unsafe_acos((T)450.0)));

CHECK(Math::atan((T)-0.1) == doctest::Approx((T)-0.0996686525));
CHECK(Math::atan((T)0.1) == doctest::Approx((T)0.0996686525));
Expand Down