Incorrect Computation

[yohad]

• unyt version: 2.9.5
• Python version: 3.12.5
• Operating System: Archlinux

Description

Incorrect result of computation:

import unyt as u

T = 10**7 * u.K
print(T**4) # 4477988020393345024 K**4

[neutrinoceros]

Wow, that's weird.
Interestingly, T**2 returns the correct result, but T**3 doesn't.
I also confirm the issue is still present as of unyt 3.0.3.

note: yt.units, albeit a thin wrapper around it, isn't unyt, so I edited the script.

[yohad]

From the following

a = 10**7
a4 = a**4
print(a4 % 2**64) # 4477988020393345024

I think that there must be some casting to 64bits somewhere that destroys the results but I hadn't had the time to find it

[neutrinoceros]

Yes, I arrived at the same conclusion.
Basically what's happening is that the base value 10**7 being an int, unyt_array(10**7, "K") (which is the result of 10**7 * u.K) uses int64 as its dtype. In fact, you can see the same wrong result by doing an explicit int64 computation

>>> import numpy as np
>>> np.int64(10**7)**4
np.int64(4477988020393345024)

This is an inherent limitation from the fact that unyt_array derives from numpy.ndarray and needs a finite-representation datatype, and I don't actually think we can fix that. However, it is actually straightfoward to work around it: use a floating point data type. If you actually want 10**7 as a literal, write it as 1e7 to get a Python float; the result is correct in a float64 computation.