Bug in ContinuedFraction rounding

[kliem]

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]

Changed keywords from continued fraction, rounding to continued_fraction, rounding

[slel]

comment:6

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

[mkoeppe]

comment:7

Moving to 9.4, as 9.3 has been released.

[kliem]

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.

[kliem]

comment:11

For a=39117938827825157072802367/91026195129981723206 there is a ZeroDivisionError.