sagemath / sage

Main repository of SageMath
https://www.sagemath.org
Other
1.42k stars 476 forks source link

Number field to float conversion gets sign wrong #37983

Open saraedum opened 5 months ago

saraedum commented 5 months ago

Steps To Reproduce

In a lengthy computation I manage to create a broken number field (element):

sage: x = loads(b'x\x9c\x9dV\x89\x7f\x13E\x14n+ ,\x08(\x82\x17j\xbdS\x95\xda\xa6\xb7\'\x95\x02\x96J\xc4\x07h<p\xdcl&\x99I\xf6\xc8\xb7;\x1b\xa8Z\xc1c\xb3\xd4\x1b\xbc\xef\xfb>\xfeFg7[HCA\xe4\x97\xdf/\x997\xf3\xde7\xdf|\xf3\xde\x9b\x9c\xe8\xb3\x02\xb3\xca\x07}\xe9V\x83A7tJ\xdcg\x15\xc9\xed\xf22\x83q\x9b;\xdcU\x0c\xa1Y\xf6M%-\x831\xc7\xacsVH\xbd\xf6$N\xbb\xbb}\x86\r\xf4\xe4\xda\xf8\x81\xf2CK\x85>\x1f\xac\x98\x96\xf2\xfcy\xa3\xca]\xeeK\x8be6\x0b\xdd\x86\xb4\xea67\xd0{\xde\x18\xdb\xf3\xeaa\x83Um\xafd\xda\x06\xfa\x8aS===\x17\xc3\x7f\xb0\x83&kr?\x90\x9e\x9b\xc7e\x11V\x11V\xcf\x19s}-\xac\xc9u*\x91\xf0\xf7\\\xd3n\x87\x1b\x94\x99)\x80\x81\xcb\x07\x08k;\xdd\x1b\x9e=\xefz\x8e4\xed\x8e!;\x0bbKW\x19\x07\xce\xbb\x82u\xb9\x8b@\xd3k\xcc\xf2\xdc\xb6,\x9eo,)\xc6\xce\x02\x93\xf61`\xe4\xc4\xdab\xaf\x96f\x1e\xeb#l(,*\\A\xd8x\x04\x9bV<\xa3\x91^\xe4\x19\x0b\x9b\x8b}:vG\x1eWF\xb8\x8a\xb0ElN\xc1\x86pu\x84\xad\x84m\x99=\x8ck"\\K\xb8\x8e/\xea\r\xae\'\xdc\x90N\x9b\xd8\x1e\xe1\xc6\xee\x1b,\xcf\xbb\xa6\xa3o\xdb\xb2\xcd 0\x96[\xb8)W\xec\xd7\x91\xd3v\x95\x97|SZ\xc43\xa5\xd9Q\xa9\x04\xb3L\xc5\xab\xfa\xfaqs\'{\xa0d\xfa\xc6\xb91\x06\xfa-\xc6J\xa1\xb4\x95t\x19\xd3y\xa6L\xa5|\x03\xb7ts\n]\x89P\x1f\x9c7|\x1e\xe8\xc4M\x05\xd0\xaa\xfa\xbc\x1cZ:\x0fom\x07d\xdbK\x1e,\r\xe7\x8d}\x9etw-\x19\xb8\xad\xed\xe8\xc8\xc0\xd2.\xb6m$_\xcc\xe1Jx\x9a\xce\xed\xa2\xfbs.\xb0cV\x1d3`\xa6[ff\xb9,\x95lr\xd6\x9e3\xf6\xa7?\xd3ny:[hO\x18\xb8c`\x01w\xc6\xb8\x8b\x90+\xae\xd7\xf2\xb5\xc52\x8fI\xcf\xc1@q\x83\x9e\x99\x91\xfa\xb0\xb2\x14&Q\xb8;\xc6=\x84{\xc5@qK\xa2u\x066\x1d\x04\x9e%\xcd\xd4cG\x8cA\xc2}]\x1e\xbb<\xc7\tU\xdbc(\xc60!\xaf=6vx\x1cv\xa52m\x8c\xc4\x18%\x8c\xe9\xc5\x84M\'\xf0x\x8c\t\xc2\xa4^\xd9\xd4\x116\xeb&\x85\xc81\x15\xe3~\xc2\x03zu\x8d^\xcd\xc0\x1e\x8c\xf1\x10\xe1a=\xb96=HS&%\x8bGb\xec$Lg{tR{4\xc6.\xc2\xcc\x8a\xea\x06\\\x05I\x16\x19\x07\xf5\xc0\xc0n-\xdc\x9e\x18{\t\x8fe\xf8\xb3nE\xea}9fc\xec#\xcc\xb5\xf0x\x84\xfd\x0b(\xc4x\x82p \x15\xb3a\xfaIsKS\x16O\xb6@\x84\x83-\x1c*\x14D\xbf\xc2a\xc2S\xba\'<\xfd\x1f\tj\xa0\xb8\x82Ga\xf71\xc5\xdd\xe4|Y\x075\xf0\xcc\x05\x80\xf6&\x8d\xd3L\x1a\x00\x9e\x15\xbd9\xb1Z\xac\xd1\xf5.\xd6\xe5\x84\xae\xf0\xe7r\xe2L\xf9>\x1f\xe1\x08\xe1\x05\xb1Y\xe8\xc2e\x84\x17\xf5H\x97\xacI(\xa5%k\x11\xcaB\x17+/\x88\xed\x85\xc5\x82Be\x01\xd5\xe2j\x1dl\tn\xd5!\x16\x03\x05I\xa8\xad\xc4\x99<O\xf3\xac_\x80\xe7\xd9\xaet\xc87\xad:\xd7|\xed\x8c\xa5\xd3\xc1\xd2\x8d\xe0\x11\x1am\x96 \xf8m\x96\x01A\xa5,CB3\xc2Q\xc2\xb1e\xbdK\xab\xc9\x9cFE\xeaw\xc8\xd2\x86\xe2\x8c%\n\xcf\xba\x8a\xfb\xcd\xac\x15,=HY\xc3\xd7\xcf\xd1\xfc\xff\xc6X\n\x1e2\xf0\xd2\xfe\x9e\xde\xc5\x16^&\xbc\xb2\x12\x8c\xdf\x05\xd3~r:\xe2\x17.&\xaa{\xd3W\x93M\x8b\xab\xb4VT\x989\x8c\xe31N\x10^+N$\rx\xb0<9\x1aZ\x8df\xcd\xad\x8dV\x85\xe5\xfa\xf9r\x8d\xdbc%5\xe1\xa9\xb2\xaa\x96\xec\xfc\x94o\x0f\x8fs\x8c)\xab\xea\x0c\xd7&\xc6\'\x9d\x9dCx}\xae?\xc6\x1b\x847\xc5\x82X\x06?\x83(F\x8b\x10_2\xbc\xad\xe1O&\xf0\x8b\x84\xb7b\xbcMxg\xae7\xc6\xbb\x84\xf7Zx\x9f\xf0\x81\xe8\xc9\x89\xda\xa5\xfe\xe5\xb8\xd0\x7f\r\x16\xc0\xd7\ty\xaa\x13[?\xad\xbc\xaa\xf3.}\xdf\x96\x0c\x9cN\x12\xedC\xc2G\xe2t\x92h\x1f\x13>i\x8f>%|\xa6\xf09\xe1\x8b\x16\xbe$|\x15\xe1k\xc27-|\x9b\x95\xc7w\x0b\xf8^\x88\x93\xba(~ \xfc(N\x89\xd3\xc5\xadZ\xa7\x91\x86\xeb\xf0\xd0\x1aj\xd6\xa6\x94\xa3x\xad\xec7\xc6\x04~\x8a\xf03\xe1\x17\xed\xb3-\xcds\'o\xdb\xf5\x91\xdax\x98\x9f\xb2\xea\xc1\xe4\xc8pm$\xc0\xaf\x11~#\xfc\x9e\x01U\xc6\xc3\x86\xdf\xa8\x8e\x84b\xb4\xe2\xe7=9Z\x9d\xc0$\xfe\x88\xf0\'\xe1/\x85\xbf\t\xff\x0c\xfe\x0b\xd5\x9d{1')
sage: x.parent()
Number Field in a with defining polynomial y^2 - 2 with a = 1.414213562373095?
sage: x
-549013069592069/3105686915748774*a + 4822645609326531104236812231857/19290582437306264907490165006152
sage: x > 0
True
sage: float(x) < 0
True
sage: RealField(63)(x) < 0
True
sage: RealField(64)(x) > 0
True

