efficient algorithm to compute continued fraction of a sum of continued fractions

[sstarosta]

Implement a method which, taking continued fraction x as an input, computes the continued fraction of (a*x+b)/(c*x+d) using Gosper's algorithm. Call this method with proper arguments upon

sage: 3*x
sage: x+1

(overlaps with the following TODO item of sage/rings/continued_fraction.py: Add Gosper’s algorithm to compute the continued fraction of (ax + b)/(cx + d) knowing the one of x (see Gosper (1972, ​http://www.inwap.com/pdp10/hbaker/hakmem/cf.html), Knuth (1998, TAOCP vol 2, Exercise 4.5.3.15), Fowler (1999). See also Liardet, P. and Stambul, P. “Algebraic Computation with Continued Fractions.” J. Number Th. 73, 92-121, 1998.)

[git]

Branch pushed to git repo; I updated commit sha1. Last 10 new commits:

​171498b	Solved the big bug with the negative value of infinite continued fraction! :)
​7836c3e	Don't know what changed. Backing up.
​7c6b445	Added a few lines to check for problem with singularity.
​56c76e2	Latest version... Don't know what changed :)
​3ab71ab	Done some code refactoring, tried to make the code simpler and more readable. I hope I didn't break it.
​2bd5194	Of course I broke something. Bug for infinite fractions fixed and code shortened.
​0dece4c	Added a test and fixed associated bugs (hopefully).
​e246b63	Added a test and fixed associated bugs (hopefully).
​f9714f8	Added check for arguments' validity.
​3fcc597	Moved Gosper iterator.

[git]

Branch pushed to git repo; I updated commit sha1. New commits:

​51f21d5

Better check coefficients using Integer.

[vdelecroix]

Hello,

Some remarks:

1. You would better edit the description of this ticket to explain what you want to do.

2. Your branch does not merge cleanly with the last beta. The reason is that you added a line in lazy_list.pyx: from sage.rings.all import Infinity. It is always better to base your work on the latest available version of the source code. You can consult the ​developer manual.

3. The commit message should reflect what is inside the commit. Not your mood at the time you did the commit.

4. The file with gosper algorithm should not be in sage/misc/. Put it in the same directory as continued_fraction.py. And you would better call it continued_fraction_gosper.py.

5. In continued_fraction.py you added the following in the comparison code

   @@ -530,6 +530,8 @@ class ContinuedFraction_base(SageObject):
                if a == ZZ_0 and b == ZZ_0 and i:  # rational case
                    return 0
                i += 1
   +            if i == 300:
   +                return 0
   

   I understand why you did it but it is definitely not a solution. What you can do is to raise a RuntimeError after a certain limit.

6. All methods must be documented as explained in details in the ​developer manual.

Vincent

[vdelecroix]

And last but not the least: it would be very nice to have Gosper algorithm within Sage!!

[mirgee]

Hi,

thanks for the great input.

As for 5.: We need a way to compare instances of _infinite, and even though this is not a valid test, if enough terms are the same, the two continued fractions could be regarded as equal for most practical purposes. Another option would be to compare values, as in _real, but the problem is that

1. value() method is by documentation optional,

2. in case of _infinite, this would delegate the computation of value to RLF, which itself may use comparison,

3. this solution would effectively do the same thing - compute the value of the operands with certain precision and then compare the intervals.

What do you think?

[git]

Branch pushed to git repo; I updated commit sha1. New commits:

​fa1173a

Moved and renamed the continued_fraction_gosper. Added some documentation.

[vdelecroix]

Replying to mirgee:

Hi,

thanks for the great input.

As for 5.: We need a way to compare instances of _infinite, and even though this is not a valid test, if enough terms are the same, the two continued fractions could be regarded as equal for most practical purposes. Another option would be to compare values, as in _real, but the problem is that

1. value() method is by documentation optional,

2. in case of _infinite, this would delegate the computation of value to RLF, which itself may use comparison,

3. this solution would effectively do the same thing - compute the value of the operands with certain precision and then compare the intervals.

What do you think?

Everything would be better than an silent approximate comparison. If equality can not be decided just raise an error (see ​Python exceptions). I think that in that case an ArithmeticError would make sense.

If you want to check equality up to some precision, you should do

sage: abs(cf1 - cf2) < 2**-100

or

sage: R = RealIntervalField(256)
sage: R(cf1).overlaps(R(cf2))

[vdelecroix]

For the sake of this ticket I would not implement something like 1+cf or 3*cf which needs some work. Having a method apply_homography would be enough for a first ticket. In order to include this first part in Sage you need to

• make it compatible with the last version of Sage source code (see 2. in 4). By the way why did you make this change in lazy_list.pyx?

• add doctests to each method or function

[mirgee]

The change in lazy_list() is unnecessary and can be deleted. Personally I don't think it is necessary to add doctests for the simple helper functions, but I will comply with the standards and add them, no problem. I will also look into the the comparison.

[git]

Branch pushed to git repo; I updated commit sha1. New commits:

​fa492af	Added more doctests in continued_fraction_gosper.py.
​fea8eb9	Removed unnecessary comments.

[git]

Branch pushed to git repo; I updated commit sha1. New commits:

​92b8aa0

Added more unit tests, replaced overwritten method quotients() with get preperiod length, so that all doctests pass. Added a few more doctests. Now all methods are accompanied with doctests.

[git]

Branch pushed to git repo; I updated commit sha1. New commits:

​e063f2f

Using a set of tuples instead of list of dictionaries for state tracking, should improve performance. Using lazy_list instead of Word for instances of infinite continued fractions.

[chapoton]

rebase, squashed, py3-compatible

---

New commits:

​3308c92

Implement Gosper algorithm for homographic action on continued fractions

[git]

Branch pushed to git repo; I updated commit sha1. New commits:

​031088f

fix for py2: also define next

[git]

Branch pushed to git repo; I updated commit sha1. New commits:

​aeb6da1

Merge branch 'public/ticket/19120' in 9.2.b12

[git]

Branch pushed to git repo; I updated commit sha1. New commits:

​fe2cd08

fix one returns

[vdelecroix]

not suitable yet. I will provide a commit.

[vdelecroix]

Author? (name in the git log is John Doe <miroslav.kovar@concur.com>)

[mirgee]

Replying to vdelecroix:

Author? (name in the git log is John Doe <miroslav.kovar@concur.com>)

I am the author: Miroslav Kovar (miroslavkovar@…)

[git]

Branch pushed to git repo; I updated commit sha1. New commits:

​81e0693

fixes, documentation, copyright, optimization

[vdelecroix]

Replying to mirgee:

Replying to vdelecroix:

Author? (name in the git log is John Doe <miroslav.kovar@concur.com>)

I am the author: Miroslav Kovar (miroslavkovar@…)

Thanks. I put your name and mail in the copyright header of the file.

[chapoton]

patchbot plugins complain

[chapoton]

Le reference [LS1998] a une allure bizarre, avec un P isolé

[git]

Branch pushed to git repo; I updated commit sha1. New commits:

​3323f09

enoding, import, references

[chapoton]

Merci. Pyflakes est encore pas content:

src/sage/rings/continued_fraction_gosper.py:37:1 'sage.rings.real_mpfr.RR' imported but unused

[git]

Branch pushed to git repo; I updated commit sha1. New commits:

​7161eaf

remove unused import

[vdelecroix]

There is something wrong in the way infinity is handled. In the situation where we apply x -> 1 / (qx - p) to x=p/q Gosper iterator is empty. But currently continued_fraction([]) is wrongly 0.

[git]

Branch pushed to git repo; I updated commit sha1. New commits:

​8cafbb4

fix infinities in continued fraction and do not forward values

[vdelecroix]

fixed

[chapoton]

ok, this seems to be ready for a positive review ?