Fix moon orbit generation

[sam52796]

Fixes #5885.

The procedural generation was going wrong when adding moons to planets whose masses were small compared to their star's mass, because calculating the Hill Radius involves a cube root of the mass ratio, and the fixed point algorithm couldn't handle it. The system generator therefore calculated a far too large orbital slice, and filled it with moons at huge orbits. Doing a bit of parameter scaling fixes this issue, but results in different numbers of moons being generated, and hence probably qualifies for a save bump.

The second commit is to fix binary star orbits in custom systems. The procedural generation appears to get binary star systems correct, but if they are added in the system designer, or by hand, parameters can be left to random chance, and then the stars don't line up in their orbits correctly. Binary stars should always be on opposite sides of the barycentre, i.e. at the furthest possible distance allowed by their orbits, such that a line drawn between the two stars passes through the barycentre. Without this fix, 55 Cancri (for example) has both stars on the same side of the barycentre, which is gravitationally impossible.

While the second commit isn't save-breaking in a technical sense, it could lead to players being stranded if the star they are heading towards gets re-oriented to a different point in the orbit, so it seems like a good idea to lump it in with other save-breaking changes.

[sam52796]

I realised that a better fix for the cube root issue was to do the scaling in the cube root function itself, so I've edited the first commit so that it should produce a better result for any caller, instead of relying on the caller to scale where appropriate.

While testing I found another issue which is that moons can be created with 0 kg mass, e.g. the first moon in Gliese 1, so I added a third commit to abort creating a planet/moon if it's going to have a mass of 0. This is probably also save breaking, since now some moons will no longer appear. Now that I think about it, this commit may belong more with my previous "Fix moon populations" pull request, so I can move it there if that's preferred.

[fluffyfreak]

this might be better moved to a common header such as StarSystemGenerator.h which is visible to both, instead of duplicating

[sam52796]

Since there is nothing star system generation specific about the value of pi, I thought it would be better in a more generic include file. It took all of my self-control not to move the definition of pi to the file called Pi.h, but I settled eventually on fixed.h.

[sturnclaw]

As a future-proofing method, I'd request the constant be defined as fixedf::PI instead, as that allows it to have a different definition based on the precision of the fixed-point template specialization instead.

template<...>
class fixedf {
    static const fixedf PI;
};

template<>
constexpr fixed fixed::PI = fixed(..., ...);

template<>
constexpr fixed fixedf<48>::PI = fixedf<48>(...);

[sturnclaw]

I do not have time at current to fully dig into the codepath and analyse the numerical inaccuracy at play here (and the effects of this change on its stability) - however, I'd strongly prefer if the rescaling version of this function be moved into a separate function, e.g. fixedf::CubeRootOfRescale. The automatic rescaling may be undesired behavior in codepaths that operate on values with known input bounds, and I'd prefer it to be opt-in.

That being said, I appreciate your tackling of this issue and I'm looking forward to reviewing it in-depth!

[sam52796]

Are you aware that the cube root function is only used in one other place, and right next to that place there is a comment implying that underflow is a persistent problem for the function? In fact, that comment is what gave me the idea to use the bit shifting approach in the first place, since it explains how to use bit shifting to fix the problem.

I say this because I think it's better not to clutter up the codebase with duplicate library functions where the caller is expected to know which version to call. I realise that sometimes it's unavoidable when you really don't want to risk breaking things, but this doesn't seem to be one of those times.

[sturnclaw]

Are you aware that the cube root function is only used in one other place, and right next to that place there is a comment implying that underflow is a persistent problem for the function?

Yes - I wrote that comment.

I say this because I think it's better not to clutter up the codebase with duplicate library functions where the caller is expected to know which version to call.

I do agree that duplicating functionality is to be avoided. If this function has noticeable loss-of-precision relative to the original, then there is a material difference between the two functions and thus a reason to separate it as a separate entry point.

After having built a test harness for both versions of the function, the proposed implementation has a worst-case error magnitude increase of 20% when the while (a > fixedf(1024, 1)) branch is taken, due to discarding up to 24 bits of fractional precision.

The good news is that replacing that branch with if (a > fixedf(1024, 1)) x = a * fixedf(2, 3); as the initial estimate experimentally fixes a noticed issue where the x * x term can catastrophically overflow (with an input value of 0x4299E58D00000000L), while retaining the same error magnitude as the original version of the function.

When that is resolved, I agree with your intent that the implementation is a direct upgrade in terms of numerical accuracy and stability and retract my initial recommendation.

[sturnclaw]

For clarity, the while (a > fixedf(1024, 1)) branch did not present any reduction in computation error (with regards to a double-precision floating-point cube root using std::pow) compared to computing the N-R method cube root of the original value directly. Additionally, this increase in error metric was observed in both 16.48 fixed-point as well as 32.32 fixed-point representations.

The only noticed outliers as mentioned above were specific values triggering an overflow that caused the first iteration of the method to wildly overshoot beyond the ability of 47 more iterations to converge. Changing the initial estimate of x to a rough estimate of the first iteration result allows large input values to nominally remain within the allowed limits of the fixed-point representation.

[sturnclaw]

After having dug into the behavior of the function across a wide range of input values, I'm happy to find the precision of the result is much improved as long as only a scant few fractional bits of precision are discarded.

I've left suggested changes based on my experimental findings - with all changes applied, the algorithm has effectively uniform accuracy across the entire input range in both 32.32 and 16.48 precision modes (until the value becomes too small for the existing code to convert to double-precision floating-point, at which case I cannot compute a reference value).

This typically represents an order of magnitude (or more) increase in accuracy at high input values (1,000+), and at extremely low values the shift trick is successfully gaining significant amounts of precision.

My original request to split the function was based on my intuition that the shift at high values would discard significant amounts of the fractional bits, and I'm happy to retract that feedback now that a way to avoid the loss of precision has been found.

[sam52796]

I've made all of @sturnclaw's requested changes, I think. C++ templates and fixed point math are both things that I prefer to think about as little as possible, so the suggestions were mainly copied verbatim. It compiles though! And also seems to work in all the cases I tested.

Note that I moved the PI constant template change into its own commit because it felt like its own thang.

Also note that this pull request also fixes #6302.

[sturnclaw]

This looks good! I've left a few nitpick/sanity-check comments, but I consider this ready for squash. I'll merge it along with the next batch of save-bump PRs once we enter the full feature-development cycle for the year; I'd like to cut an interim save-compatible release with our current accumulated bugfixes + the OpenAL audio backend sometime around the end of the month.

[sturnclaw]

Nitpick: if this is duplicated, the file it is duplicated from should be updated to also use the new fixed::PI values.

[sturnclaw]

I believe this can be more succinctly expressed as star_A->GetSuperType() == SystemBodyType::SUPERTYPE_STAR etc.

In general the fixup-at-load approach is good. Ideally an extension to the editor for creating binary (and trinary) orbits from multiple bodies would be nice, though I suspect that would involve non-trivial engineering to track and express that a set of bodies are intended to be in a binary orbit.

[sturnclaw]

Sanity-check: should this be a <= fixedf(0, 1) instead? I believe computing the N-th root of a negative number is invalid by definition until you start involving complex numbers...