hsinhaoyu / cont_frac

Bill Gosper's continued fraction arithmetic implemented in Python

In 1972, mathematician Bill Gosper described an ingenious method for performing arithmetic on continued fractions (see here for his original publication). This repo provides a Python implementation, and an article to explain the logic.

About the source code

The code was written as an exercise in literate programming. The Python code was written as parts of an article, using the orgmode markup, which is popular among Emacs users. The main article is in the org/ directory. The Python code blocks can be extracted from the .org files, and assembled into normal Python code in src/. In Emacs, this process is done with the org-html-export-to-html command. On a computer with Emacs installed, you can also run make tangle.

cont_frac's Issues

add plotting

run all demo blocks again

Apply the same tests to two functions

Appendix is not well-documented

Add CI/CD to deploy document to Github Pages

add a new intro section

use local css instead of the ones on GitHub

remove Mathematica code

Numpy type checking

Write a new blog post

No syntax highlight in the exported HTML

Figure out how to set up Github pages

make html doesn't work

all that's needed is (setq org-src-fontify-natively t) in init file?

typing is messy and spotty in appendix

Undetected Overflow

There are situations where a local integer value >= 2^63 is not correctly handled. This is demonstrated by the following sequence, which attempts to calculate the square-root of 5:

$ git clone https://github.com/hsinhaoyu/cont_frac.git hsinhaoyu_cont_frac
Cloning into 'hsinhaoyu_cont_frac'...
remote: Enumerating objects: 601, done.
remote: Counting objects: 100% (601/601), done.
remote: Compressing objects: 100% (333/333), done.
remote: Total 601 (delta 302), reused 513 (delta 217), pack-reused 0
Receiving objects: 100% (601/601), 3.93 MiB | 4.99 MiB/s, done.
Resolving deltas: 100% (302/302), done.
$ cd hsinhaoyu_cont_frac/src
$ git log -n1
commit 1d6f920 (HEAD -> main, origin/main, origin/HEAD)
Author: Hsin-Hao Yu [email protected]
Date: Sat Sep 3 11:23:36 2022 +1000

fixed typo

$ python3
Python 3.9.16 (main, Mar 8 2023, 22:47:22)
[GCC 11.3.0] on cygwin
Type "help", "copyright", "credits" or "license" for more information.

from itertools import islice
import cont_frac_io
def cf_sqrt5():
... yield 2
... while True:
... yield 4
...
for idx, n in enumerate(list(cont_frac_io.cf_convergents1(islice(cf_sqrt5(), 50)))):
... if idx >= 25:
... print(f'{idx:4d} {str(n):50s} ==> {n[0]/n[1]}')
...
25 10000136862780489/4472197161895732 ==> 2.2360679774997894
26 42361259535039638/18944531186571953 ==> 2.23606797749979
27 179445175002939041/80250321908183544 ==> 2.2360679774997894
28 760141959546795802/339945818819306129 ==> 2.23606797749979
29 3220013013190122249/1440033597185408060 ==> 2.23606797749979
30 -4806550061402266818/6100080207560938369 ==> -0.7879486658953493
31 2440556841290606593/7393610353719609920 ==> 0.33008999995013266
32 4955677303760159554/-1218966524979725183 ==> -4.065474483675905
33 3816521982621693193/2517744253800709188 ==> 1.515849744016053
34 1775021160537380710/8852010490223111569 ==> 0.20052180942373035
35 -7530137448938335583/1032298067274052232 ==> -7.294537970823548
36 8547959512203141610/-5465541314390231119 ==> -1.5639730852820026
37 8214956526164679241/-2383123116577320628 ==> -3.447138953510355
38 4514297469442755342/3448710293010037985 ==> 1.3089813541579547
39 7825402330226148993/-7035026018246720304 ==> -1.1123487404210635
40 -1077581357071751918/-6244649706267291615 ==> 0.1725607372324285
41 3515076901939141321/4879863304103216468 ==> 0.7203228211297437
42 -5464017823024738250/-5171940563563977359 ==> 1.0564734369761377
43 105749683549739937/2638845123556858648 ==> 0.04007422891389759
44 -5041019088825778502/5383439930663457233 ==> -0.9363936727728142
45 -1611582598043822455/5725860772501135964 ==> -0.2814568258074952
46 6959394592708483294/-8606605126751102143 ==> -0.8086108855020258
47 7779251699080559105/8192928412915830624 ==> 0.9495080766013864
48 1182913241611616482/5718364451202668737 ==> 0.20686216342206554
49 -5935839408182526583/-5827101929692597660 ==> 1.0186606446569686

Note that for iterations <= 29, a value close to the approximate correct value of 2.23606 is produced. Later iterations should produce more accurate results but instead produce values that are clearly incorrect. For instance, iteration 30 produced a value of -0.7879486658953493, which is clearly nowhere close to the approximate correct value of 2.23606. Note that in iteration 29 the numerator is equal to 3220013013190122249 (0x2CAFC8C66F551B09), which is less than 2^63. On the next iteration the numerator should grow to an even greater value, but instead becomes -4806550061402266818. It appears that during iteration 30, the numerator overflowed and became negative.

At first, this behavior was quite surprising in that by default Python automatically promotes an int to a long and then finally an arbitrary precision BigInteger. It appears for the issue demonstrated above that the issue is with the following statement:

res = np.matmul(res, h(a))

within the function cf_convergents1_ of cont_frac.py. It appears the issue is that the matmul function of numpy is unable to handle integer value >= 2^63.

It appears that this code makes significant use of the numpy module, which appears to be unable to promote integer values to an arbitrary precision type. Further, these overflow conditions are not detected/reported, which significantly limits the reliable use of this library.

use git hook to make sure that html and src are up to date

write README