double float parsing broken

[h-milz]

CircuitPython version

Adafruit CircuitPython 9.0.0-beta.2-25-g95f24318e4-dirty on 2024-03-11; Adafruit Feather M4 Express with samd51j19

I'm still trying to get double float to work. Linking against libm works fine. When MICROPY_FLOAT_IMPL_DOUBLE is 
off, I can do: 

>>> import math as m 
>>> m.sin(4)
-0.756802
>>> m.sin(4.0)
-0.756802
>>> m.sqrt(42)
6.48074

When I activate _DOUBLE and _CMATH, it looks like this: 

>>> import math as m
>>> import cmath as c
>>> m.sin(4)
1.839006457345601e-314
>>> m.sin(4.0)
0.0
>>> m.sqrt(42)
6.146221405944527e-315
>>> c.sin(4)
(-0.7568024953079282+-0j)
>>> c.sin(4.0)
0j
>>> c.sqrt(42)
(6.48074069840786+0j)
>>> c.sqrt(-42)
(3.968309176204584e-16+6.48074069840786j)

where the last real part appears to be a rounding issue. 

Apparently, double float parsing is broken. Integer parsing seems to work ok for cmath entry. Watch the 
difference between c.sin(4), where the real part of the result is correct, and where real mode sin from libm.a 
is invoked under the hood, and m.sin(4) which takes input from the command line, calls the same real 
mode sin, and obviously returns gibberish. 

When the argument is a float (4.0), it is apparently parsed wrongly - c.sin internally calls m.sin as you can tell. 

Similarly for sqrt - m.sqrt returns gibberish, but c.sqrt works fine for both signs. 

I dug through parsenum.c, objfloat.c, formatfloat.c and a couple of more but I can't find anything suspicious. 

(As a sidenote, c.sqrt when called with a complex argument returns nonsense as well: 

>>> c.sqrt(1.234+0.456j)
(0.4774934554525329+0.4774934554525329j)

The correct result should be something like 1.129+0.202j. No idea if this is related to num parsing, but when 
I do the same with _DOUBLE disabled, it works fine.)


### Code/REPL

```python
none

Behavior

see above

Description

see above

Additional information

No response

[h-milz]

Having dug a bit more, I can see a pattern. All functions returning their results via mp_obj_new_float() seem to exhibit this behaviour, which is why complex functions returning complex results appear to work fine whereas cmath functions returning reals (like phase() or polar()) return nonsense.

It appears that float64 results get masked out by 0xffffffff at some point:

c.sin(27)
(0.956375928404503+-0j)
m.sin(27)
1.46631046621851e-314

and decoded according to IEEE754 this contains:

Float64 number: 1.466310466e-314
Sign (hex): 0
Exponent (hex): 000
Mantissa (hex): 00000b0e5ba30

Float64 number: 0.956375928404503
Sign (hex): 0
Exponent (hex): 3fe
Mantissa (hex): e9aa1b0e5ba30

See what happens?

This also happens internally:

c.rect (1.4142135623730951, 0.7853981633974483)
(8.496024063526605e-315+0j)

The result should be (1+1j), but c.rect uses m.sin and m.cos internally, which in turn return their results using mp_obj_new_float(). (8.496024063526605e-315 is a masked out sqrt(2)).

mp_obj_new_float() is declared in py/obj.h for float and double separately (although the various ifdef levels leave me wondering ... ), and there is a function declaration in py/objfloat.c.

So presumably num parsing is not the (only) problem, but output formatting of doubles (although float_make_new() also uses mp_obj_new_float() in some cases).

At this point I'm lost for now. I'd like to add some debugging printfs but failed to find out how this works. Hints appreciated.

[tannewt]

Sounds like you may be hitting the object representation: https://github.com/adafruit/circuitpython/blob/main/py/mpconfig.h#L121

CircuitPython uses C which packs floats into the 32 bit pointers instead of allocating them separately to the heap.

[h-milz]

@tannewt Thank you for pointing this out! I see where this is going to.

When I set MICROPY_OBJ_REPR to REPR_D in mpconfigboard.h, the (mp_obj_t) casts in circuitpy_mpconfig.h:230 lead to an illegal pointer to int error. I replaced them with the (cleaner) MP_OBJ_FROM_PTR macro, which compiles cleanly with REPR_C, and removes the illegal pointer to int error with REPR_D, but now throws:

../../py/obj.h:292:28: error: initializer element is not constant
292 | #define MP_OBJ_FROM_PTR(p) ((mp_obj_t)((uintptr_t)(p)))

which is correct because at compile time, the value of p is not known. With REPR_C, the 32-bit pointer p can be allocated right away, but with REPR_D, the then 64-bit pointer needs to be allocated on the heap, and hence the address is unknown until the final linker stage. Bummer.

Too bad, because float64 and REPR_D references are all over the place, suggesting that one could use them, but apparently this is not (no longer?) the case. Any idea how to get this done?

[tannewt]

Too bad, because float64 and REPR_D references are all over the place, suggesting that one could use them, but apparently this is not (no longer?) the case. Any idea how to get this done?

Why do you want double support? What port are you targeting?

MicroPython has lots of settings but we only use a subset of them in CircuitPython. I'm not sure if they have any builds that do double. I only see references to it in the unix port. I know REPR_D is generally used for 64-bit builds though.

[dhalbert]

Use ...OBJ_REPR_A and MICROPY_FLOAT_IMPL_DOUBLE. That shold allocate 64-bit floats on the heap. As @tannewt mentioned, OBJ_REPR_D is used on 64-bit architectures, and we don't know that it would function properly on at 32-bit machine.

Note that SAMD51 does not have hw support for 64-bit floats, so it will be noticeably slower.

[h-milz]

@tannewt The use case is a small scientific calculator based on a Keyboard Featherwing and a M4 Express or RP2040. While 6-7 digits display resolution would be acceptable, follow-up errors due to rounding would not. So I experimented with Jeff's udecimal and utrig libraries and wrote hyperbolic functions and complex decimal math extensions, which works fine so far. But memory is extremely tight even on the RP2040, and I would still like to squeeze e.g. Fractions in. And we haven't talked about numerical integration yet. Maybe an ESP32-S2 board would work.

@dhalbert Thank you, this works fine. Comparing the internal libm with the GNU libm, I noticed that the internal libm sometimes has rounding errors in the last 2 or 3 digits, so I'll go for the GNU libm I suppose. Concerning efficiency, yes, having the floats on the heap will definitely be slower than 32 bit floats in the pointers, but still this custom solution is way more efficient than calculating, say, the cosh of a complex argument in Python decimal. :-)