Bug in ContinuedFraction rounding
Opened this issue · 9 comments
kliem commented
Rounding with RNDD and RNDZ with 17 bits of precision doesn't work for the continued fraction of 1761/1024:
sage: fields = []
sage: for rnd in ['RNDD', 'RNDZ']:
....: fields.append(RealField(prec=17, rnd=rnd))
sage: a = 1761/1024
sage: cf = continued_fraction(a)
sage: for R in fields:
....: if R(cf) == R(a):
....: print('this worked')
This contradicts a doctest in src/sage/rings/continued_fraction.py which claims pretty much that this always works (running a loop on 3000 random values with assertion error; but the values were never really random).
Looking at this example in more detail:
sage: Rd, Rz = fields
sage: ad, az, bd, bz = Rd(a), Rz(a), Rd(b), Rz(b)
sage: ad, az, bd, bz
(1.719, 1.719, 1.719, 1.719)
sage: [x.sign_mantissa_exponent() for x in (ad, az, bd, bz)]
[(1, 112704, -16), (1, 112704, -16), (1, 112703, -16), (1, 112703, -16)]
On top of that, this test can result in the following errors:
sage -t --long --warn-long 62.5 --random-seed=3013 src/sage/rings/continued_fraction.py
**********************************************************************
File "src/sage/rings/continued_fraction.py", line 649, in sage.rings.continued_fraction.ContinuedFraction_base._mpfr_
Failed example:
for n in range(3000): # long time
a = QQ.random_element(num_bound=2^(n%100))
if a.denominator() % 8 == 0: # not precices enough # :trac:`29957`
continue
cf = continued_fraction(a)
for R in fields:
assert R(cf) == R(a)
Exception raised:
Traceback (most recent call last):
File "/amd/compute/mwagerin/git/sage_compute/sage2/local/lib/python3.9/site-packages/sage/doctest/forker.py", line 718, in _run
self.compile_and_execute(example, compiler, test.globs)
File "/amd/compute/mwagerin/git/sage_compute/sage2/local/lib/python3.9/site-packages/sage/doctest/forker.py", line 1137, in compile_and_execute
exec(compiled, globs)
File "<doctest sage.rings.continued_fraction.ContinuedFraction_base._mpfr_[25]>", line 7, in <module>
assert R(cf) == R(a)
File "sage/structure/parent.pyx", line 898, in sage.structure.parent.Parent.__call__ (build/cythonized/sage/structure/parent.c:9338)
return mor._call_(x)
File "sage/structure/coerce_maps.pyx", line 287, in sage.structure.coerce_maps.NamedConvertMap._call_ (build/cythonized/sage/structure/coerce_maps.c:6042)
cdef Element e = method(C)
File "/amd/compute/mwagerin/git/sage_compute/sage2/local/lib/python3.9/site-packages/sage/rings/continued_fraction.py", line 704, in _mpfr_
return R(sgn * m_even) >> N
File "sage/rings/real_mpfr.pyx", line 2710, in sage.rings.real_mpfr.RealNumber.__rshift__ (build/cythonized/sage/rings/real_mpfr.c:21059)
return x._rshift_(Integer(y))
File "sage/rings/real_mpfr.pyx", line 2691, in sage.rings.real_mpfr.RealNumber._rshift_ (build/cythonized/sage/rings/real_mpfr.c:20906)
mpfr_div_2exp(x.value, self.value, n, (<RealField_class>self._parent).rnd)
OverflowError: can't convert negative value to unsigned long
**********************************************************************
sage -t --long --warn-long 62.8 --random-seed=3137 src/sage/rings/continued_fraction.py
**********************************************************************
File "src/sage/rings/continued_fraction.py", line 649, in sage.rings.continued_fraction.ContinuedFraction_base._mpfr_
Failed example:
for n in range(3000): # long time
a = QQ.random_element(num_bound=2^(n%100))
if a.denominator() % 8 == 0: # not precices enough # :trac:`29957`
continue
cf = continued_fraction(a)
for R in fields:
assert R(cf) == R(a)
Exception raised:
Traceback (most recent call last):
File "/amd/compute/mwagerin/git/sage_compute/sage2/local/lib/python3.9/site-packages/sage/doctest/forker.py", line 718, in _run
self.compile_and_execute(example, compiler, test.globs)
File "/amd/compute/mwagerin/git/sage_compute/sage2/local/lib/python3.9/site-packages/sage/doctest/forker.py", line 1137, in compile_and_execute
exec(compiled, globs)
File "<doctest sage.rings.continued_fraction.ContinuedFraction_base._mpfr_[25]>", line 7, in <module>
assert R(cf) == R(a)
File "sage/structure/parent.pyx", line 898, in sage.structure.parent.Parent.__call__ (build/cythonized/sage/structure/parent.c:9338)
return mor._call_(x)
File "sage/structure/coerce_maps.pyx", line 287, in sage.structure.coerce_maps.NamedConvertMap._call_ (build/cythonized/sage/structure/coerce_maps.c:6042)
cdef Element e = method(C)
File "/amd/compute/mwagerin/git/sage_compute/sage2/local/lib/python3.9/site-packages/sage/rings/continued_fraction.py", line 709, in _mpfr_
assert m_even / (ZZ_1 << N) <= p_even/q_even
File "sage/rings/integer.pyx", line 2040, in sage.rings.integer.Integer.__truediv__ (build/cythonized/sage/rings/integer.c:14337)
raise ZeroDivisionError("rational division by zero")
ZeroDivisionError: rational division by zero
The overflow can be reproduced like this:
fields = [RealField(prec=17, rnd=rnd) for rnd in ['RNDN', 'RNDD', 'RNDU', 'RNDZ', 'RNDA']]
a = -47866071760720267/173220919737
cf = continued_fraction(a)
[R(cf) == R(a) for R in fields]
The zero division error like this:
fields = [RealField(prec=17, rnd=rnd) for rnd in ['RNDN', 'RNDD', 'RNDU', 'RNDZ', 'RNDA']]
a = -4330659901673730869863039591/17079070615116212716183
cf = continued_fraction(a)
[R(cf) == R(a) for R in fields]
In #29979, a doctest was marked not tested because of this.
Component: number theory
Keywords: continued_fraction, rounding
Issue created by migration from https://trac.sagemath.org/ticket/29957
fchapoton commented
Changed keywords from continued fraction, rounding to continued_fraction, rounding
slel commented
Description changed:
---
+++
@@ -11,4 +11,4 @@
....: print('this worked')
```
-This contradicts a doctest in `src/sage/rings/continued_fraction.py` which claims pretty much that this always works (running a loop on 3000 random values with assertion error; but the values where never really random).
+This contradicts a doctest in `src/sage/rings/continued_fraction.py` which claims pretty much that this always works (running a loop on 3000 random values with assertion error; but the values were never really random).slel commented
Description changed:
---
+++
@@ -4,7 +4,7 @@
sage: fields = []
sage: for rnd in ['RNDD', 'RNDZ']:
....: fields.append(RealField(prec=17, rnd=rnd))
-sage: a = 1761/1024
+sage: a = 1761/1024
sage: cf = continued_fraction(a)
sage: for R in fields:
....: if R(cf) == R(a):
@@ -12,3 +12,14 @@
```
This contradicts a doctest in `src/sage/rings/continued_fraction.py` which claims pretty much that this always works (running a loop on 3000 random values with assertion error; but the values were never really random).
+
+Looking at this example in more detail:
+
+```
+sage: Rd, Rz = fields
+sage: ad, az, bd, bz = Rd(a), Rz(a), Rd(b), Rz(b)
+sage: ad, az, bd, bz
+(1.719, 1.719, 1.719, 1.719)
+sage: [x.sign_mantissa_exponent() for x in (ad, az, bd, bz)]
+[(1, 112704, -16), (1, 112704, -16), (1, 112703, -16), (1, 112703, -16)]
+```slel commented
comment:6
Is the element constructor of RealField(prec=17, rnd='RNDZ') buggy?
mkoeppe commented
comment:7
Moving to 9.4, as 9.3 has been released.
kliem commented
comment:8
Turns out there are failures for 8 a divisor of the denominator.
Some examples:
32771/8, 16391/16, 8195/32, 4115/64, 2051/128, 1091/256, 515/512, 259/1024, 141/2048, 79/4096, 35/8192, 25/16384, 23/32768, 7/65536, 3/131072, 5/262144, 3/524288.
mwageringel commented
Description changed:
---
+++
@@ -23,3 +23,90 @@
sage: [x.sign_mantissa_exponent() for x in (ad, az, bd, bz)]
[(1, 112704, -16), (1, 112704, -16), (1, 112703, -16), (1, 112703, -16)]
```
+
+---
+On top of that, this test can result in the following errors:
+
+```
+sage -t --long --warn-long 62.5 --random-seed=3013 src/sage/rings/continued_fraction.py
+**********************************************************************
+File "src/sage/rings/continued_fraction.py", line 649, in sage.rings.continued_fraction.ContinuedFraction_base._mpfr_
+Failed example:
+ for n in range(3000): # long time
+ a = QQ.random_element(num_bound=2^(n%100))
+ if a.denominator() % 8 == 0: # not precices enough # :trac:`29957`
+ continue
+ cf = continued_fraction(a)
+ for R in fields:
+ assert R(cf) == R(a)
+Exception raised:
+ Traceback (most recent call last):
+ File "/amd/compute/mwagerin/git/sage_compute/sage2/local/lib/python3.9/site-packages/sage/doctest/forker.py", line 718, in _run
+ self.compile_and_execute(example, compiler, test.globs)
+ File "/amd/compute/mwagerin/git/sage_compute/sage2/local/lib/python3.9/site-packages/sage/doctest/forker.py", line 1137, in compile_and_execute
+ exec(compiled, globs)
+ File "<doctest sage.rings.continued_fraction.ContinuedFraction_base._mpfr_[25]>", line 7, in <module>
+ assert R(cf) == R(a)
+ File "sage/structure/parent.pyx", line 898, in sage.structure.parent.Parent.__call__ (build/cythonized/sage/structure/parent.c:9338)
+ return mor._call_(x)
+ File "sage/structure/coerce_maps.pyx", line 287, in sage.structure.coerce_maps.NamedConvertMap._call_ (build/cythonized/sage/structure/coerce_maps.c:6042)
+ cdef Element e = method(C)
+ File "/amd/compute/mwagerin/git/sage_compute/sage2/local/lib/python3.9/site-packages/sage/rings/continued_fraction.py", line 704, in _mpfr_
+ return R(sgn * m_even) >> N
+ File "sage/rings/real_mpfr.pyx", line 2710, in sage.rings.real_mpfr.RealNumber.__rshift__ (build/cythonized/sage/rings/real_mpfr.c:21059)
+ return x._rshift_(Integer(y))
+ File "sage/rings/real_mpfr.pyx", line 2691, in sage.rings.real_mpfr.RealNumber._rshift_ (build/cythonized/sage/rings/real_mpfr.c:20906)
+ mpfr_div_2exp(x.value, self.value, n, (<RealField_class>self._parent).rnd)
+ OverflowError: can't convert negative value to unsigned long
+**********************************************************************
+```
+
+```
+sage -t --long --warn-long 62.8 --random-seed=3137 src/sage/rings/continued_fraction.py
+**********************************************************************
+File "src/sage/rings/continued_fraction.py", line 649, in sage.rings.continued_fraction.ContinuedFraction_base._mpfr_
+Failed example:
+ for n in range(3000): # long time
+ a = QQ.random_element(num_bound=2^(n%100))
+ if a.denominator() % 8 == 0: # not precices enough # :trac:`29957`
+ continue
+ cf = continued_fraction(a)
+ for R in fields:
+ assert R(cf) == R(a)
+Exception raised:
+ Traceback (most recent call last):
+ File "/amd/compute/mwagerin/git/sage_compute/sage2/local/lib/python3.9/site-packages/sage/doctest/forker.py", line 718, in _run
+ self.compile_and_execute(example, compiler, test.globs)
+ File "/amd/compute/mwagerin/git/sage_compute/sage2/local/lib/python3.9/site-packages/sage/doctest/forker.py", line 1137, in compile_and_execute
+ exec(compiled, globs)
+ File "<doctest sage.rings.continued_fraction.ContinuedFraction_base._mpfr_[25]>", line 7, in <module>
+ assert R(cf) == R(a)
+ File "sage/structure/parent.pyx", line 898, in sage.structure.parent.Parent.__call__ (build/cythonized/sage/structure/parent.c:9338)
+ return mor._call_(x)
+ File "sage/structure/coerce_maps.pyx", line 287, in sage.structure.coerce_maps.NamedConvertMap._call_ (build/cythonized/sage/structure/coerce_maps.c:6042)
+ cdef Element e = method(C)
+ File "/amd/compute/mwagerin/git/sage_compute/sage2/local/lib/python3.9/site-packages/sage/rings/continued_fraction.py", line 709, in _mpfr_
+ assert m_even / (ZZ_1 << N) <= p_even/q_even
+ File "sage/rings/integer.pyx", line 2040, in sage.rings.integer.Integer.__truediv__ (build/cythonized/sage/rings/integer.c:14337)
+ raise ZeroDivisionError("rational division by zero")
+ ZeroDivisionError: rational division by zero
+```
+
+The overflow can be reproduced like this:
+
+```
+fields = [RealField(prec=17, rnd=rnd) for rnd in ['RNDN', 'RNDD', 'RNDU', 'RNDZ', 'RNDA']]
+a = -47866071760720267/173220919737
+cf = continued_fraction(a)
+[R(cf) == R(a) for R in fields]
+```
+
+The zero division error like this:
+
+```
+fields = [RealField(prec=17, rnd=rnd) for rnd in ['RNDN', 'RNDD', 'RNDU', 'RNDZ', 'RNDA']]
+a = -4330659901673730869863039591/17079070615116212716183
+cf = continued_fraction(a)
+[R(cf) == R(a) for R in fields]
+```
+kliem commented
comment:11
For a=39117938827825157072802367/91026195129981723206 there is a ZeroDivisionError.
mwageringel commented
Description changed:
---
+++
@@ -110,3 +110,4 @@
[R(cf) == R(a) for R in fields]
```
+In #29979, a doctest was marked `not tested` because of this.