Bug in ContinuedFraction rounding

kliem commented 4 years ago

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.

CC: @slel @videlec

Component: number theory

Keywords: continued_fraction, rounding

Issue created by migration from https://trac.sagemath.org/ticket/29957

fchapoton commented 4 years ago

Changed keywords from continued fraction, rounding to continued_fraction, rounding

slel commented 3 years ago

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 3 years ago

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 3 years ago

comment:6

Is the element constructor of RealField(prec=17, rnd='RNDZ') buggy?

mkoeppe commented 3 years ago

comment:7

Moving to 9.4, as 9.3 has been released.

kliem commented 3 years ago

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 3 years ago

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
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, (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
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 3 years ago

comment:11

For a=39117938827825157072802367/91026195129981723206 there is a ZeroDivisionError.

mwageringel commented 3 years ago

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.

sagemath / sage

Bug in ContinuedFraction rounding #29957