Note that somehow the number field is broken. Copying the element over into a "fresh copy" of the quadratic number field, fixes things:

sage: R.<a> = QuadraticField(2)
sage: y = R(str(x))
sage: x == y
sage: float(y) > 0
True

Expected Behavior

The float of a positive number should be non-negative.

Actual Behavior

The float of a positive number is negative.

Additional Information

I couldn't find an issue like that reported. I am not sure how complex embeddings work which seem to be the underlying machinery here. Maybe it's another instance of interval arithmetic not doing the right thing?

Environment

I tried this with SageMath 10.3 from ArchLinux and SageMath 10.2 from conda-forge.

Checklist

saraedum commented 5 months ago

@videlec this is biting me in the harmonic differentials code of sage-flatsurf.

saraedum commented 5 months ago

Here's a simplified reproducer without number fields:

sage: x = loads(b'x\x9c\x9dTi\x97\x13E\x14\xcd\x0c\xa0\xd0\xca"2\xb8+\xa0@@\x89\xd9&\x0bn\x803 \x06\xc2X\x02\xb6\x0b\x96\x9dJ\xa5\xab\x93\xea\xe5V\xaa\xe78\x9e\xe39~\xe9\xe4\xe4_[\xddI\x98\xcc\x18\xe6\x8c|\xabWu\xdf}\xef\xdd\xbaU\xff\xac\xb2\xa1\xe3\xf2\x82\xf2\x02wXP\xdc\x91T:\x7f\xedX\xbe3\xe0\x94K\xee\xf3@[\xc8-\x07\x11\xb3z`\x16w=.\xbb\x16V\xae\x11\xec\xa1\x03:\x8e\xb2nK\x97w\x94\xe3\xb1\x14n\xe1\xc8\x12D{\xf3O\xcd\x83\xa1\x17\x06\x9b\xf3\x92G\x0f \xba\xc7\x03\xae\x1c\x1d*\x0b\xc7\xa6\xb0\xa1V1\xd3\xb1\xe2\x85\x9e\xc3\xcc\xc1\x8e\xe5\xa6\x18\x8f\xd1YL\xe3 \xf2\xd8@r\x0b\xaf\xe4_\x94#\xc3p\x10G\xd4\x95a\'\xed\xf4U\xbb\x99\xcb\xe5\x16\xda\x08b\xbf\xc3\x15\xed\xa5\xe3\xee\r\xdaY\x90\xe9@\xb7\xb9J\')\xe3x\x82\x13\x04V\xcbj\xad\x8e\xf0Z~\x8f\x86\x8e6\x10\xa3c\x96n\x91Y8\x13\xf2u#\xe4\xc9Ex\x14\xca\x9d \xf4=G.,\xe9.\x89\xf4\x8cd[/<\xc1\xa9\xfc!\xd8\xcc\x19ea0\x95\xc5h;W\x8c\xee\x12\x13\x83\xb1p:/N\xda+F\x9a\x1d\x9cI\xf0F{\xa2q\x96\xe0\xcdg8\xb7t\xc6\xa9\x97\x9eGX\xb3WM\xee\x8d2\xce\'x\x8b\xe0m\xb1\x96\x91\x15\xf1N\x82w\t\xde\x9b\xc5%\xbc\x9f\xe0\x03\x82\x0f\xf9\xc4\x14\xf8\x88\xe0B\xb6\xed\xe0b\x82Kmq\xb1=ik|\xfc7>\xb1\x8f\x99m&8\x1b\xe0\xf2d\xa8q\x85\xe0\xea2\x8f\x9104R\xe4\x0f\xf0\xd5\xee\xa0\x8f\x95\xc3\x06\xdc\xf8\xeb\x9a8\x95\x17f\xb4\xeby\xf1\xbc\xefO\x13|FpC\xac\t\xd3q\x81\xe0s\xb32\xbd\x16\tJY\xafe\x82J\x82*\xc1\xfa\x7f\x9e\x8d\x1f\xf5<\x8bRf\x02\xcd)M_\xc4\xfd@s\xb5=\xbb\xfb\x99\xfd\xe7\x1e*Y\xa8\xfdo\x8eyr\xd1B\xfdane2B\x83\xa0\xb9\x8cF\xed\xa3\x99\xbax!\xff\xe6a\xb2\xf6\x17\xfd"-j\x1f5Z\x91\xf6\xc6\x13|9\xc6W\x04_\xdb\xf5\xf4N\x0b\xddF5f\xd1v?\xe8W]\xc1\x02U\xee\xf6\xb9\\\xef\xe8z\xa8\xbb\xda\xed\xc8rS\xc9R\x8dc]3\xd7/\xf5\xeb\xb5\x86\x7f\xab\x88oZ\x17\xc6\xb8Ep[\xdc\x14{\xe87pg\x8co\t6^\x9a^\x1a\xfa\xcd\x94\xfe.\xc1\xbd1\xbe#\xb8\xdfZ\x19\xe3{\x82\xd6\x08\x0f\x08\x1e\xb2\xc3|\x00\xf3\xaf\x92"v\xba\xa9\xd5\x99\x11)\xf3\xfd\xc2\xcf\xb0\xb9\x1fc\xae\xb7\x9d\x17W_\xb6\xc2A\xd4t\x08e\xdc\xfeh\x91\xdb|\x05\xdc5\xa6\xce\xfa\x9a\x07\xd8\xb2\xd7\x8ct\x95(\xf0y\xcc\x8a\xdb\xfd\xa6\xf65\xefwU\xb4.\xf0C\x02B\xf0\xa3\xd8\xb2\xcfg\xd6\xf7\xcbR\x0e*\xfdZ\\n\xb2\xc1\xb0Q)\xf5+C<N\xf0\x84\xe0\xa9\x98\x12\xf5jq\xa4"\xb7\x12\x8bjO\x95C\xaf\xea\xd6\xd1\xc0O\tl\x82\x9f5~!\xf8u\x84\xdf\x08\x9e%\xf8\x9d\x80\x8e\xf0\x07\x81S\xf8\x17T\xa3M\xa8')
sage: x.parent()
Real Lazy Field
sage: x > 0
True
sage: float(x) < 0
True
saraedum commented 5 months ago

