are numpys dtypes kept even for non-numpy multiplications/divisions/sums?

Question

I am working on a simple simulator for a radio frequency application and have to deal with very low complex numbers. During the process I have a Matrix like np.array([[A,B],[C,D]], dtype=np.clongdouble) which ensures the necessary "resolution(?)". However, I have to do stuff like

den = A+B/z0+C*z0+D
s11 = A+B/z0-C*z0-D)/den
s12 = 2*(A*D-B*C)/den
s21 = 2/den
s22 = (-A+B/z0-C*z0+D)/den

I think Z0 is of type double since it's calculated without numpy.

Now I wonder: Do I have to do the calculations of den etc with numpy to achieve/keep resolution or are the 'normal' calculations sufficient enough?

Python's floating-point numbers are usually 64-bit floating-point numbers, nearly equivalent to np. — D.L, Aug 16 '22 at 17:40
Terminology: I think what you mean by "resolution" is what we call "precision" — Code-Apprentice, Aug 16 '22 at 17:46

score 1 · Answer 1 · edited Aug 18 '22 at 12:49

1

Python’s floating-point numbers are usually 64-bit floating-point numbers, nearly equivalent to np.float64.

we can actually see this from the basic types in the docs here: https://numpy.org/doc/stable/user/basics.types.html

So Python alone is satisfactory for accuracy.

edited Aug 18 '22 at 12:49

norok2

25,683
4
73
99

answered Aug 16 '22 at 17:44

D.L

4,339
5
22
45

2

Nearly equivalent? I'd say exactly equivalent, unless provided with evidence to the contrary. – Mark Ransom Aug 16 '22 at 17:51
"So Python alone is satisfactory for acuracy." why do you say so? The dtype of the matrices involved in the calculation use twice the number of bits of Python's floats. – norok2 Aug 17 '22 at 08:33
@norok2, it is **from the numpy docs** (right at the bottom) which i have provided a link to. – D.L Aug 17 '22 at 17:35
in the docs I read clongdouble to be two extended-float which is 128 bit. – norok2 Aug 17 '22 at 17:41
@MarkRansom, you may well be right... although, there must have be a reason that the word "nearly" was specifically used. – D.L Aug 17 '22 at 17:46
Probably just so nobody would file a bug if they noticed a difference. – Mark Ransom Aug 17 '22 at 19:16
1

@MarkRansom: The non-equivalence they're referring to is most likely in things like edge case behavior, methods supported, and broadcasting behavior, not in floating point representation format. For example, `1.0/0.0` will raise a ZeroDivisionError, while `numpy.float64(1.0)/numpy.float64(0.0)` will issue a warning and return an infinity value under default settings. – user2357112 Aug 18 '22 at 12:38

score 1 · Accepted Answer · edited Aug 18 '22 at 17:42

Python's complex numbers are represented "as a pair of machine-level double precision floating point numbers". This is typically 2 x 64 = 128 bits.

NumPy's clongdouble are represented as a pair of machine-level "extended-precision floating-point numbers". This is typically 2 x 80 = 160 bits.

