Normalizing a vector with great accuracy and big numbers using fixed-point arithmetic

Question

Why do I need this?

I'm creating a game about space. For a space game to work it'll need big and accurate numbers. Floating-point numbers are not really great for such application because if you get pretty far away in the world the precision will be worse, so the physics wont be the same, and so on.

What's the problem

In space there generally tends to be planets. So to create one I have to generate the spheres (planets) mesh. The problem is if I want a planet like Jupiter (with the radius of ~69911km), when I normalize a point and multiply it by the radius of the planet, and therefore "put" it on the surface of the planet, I don't have enough precision (the mesh isn't exactly round, there's an error of about 10-15m). Here's some footage of that error (the flickering of the planet is due to video compression, and the side length of the smallest square in the mesh grid is about ~10m). The problem isn't in the numbers, it's in the method I use to normalize the point.

In the video I am using 64-bit fixed-point numbers with 20 bits of precision. And I'm normalizing the points just by multiplying them with the fast inverse square root (I'm using just one iteration of the Newtons method, with two its better but still not nearly enough). In which I convert the fixed-point numbers to double precision floats just for the fast inverse square root.

My thoughts and attempts

I think the only solution to get such accurate vector normalization is to use an iterative method:

1. Calculate the point normally (normalize it and multiply by the planets radius)
2. Nudge it closer and closer to the radius until the error is small enough (I would like for it to be 0.01m)

The 2. step is the hard one. I don't have the math skills or any experience in this field of computer science to figure it out. We can iteratively increase or decrease a vector by some amount, but what amount, and how do we know we wont overshoot it. I thought about trying the Newtons method again but on the actual coordinates, not the inverse square root of the length. Which includes 128-bit division to preserve the precision needed, which cant be done efficiently (remember i have to do this maybe a million times, for each vertex of the planet mesh)

Any thoughts on this?

And the scale doesn't have to stop there, what if i wanted to make a star, the radius of the sun is 10 times bigger than that of Jupiter. And there are bigger stars.

I'm probably not the first one to ponder this question since many people have tried making the earth at a real scale (remember that earths radius is 10 time smaller than Jupiters). And i probably wont be the last one to try this so there will be an answer sooner or later.

Update 1, after some comments

1. Yes im using a 3D vector with 64bit numbers with 20 bits of precision as numbers. And I am obviously using a local coordinate system to the planet for creating the points (I need to for the normalization to work as needed).
2. Using floating-point numbers is my last hope. They're really not made for this kind of application because of the precision difference at different scales.
3. The scales I want to make this work for the celestial bodies (stars, planets) is at a minimum the size of the Sun. For the position of objects I am using two 64-bit ints (effectively one 128-bit int) so the world in the game can be bigger than our actual observable universe if I want to make it (thanks to fixed-point numbers). But I don't need 128-bit ints for the positions of vertices of planets, 64-bit is good enough (though I can use 128-bit ints in the process to calculate the values). And as I mentioned, it would be best for the mesh vertices to have an error less than 0.01 meters, where 1 unit of a coordinate is 1 meter (imagine walking on the surface of the planet, anything bigger than 0.01m could be noticeable).

The reason I'm pushing fixed-point numbers so much is because I can't use the position of objects as floating-point numbers because of physics. And lets say we have a big planet, bigger than Jupiter. Because I need to calculate the vertices relative to Jupiters position and then subtract the cameras position from them in the shaders (which are mostly limited to 32-bit floats), the errors would just add up and it would be noticable. There are workarounds but the core problem will always come up with using floats.

TL;DR

I need a way to calculate a point on a sphere with a radius that can be as big as 1000000000 = 10^9, but with an error less than 0.01, given the angles at which the point needs to rest. And this method needs to be relatively efficient.

Answer 1

My understanding is that you are

1. Making a 3D vector, represented in units of meters as 64-bit fixed-point numbers with 20 bits of precision.
2. Normalizing the vector to unit length, scaling by a fast inverse square root approximation.
3. Scaling the vector by Jupiter's radius 69911 km.

The challenge is that any error from inv sqrt approximation or round-off in the normalized vector is magnified when scaling by the planet radius. To achieve 0.01 m accuracy, the unit vector's length needs to accurate with error less than

0.01 m / 69911 km = 1.43×10^-10.

To get there, the inv sqrt needs to be accurate to 11 digits or better, and the unit vector needs 32 fractional bits at a minimum (for round-off error of 2^-33 = 1.16×10^-10).

I suggest computing the surface points in a local coordinate system with origin at the center of the planet represented with 64-bit floats, then translating the points to the universe coordinate system. 64-bit floats are accurate to ~16 digits, enough to get the desired accuracy. This seems like the easiest and most efficient solution assuming 64-bit floats are an option.

Alternatively for an iterative "nudging" approach, you can do this:

R = planet radius
x, y, z = surface point with correct angle but possibly wrong length
for 1, 2, ...
  scale = 0.5 * (R² / (x² + y² + z²) + 1)
  x *= scale
  y *= scale
  z *= scale

This is a fixed point iteration that converges rapidly. An example run, starting with a very inaccurate surface point x, y, z = [6991.1, 55928.8, -6991.1]:

#  x        y          z        Radius
-------------------------------------------------
0  6991.10  55928.80  -6991.10  56795.96489065047
1  8791.84  70334.70  -8791.84  71425.22857460588
2  8607.42  68859.40  -8607.42  69927.05096841766
3  8605.45  68843.60  -8605.45  69911.00184215968

A further thought:

remember i have to do this maybe a million times, for each vertex of the planet mesh

It looks like from the video that you are already applying a level of detail scheme to reduce the number of vertices when far away. The hierarchical mesh can be leveraged to reduce vertices when near the surface as well: skip finer-stage vertices when their parents are off screen based on frustum culling.