Note that float(x) is implemented as float(x._value) and ._value is an algebraic real. So maybe this is related to #37927?

saraedum commented 5 months ago

Then again, ._value is just a ANExtensionElement here. So maybe that's unrelated.

videlec commented 5 months ago

Note that float(x) is implemented as float(x._value) and ._value is an algebraic real. So maybe this is related to #37927?

I do not think so (I am following your second definition of x). Even though the enclosure is super close to zero and the minimal polynomial gigantic it looks correct

sage: p = x._value.minpoly()
sage: a = min(p.roots(RealIntervalField(2048), False))
sage: a.intersection(x._value._value)  # one would get an error if there was no intersection
6.62...e-30
videlec commented 5 months ago

Though the enclosure contains zero

sage: x._value._value.contains_zero()
True

which makes it impossible (without further refinement) to determine the sign.

videlec commented 5 months ago

In the description you wrote

sage: RealField(63) < 0
True
sage: RealField(64) > 0
True

Maybe you wanted instead

sage: RealField(63)(x) < 0
True
sage: RealField(64)(x) > 0
True
videlec commented 5 months ago

For the second x you defined, I get a different behavior

sage: float(x) > 0
False
sage: x > 0
True
sage: float(x) > 0
True
videlec commented 5 months ago

And in my case, the real lazy field is not involved

