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

Fix for a corner-case numerical problem and corresponding unit test. #181

Open
wants to merge 3 commits into
base: master
Choose a base branch
from

Conversation

EdFuentetaja
Copy link

@EdFuentetaja EdFuentetaja commented Jan 16, 2025

Thank you for this excellent library.

I noticed this problem decoding some ADS-B messages. The local airborne decoding failed is some particular cases.

There seems to be some oddities with the Python modulus operator. E.g.:

>>> 30.508474576271183 % 6.101694915254237
6.101694915254235
>>> import math
>>> def mod(x, y):
...     return x - y * math.floor(x / y)
...
>>> mod(30.508474576271183, 6.101694915254237)
0.0

Checking the standard (1090 MOPS, Vol.1 DO-260C), section 1.7.5 in Appendix A, shows an equivalent and simpler expression that doesn't require use of a modulus function. I think this should be the preferred approach since it's more accurate and requires less computations.

lat_ref = 30.508474576271183 # Close to (360.0/59.0)*5
lon_ref = 7.2*5.0+3e-15
pos = adsb.airborne_position_with_ref(
"8D06A15358BF17FF7D4A84", lat_ref, lon_ref
Copy link
Collaborator

Choose a reason for hiding this comment

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

The correct message here (with CRC) would be 8d06a15358bf17ff7d4a84b47b95

Copy link
Author

Choose a reason for hiding this comment

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

Thanks for this. I added the CRC to the message.

@xoolive
Copy link
Collaborator

xoolive commented Jan 16, 2025

Thanks for pointing that out!
I have a related question: are there any required changes in the decoding of surface messages as well?

@EdFuentetaja
Copy link
Author

Thanks for pointing that out! I have a related question: are there any required changes in the decoding of surface messages as well?

You are right, it's the same situation. I just applied the same fix there.

In general, given the odd behavior of the % operator, I would replace it everywhere by its definition:

def mod(x,y):
    return x - y*math.floor(x / y)

I'd like to know more about how the % operator is implemented in Python to understand where the numerical error is coming from.

@EdFuentetaja
Copy link
Author

EdFuentetaja commented Jan 17, 2025

Python's % operator is mostly a call to the fmod C function. On x86 with -ffast-math this is implemented with a call to the fprem CPU instruction, which is affected by the same numerical problem. I think this is a very peculiar corner case.
Numerically accurate computation is hard. I'd follow the recommendation of the standard and implement MOD as stated in A.1.7.2.c:

MOD(x, y) = x - y floor(x/y) where y ≠ 0.

@xoolive
Copy link
Collaborator

xoolive commented Jan 17, 2025

Yes nice comment. We should look into more occurrences of calculations of j and m in this file.

FYI I did the same in my Rust version following your comment (current version here but opened issue and PR accordingly) and it seems the modulo in Rust behaves differently but has a similar issue.

@junzis
Copy link
Owner

junzis commented Jan 17, 2025

Thank you, @EdFuentetaja. Awesome bug fix! we will fix the surface position functions too.

@EdFuentetaja
Copy link
Author

Thank you @junzis and @xoolive. Is there anything else I could do you help this PR get merged? Thank you for the positive feedback.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

4 participants