Different Result for Float vs. Array

[JoshPaterson]

When I calculate an altitude using a float for the jd, and compare it to an altitude calculated when that same float is put into an array, I get a different result.

The difference is only in the last digit, which in itself is not causing me any problems, but I don't know if it's a symptom of a bigger problem.

When I run this:

from skyfield.api import Loader, Topos
from numpy import array

data_directory = ''
load = Loader(data_directory)
ts = load.timescale()

ephem = load('de430t.bsp')
sun = ephem['sun']
earth = ephem['earth']

greenwich = earth + Topos('55 N', '0 W')

jd_float = 2458275.9999190294
jd_array = array([jd_float])

float_result = greenwich.at(ts.tt(jd=jd_float)).observe(sun).apparent().altaz()[0].degrees
print(float_result)

array_result = greenwich.at(ts.tt(jd=jd_array)).observe(sun).apparent().altaz()[0].degrees[0]
print(array_result)

I get this:

57.67041953992528
57.67041953992529

[brandon-rhodes]

Differences in the last digit of a floating point calculation are to be expected and, happily, are not a problem. Even very slight differences in how a computation is performed — on some processors, even the choice of which CPU register they get stored in — will change at least the last digit, or even the last two digits, of a floating point operation.

On my ThinkPad, for example, I just ran your script and got

57.670419539925305
57.670419539925305

The difference matters only if the two numbers are of actually physically different significance. When looking at a rounding difference, always compute the physical meaning:

57.670419539925305 - 57.67041953992528 = 2.8421709e-14

2.8421709e-14 degrees * 3600 = 1.0231815e-10 arcseconds

If I'm doing the math correctly here, the difference is around 1e-10 arcseconds. Since the current most precise radiotelescopes can only distinguish down to 1e-4 arcseconds — and visible light telescopes are far less sensitive — the rounding differences between our machines, and between Python floats and NumPy floating-point arrays that you observed, are not of physical significance.

The same thing is not quite true of Skyfield dates! You can read this if you're curious about how we do run into some limitations there:

http://aa.usno.navy.mil/software/novas/USNOAA-TN2011-02.pdf

[brandon-rhodes]

Oh — and, if you want to read more about floating-point numbers themselves and the crazy little details of rounding, I recently ran across a professor who seems to be rather famous in that area:

http://people.eecs.berkeley.edu/~wkahan/

In particular search that page for the link to his paper "How Futile are Mindless Assessments of Roundoff in Floating-Point Computation?" which, if I'm remembering correctly, was recently cited by a Skyfield contributor advocating for a different algorithm in one of the Skyfield routines involving angles! It's a quite interesting corner of the computational numerics field.

[JoshPaterson]

Haha, I’m actually the contributor who cited that! Looks like I need to give it another read through ☺️ Thanks for the detailed response!

[brandon-rhodes]

I kind of wondered if you had been involved in that citation, but I figured it would be faster for you to tell me if you had than for me to go looking for it since I didn't remember if it was in an issue or a pull request :)