np.float128 on x86-64 is 80-bit extended precision, not quad
NumPy 2.5.4 / SciPy 1.18.1 on x86-64 Linux. np.float128 exists and takes 128 bits of storage, but finfo(np.longdouble) shows nmant=63 and eps=1.08e-19. That is x87 extended precision, about 3 more decimal digits than float64 (eps 2.2e-16), nowhere near IEEE quad (eps 1.9e-34): 1 + 2**-64 == 1 in it. Three traps. np.longdouble(0.1) widens the double 0.1 and gains nothing; it equals 0.1, while np.longdouble('0.1') (a string) is closer. np.linalg rejects it: inv raises TypeError 'array type float128 is unsupported in linalg'. scipy.special is patchy: expit and logsumexp return float128, but erf and gamma raise TypeError for float128 input. To check a float64 result, use mpmath or fractions, or rearrange the algorithm; longdouble is not a portable reference precision.
Every point reproduced on x86-64 Linux, NumPy 2.5.4 / SciPy 1.18.1 / CPython 3.14.10. finfo(float128): nmant=63, eps=1.084202172485504434e-19, itemsize 16, and np.longdouble is np.float128. 1 + 2**-63 > 1 is True and 1 + 2**-64 > 1 is False. np.longdouble(0.1) == 0.1 is True, while np.longdouble('0.1') is closer to 1/10 computed in longdouble. np.linalg.inv on a float128 eye(2) raises TypeError 'array type float128 is unsupported in linalg'. scipy.special: expit and logsumexp return float128, erf and gamma raise "ufunc ... not supported for the input types".