However, what NumPy provides for long double (as per @WarrenWeckesser's comment) is platform and compiler dependent.

Therefore, the datatype you are using from NumPy may (or may not) have higher precision.

To force a number to have a specific precision, you can just cast it with np.array(..., dtype=...):

import numpy as np


a = 1.0j
b = np.array(a, dtype=np.clongdouble)
print(b.dtype, b.nbytes)
# complex256 32

Note that the names here are a bit misleading as the number of bits refer to the alignment and not necessarily to the number of bits actually used for the data.

However, you do not need to do any specific casting for the above formulae to work, as the result of an operation with a higher precision and a lower precision operand would result in a higher precision result:

c = np.array(a)
print(c.dtype, c.nbytes)
# complex128 16

d = b + c
print(d.dtype, d.nbytes)
# complex256 32

e = b + a
print(e.dtype, e.nbytes)
# complex256 32

print(d == e)
# True

Eventually, you may want to compute z0 with NumPy's higher precision data types, if that precision is available and it is critical to your result.

Note that this table would be your starting point to see how NumPy maps Python data types.

Note on Floating Point information

Unfortunately, you cannot easily get the information on the number of bits for the different data types that Python uses.

You can get some information with sys.getsizeof() but this will have information on the memory usage of the Python object, which includes Python's book-keeping information.

However, you can get some information on Python's float with sys.float_info:

import sys


a = 10+1.0j
print(sys.float_info)
# sys.float_info(max=1.7976931348623157e+308, max_exp=1024, max_10_exp=308, min=2.2250738585072014e-308, min_exp=-1021, min_10_exp=-307, dig=15, mant_dig=53, epsilon=2.220446049250313e-16, radix=2, rounds=1)

Similarly, you can get NumPy floating-point information with np.finfo() (note that when feeding complex datatypes, you get the information on the underlying floats):

for dtype in (float, np.float_, np.double, np.cdouble, np.longdouble, np.clongdouble):
    print(np.finfo(dtype))

Machine parameters for float64
---------------------------------------------------------------
precision =  15   resolution = 1.0000000000000001e-15
machep =    -52   eps =        2.2204460492503131e-16
negep =     -53   epsneg =     1.1102230246251565e-16
minexp =  -1022   tiny =       2.2250738585072014e-308
maxexp =   1024   max =        1.7976931348623157e+308
nexp =       11   min =        -max
---------------------------------------------------------------

Machine parameters for float64
---------------------------------------------------------------
precision =  15   resolution = 1.0000000000000001e-15
machep =    -52   eps =        2.2204460492503131e-16
negep =     -53   epsneg =     1.1102230246251565e-16
minexp =  -1022   tiny =       2.2250738585072014e-308
maxexp =   1024   max =        1.7976931348623157e+308
nexp =       11   min =        -max
---------------------------------------------------------------

Machine parameters for float64
---------------------------------------------------------------
precision =  15   resolution = 1.0000000000000001e-15
machep =    -52   eps =        2.2204460492503131e-16
negep =     -53   epsneg =     1.1102230246251565e-16
minexp =  -1022   tiny =       2.2250738585072014e-308
maxexp =   1024   max =        1.7976931348623157e+308
nexp =       11   min =        -max
---------------------------------------------------------------

Machine parameters for float64
---------------------------------------------------------------
precision =  15   resolution = 1.0000000000000001e-15
machep =    -52   eps =        2.2204460492503131e-16
negep =     -53   epsneg =     1.1102230246251565e-16
minexp =  -1022   tiny =       2.2250738585072014e-308
maxexp =   1024   max =        1.7976931348623157e+308
nexp =       11   min =        -max
---------------------------------------------------------------

Machine parameters for float128
---------------------------------------------------------------
precision =  18   resolution = 1e-18
machep =    -63   eps =        1.084202172485504434e-19
negep =     -64   epsneg =     5.42101086242752217e-20
minexp = -16382   tiny =       3.3621031431120935063e-4932
maxexp =  16384   max =        1.189731495357231765e+4932
nexp =       15   min =        -max
---------------------------------------------------------------

Machine parameters for float128
---------------------------------------------------------------
precision =  18   resolution = 1e-18
machep =    -63   eps =        1.084202172485504434e-19
negep =     -64   epsneg =     5.42101086242752217e-20
minexp = -16382   tiny =       3.3621031431120935063e-4932
maxexp =  16384   max =        1.189731495357231765e+4932
nexp =       15   min =        -max
---------------------------------------------------------------

Unfortunately, the actual precision of the `long double` type is [platform and compiler dependent](https://en.wikipedia.org/wiki/Long_double#Implementations). — Warren Weckesser, Aug 18 '22 at 12:58
NumPy's `float128` is not platform-independent; see, for example, https://stackoverflow.com/a/17023995/1217358 — Warren Weckesser, Aug 18 '22 at 13:40
Thanks for pointing out this. I have updated my answer to be more precise. — norok2, Aug 18 '22 at 14:09

are numpys dtypes kept even for non-numpy multiplications/divisions/sums?

2 Answers2

Note on Floating Point information