sage: float(x._value) > 0
False
sage: x._value > 0
True
sage: float(x._value) > 0
True

If we want the sign of float(x) to be correct one needs a smarter AlgebraicReal.real_number. Though I am not sure we want to slow down such method.

videlec commented 5 months ago

The value changes quite a bit in between

sage: float(x)
-6.776263578034403e-21
sage: x > 0
True
sage: float(x)
6.629999776684485e-30
videlec commented 5 months ago

It seems to me that the description of the method AlgebraicReal.real_number is wrong here

    def real_number(self, field):
        """
        Given a :class:`RealField`, compute a good approximation to ``self`` in
        that field. The approximation will be off by at most two
        ulp's, except for numbers which are very close to 0, which
        will have an absolute error at most
        ``2**(-(field.prec()-1))``. Also, the rounding mode of the
        field is respected.
        ...
        """
videlec commented 5 months ago

My conclusion (for the time being):

saraedum commented 5 months ago

@videlec I see. I expected it to be the best approximation as a Python float but I see that there seems to be simply no documentation about what __float__ actually does.

What's weird is that there is a coercion from my original x's parent to RDF.

sage: RDF.coerce(x) < 0
True

Now, coercion to inexact rings is a weird animal anyway. But this feels wrong then.

sage: RDF.coerce(x) == RDF.coerce(QuadraticField(2, embedding=1).coerce(x))
False
videlec commented 5 months ago

It is true that the existence of a coercion suggest that RDF(x) should not depend on the approximation carried by x. In particular, the behavior of https://github.com/sagemath/sage/issues/37983#issuecomment-2105982840 looks wrong in that direction.

I think that based on your example, one can build two exact numbers x and y with the same parent, that compare equal and such that RDF(x) and RDF(y) are different (this is pretty straightforward to do in QQbar using the expression trees).

saraedum commented 5 months ago

Just to clarify, I can work around this in the context of harmonic differentials. I just don't trust the sign of floats anymore.

videlec commented 5 months ago

I think that the behavior of float(x) and RealField(prec)(x) ought to be specified in sage. And as you mentioned, the coercion tend to suggest that it should be the closest element in the target.

To do so, one would need to plug in the right places the single function

Maybe such function already exists somewhere in sage?

videlec commented 5 months ago

But then come the question: what should be RBF(x) and RIF(x)? Since there are coercions AA -> RBF and AA -> RIF one could argue the same. Computing the closest interval certainly makes sense. But the closest ball seems like a weird concept to me. And beyond that, it is not completely clear to me whether this is the desirable behavior for RBF and RIF. Many times, you just want a ball or interval that is computed at the precision of RBF or RIF. This is what is done when you call functions as in RBF(1/3).zeta().

saraedum commented 5 months ago

Why should there be a coercion in these cases? (I mean, apart from not breaking existing code.)

videlec commented 5 months ago

Why should there be a coercion in these cases? (I mean, apart from not breaking existing code.)

That is a good question. The main purpose of coercion is to have $x + y$ works when $x$ and $y$ have different parents. But I do not really see a use case here.

One could also question the purpose of coercions AA -> RR and AA -> RDF.