# Diophantine Approximation¶

Simultaneous diophantine approximation

EXAMPLES:

sage: from slabbe.diophantine_approx import best_simultaneous_convergents
sage: it = best_simultaneous_convergents([e, pi])
sage: dirichlet10 = [next(it) for _ in range(10)]
sage: dirichlet10
[(3, 3, 1),
(19, 22, 7),
(1843, 2130, 678),
(51892, 59973, 19090),
(113018, 130618, 41577),
(114861, 132748, 42255),
(166753, 192721, 61345),
(446524, 516060, 164267),
(1174662, 1357589, 432134),
(3970510, 4588827, 1460669)]


The above Dirichlet simultaneous diophantine approximations appear as columns of matrices generated by multidimensional continued fraction algorithms:

sage: from slabbe.diophantine_approx import mult_cont_frac_vs_dirichlet
sage: from slabbe.mult_cont_frac import Brun,ARP,Reverse,Selmer,Cassaigne
sage: algos = [Brun(), ARP(), Reverse(), Selmer(),Cassaigne()]
sage: mult_cont_frac_vs_dirichlet([e,pi], dirichlet10, algos)
Dirichlet                     Brun       ARP        Reverse   Selmer     Cassaigne
+-----------------------------+----------+----------+---------+----------+-----------+
(3, 3, 1)                     [4, 5]     [3]        []        [8, 12]    []
(19, 22, 7)                   [9, 16]    [6, 11]    [7, 33]   [32]       [15, 25]
(1843, 2130, 678)             [22, 27]   [16, 17]   []        [44, 48]   [36, 39]
(51892, 59973, 19090)         []         []         []        [56]       []
(113018, 130618, 41577)       []         []         []        []         []
(114861, 132748, 42255)       [33, 35]   [22, 24]   []        [62, 66]   [51, 53]
(166753, 192721, 61345)       []         []         []        []         []
(446524, 516060, 164267)      [36]       [25]       []        []         [56, 57]
(1174662, 1357589, 432134)    [39, 44]   [26, 29]   []        [68]       [61, 66]
(3970510, 4588827, 1460669)   []         [28]       []        []         []


The indices in the table are the i-th matrices. For example, the first 3 Dirichlet approximations appear in the convergents of ARP algorithm, but not the 4th neither the 5th):

sage: algo = ARP()
sage: algo.n_matrix((e,pi,1), 3)
[3 3 2]
[3 4 2]
[1 1 1]
sage: algo.n_matrix((e,pi,1), 6)
[33 19  8]
[38 22  9]
[12  7  3]
sage: algo.n_matrix((e,pi,1), 16)
[1631 2498 1843]
[1885 2887 2130]
[ 600  919  678]
sage: algo.n_matrix((e,pi,1), 22)
[114861 101941  64812]
[132748 117816  74905]
[ 42255  37502  23843]
sage: algo.n_matrix((e,pi,1), 25)
[446524 331663 842999]
[516060 383312 974277]
[164267 122012 310122]
sage: algo.n_matrix((e,pi,1), 26)
[1621186  331663 1174662]
[1873649  383312 1357589]
[ 596401  122012  432134]
sage: algo.n_matrix((e,pi,1), 28)
[3970510 2680987 1174662]
[4588827 3098490 1357589]
[1460669  986280  432134]


The Dirichlet vectors can be computed from the precedent ones. So one could expect a perfect MCF algorithm based on 3x3 matrices:

sage: from slabbe.diophantine_approx import dirichlet_convergents_dependance
sage: dirichlet_convergents_dependance([e,pi], 8)
i   vi                         lin. rec.     remainder
+---+--------------------------+-------------+-----------+
0   (3, 3, 1)                  []            (3, 3, 1)
1   (19, 22, 7)                [6]           (1, 4, 1)
2   (1843, 2130, 678)          [96, 6]       (1, 0, 0)
3   (51892, 59973, 19090)      [28, 15, 1]   (0, 0, 0)
4   (113018, 130618, 41577)    [2, 5, 1]     (0, 0, 0)
5   (114861, 132748, 42255)    [1, 0, 1]     (0, 0, 0)
6   (166753, 192721, 61345)    [1, 0, 1]     (0, 0, 0)
7   (446524, 516060, 164267)   [2, 0, 1]     (0, 0, 0)


But, looking further, it is not true anymore, because v_10 is not a greedy linear combination of the previous three:

sage: dirichlet_convergents_dependance([e,pi], 18)   # long time (7.4s)
i    vi                                  lin. rec.         remainder
+----+-----------------------------------+-----------------+-----------+
0    (3, 3, 1)                           []                (3, 3, 1)
1    (19, 22, 7)                         [6]               (1, 4, 1)
2    (1843, 2130, 678)                   [96, 6]           (1, 0, 0)
3    (51892, 59973, 19090)               [28, 15, 1]       (0, 0, 0)
4    (113018, 130618, 41577)             [2, 5, 1]         (0, 0, 0)
5    (114861, 132748, 42255)             [1, 0, 1]         (0, 0, 0)
6    (166753, 192721, 61345)             [1, 0, 1]         (0, 0, 0)
7    (446524, 516060, 164267)            [2, 0, 1]         (0, 0, 0)
8    (1174662, 1357589, 432134)          [2, 1, 1]         (0, 0, 0)
9    (3970510, 4588827, 1460669)         [3, 1]            (0, 0, 0)
10   (21640489, 25010505, 7961091)       [5, 1, 1, 1]      (0, 0, 0)
11   (25610999, 29599332, 9421760)       [1, 1]            (0, 0, 0)
12   (47251488, 54609837, 17382851)      [1, 1]            (0, 0, 0)
13   (117318127, 135587768, 43158927)    [2, 0, 1, 0, 1]   (0, 0, 0)
14   (142929126, 165187100, 52580687)    [1, 0, 1]         (0, 0, 0)
15   (164569615, 190197605, 60541778)    [1, 0, 0, 0, 1]   (0, 0, 0)
16   (307498741, 355384705, 113122465)   [1, 1]            (0, 0, 0)
17   (779567097, 900967015, 286786708)   [2, 1]            (0, 0, 0)


Remark: matrices of independant lines is 1 (with Thomas Garrity, 27oct 2016). He kind of have a proof of that…

This used to give a Value Error:

sage: dirichlet_convergents_dependance([e,pi], 19) # not tested (1min 12s)
i    vi                                     lin. rec.             remainder
+----+--------------------------------------+---------------------+-----------+
0    (3, 3, 1)                              []                    (3, 3, 1)
1    (19, 22, 7)                            [6]                   (1, 4, 1)
2    (1843, 2130, 678)                      [96, 6]               (1, 0, 0)
3    (51892, 59973, 19090)                  [28, 15, 1]           (0, 0, 0)
4    (113018, 130618, 41577)                [2, 5, 1]             (0, 0, 0)
5    (114861, 132748, 42255)                [1, 0, 1]             (0, 0, 0)
6    (166753, 192721, 61345)                [1, 0, 1]             (0, 0, 0)
7    (446524, 516060, 164267)               [2, 0, 1]             (0, 0, 0)
8    (1174662, 1357589, 432134)             [2, 1, 1]             (0, 0, 0)
9    (3970510, 4588827, 1460669)            [3, 1]                (0, 0, 0)
10   (21640489, 25010505, 7961091)          [5, 1, 1, 1]          (0, 0, 0)
11   (25610999, 29599332, 9421760)          [1, 1]                (0, 0, 0)
12   (47251488, 54609837, 17382851)         [1, 1]                (0, 0, 0)
13   (117318127, 135587768, 43158927)       [2, 0, 1, 0, 1]       (0, 0, 0)
14   (142929126, 165187100, 52580687)       [1, 0, 1]             (0, 0, 0)
15   (164569615, 190197605, 60541778)       [1, 0, 0, 0, 1]       (0, 0, 0)
16   (307498741, 355384705, 113122465)      [1, 1]                (0, 0, 0)
17   (779567097, 900967015, 286786708)      [2, 1]                (0, 0, 0)
18   (8457919940, 9775049397, 3111494861)   [10, 2, 0, 0, 0, 1]   (0, 0, 0)


BENCHMARKS:

sage: it = best_simultaneous_convergents([e,pi])
sage: %time L = [next(it) for _ in range(15)]   # not tested (4.66s)
sage: it = best_simultaneous_convergents([e,pi])
sage: %time L = [next(it) for _ in range(18)]   # not tested (52s)


AUTHORS:

• Sébastien Labbé, September 22, 2016
• Sébastien Labbé, October 10, 2016
• Sébastien Labbé, October 19, 2016, Cython improvements (10000x faster)
• Sébastien Labbé, October 21, 2016, Parallelization of computations

TODO:

• In general, how many of the dirichlet approximations are found by the MCF algos?
• Code Szekeres MCF algorithm

REFERENCES:

Szekeres, G. Multidimensional continued fractions. Ann. Univ. Sci. Budapest. Eötvös Sect. Math. 13 (1970), 113–140 (1971).

Cusick, T. W. The Szekeres multidimensional continued fraction. Math. Comp. 31 (1977), no. 137, 280–317.

slabbe.diophantine_approx.best_simultaneous_convergents(v)

Return an iterator of best convergents to a vector of real number according to Dirichlet theorem on simultaneous approximations.

INPUT:

• v – list of real numbers

OUTPUT:

• iterator

EXAMPLES:

sage: from slabbe.diophantine_approx import best_simultaneous_convergents
sage: it = best_simultaneous_convergents([e, pi])
sage: [next(it) for _ in range(5)]
[(3, 3, 1),
(19, 22, 7),
(1843, 2130, 678),
(51892, 59973, 19090),
(113018, 130618, 41577)]


TESTS:

Correspondance with continued fraction when d=1:

sage: it = best_simultaneous_convergents([e])
sage: [next(it) for _ in range(10)]
[(3, 1),
(8, 3),
(11, 4),
(19, 7),
(87, 32),
(106, 39),
(193, 71),
(1264, 465),
(1457, 536),
(2721, 1001)]
sage: continued_fraction(e).convergents()[:11].list()
[2, 3, 8/3, 11/4, 19/7, 87/32, 106/39, 193/71, 1264/465, 1457/536, 2721/1001]

slabbe.diophantine_approx.best_simultaneous_convergents_upto(v, Q, start=1, verbose=False)

Return a list of all best simultaneous diophantine approximations p,q of vector v at distance |qv-p|<=1/Q such that $$1<=q<Q^d$$.

INPUT:

• v – list of real numbers
• Q – real number, Q>1
• start – integer (default: 1), starting value to check
• verbose – boolean (default: False)

EXAMPLES:

sage: from slabbe.diophantine_approx import best_simultaneous_convergents_upto
sage: best_simultaneous_convergents_upto([e,pi], 2)
[((3, 3, 1), 3.549646778303845)]
sage: best_simultaneous_convergents_upto([e,pi], 4)
[((19, 22, 7), 35.74901433260719)]
sage: best_simultaneous_convergents_upto([e,pi], 36, start=4**2)
[((1843, 2130, 678), 203.23944293852406)]
sage: best_simultaneous_convergents_upto([e,pi], 204, start=36**2)
[((51892, 59973, 19090), 266.16775098010373),
((113018, 130618, 41577), 279.18598227531174)]
sage: best_simultaneous_convergents_upto([e,pi], 280, start=204**2)
[((114861, 132748, 42255), 412.7859137824949),
((166753, 192721, 61345), 749.3634909055199)]
sage: best_simultaneous_convergents_upto([e,pi], 750, start=280**2)
[((446524, 516060, 164267), 896.4734658499202),
((1174662, 1357589, 432134), 2935.314937919726)]
sage: best_simultaneous_convergents_upto([e,pi], 2936, start=750**2)
[((3970510, 4588827, 1460669), 3654.2989332956854),
((21640489, 25010505, 7961091), 6257.014011585661)]


TESTS:

sage: best_simultaneous_convergents_upto([e,pi], 102300.1, start=10^9) # not tested (1 min)
[((8457919940, 9775049397, 3111494861), 194686.19839453633)]

slabbe.diophantine_approx.dirichlet_convergents_dependance(v, n, verbose=False)

INPUT:

• v – list of real numbers
• n – integer, number of iterations
• verbose – bool (default: False),

OUTPUT:

• table of linear combinaisons of dirichlet approximations in terms of previous dirichlet approximations

EXAMPLES:

sage: from slabbe.diophantine_approx import dirichlet_convergents_dependance
sage: dirichlet_convergents_dependance([e,pi], 4)
i   vi                      lin. rec.     remainder
+---+-----------------------+-------------+-----------+
0   (3, 3, 1)               []            (3, 3, 1)
1   (19, 22, 7)             [6]           (1, 4, 1)
2   (1843, 2130, 678)       [96, 6]       (1, 0, 0)
3   (51892, 59973, 19090)   [28, 15, 1]   (0, 0, 0)


The last 3 seems enough:

sage: dirichlet_convergents_dependance([e,pi], 8)
i   vi                         lin. rec.     remainder
+---+--------------------------+-------------+-----------+
0   (3, 3, 1)                  []            (3, 3, 1)
1   (19, 22, 7)                [6]           (1, 4, 1)
2   (1843, 2130, 678)          [96, 6]       (1, 0, 0)
3   (51892, 59973, 19090)      [28, 15, 1]   (0, 0, 0)
4   (113018, 130618, 41577)    [2, 5, 1]     (0, 0, 0)
5   (114861, 132748, 42255)    [1, 0, 1]     (0, 0, 0)
6   (166753, 192721, 61345)    [1, 0, 1]     (0, 0, 0)
7   (446524, 516060, 164267)   [2, 0, 1]     (0, 0, 0)


But not in this case:

sage: dirichlet_convergents_dependance([pi,sqrt(3)], 12)
i    vi                      lin. rec.                  remainder
+----+-----------------------+--------------------------+-----------+
0    (3, 2, 1)               []                         (3, 2, 1)
1    (22, 12, 7)             [6]                        (4, 0, 1)
2    (47, 26, 15)            [2, 1]                     (0, 0, 0)
3    (69, 38, 22)            [1, 1]                     (0, 0, 0)
4    (176, 97, 56)           [2, 0, 1, 4]               (4, 1, 1)
5    (223, 123, 71)          [1, 0, 1]                  (0, 0, 0)
6    (399, 220, 127)         [1, 1]                     (0, 0, 0)
7    (1442, 795, 459)        [3, 1, 0, 0, 0, 1]         (0, 0, 0)
8    (6390, 3523, 2034)      [4, 1, 1]                  (0, 0, 0)
9    (26603, 14667, 8468)    [4, 0, 2, 1, 0, 0, 0, 1]   (0, 0, 0)
10   (32993, 18190, 10502)   [1, 1]                     (0, 0, 0)
11   (40825, 22508, 12995)   [1, 0, 1, 1]               (0, 0, 0)


The v4 is not a lin. comb. of the previous four:

sage: dirichlet_convergents_dependance([e,pi,sqrt(3)], 5)
i   vi                                lin. rec.        remainder
+---+---------------------------------+----------------+--------------+
0   (3, 3, 2, 1)                      []               (3, 3, 2, 1)
1   (19, 22, 12, 7)                   [6]              (1, 4, 0, 1)
2   (193, 223, 123, 71)               [10, 1]          (0, 0, 1, 0)
3   (5529, 6390, 3523, 2034)          [28, 6, 3]       (2, 5, 1, 1)
4   (163067, 188461, 103904, 59989)   [29, 14, 1, 1]   (2, 4, 1, 1)

slabbe.diophantine_approx.distance_to_nearest_integer(x)
slabbe.diophantine_approx.frac(x)

Return the fractional part of real number x.

Not always perfect…

EXAMPLES:

sage: from slabbe.diophantine_approx import frac
sage: frac(3.2)
0.200000000000000
sage: frac(-3.2)
0.800000000000000
sage: frac(pi)
pi - 3
sage: frac(pi).n()
0.141592653589793


This looks suspicious…:

sage: frac(pi*10**15).n()
0.000000000000000

slabbe.diophantine_approx.mult_cont_frac_vs_dirichlet(v, dirichlet, algos)

Returns the indices i such that dirichlet approximations appears as columns of the i-th matrix obtained from mult. dim. cont. frac. algorithms.

INPUT:

• v – list of real numbers
• dirichlet – list, first dirichlet approximations
• algos – list, list of mult. cont. frac. algorithms

OUTPUT:

table

EXAMPLES:

sage: from slabbe.diophantine_approx import best_simultaneous_convergents
sage: from slabbe.diophantine_approx import mult_cont_frac_vs_dirichlet
sage: from slabbe.mult_cont_frac import Brun,ARP,Reverse,Selmer,Cassaigne
sage: v = [e, pi]
sage: it = best_simultaneous_convergents(v)
sage: dirichlet = [next(it) for _ in range(10)]
sage: algos = [Brun(), ARP(), Reverse(), Selmer(),Cassaigne()]
sage: mult_cont_frac_vs_dirichlet(v, dirichlet, algos)
Dirichlet                     Brun       ARP        Reverse   Selmer     Cassaigne
+-----------------------------+----------+----------+---------+----------+-----------+
(3, 3, 1)                     [4, 5]     [3]        []        [8, 12]    []
(19, 22, 7)                   [9, 16]    [6, 11]    [7, 33]   [32]       [15, 25]
(1843, 2130, 678)             [22, 27]   [16, 17]   []        [44, 48]   [36, 39]
(51892, 59973, 19090)         []         []         []        [56]       []
(113018, 130618, 41577)       []         []         []        []         []
(114861, 132748, 42255)       [33, 35]   [22, 24]   []        [62, 66]   [51, 53]
(166753, 192721, 61345)       []         []         []        []         []
(446524, 516060, 164267)      [36]       [25]       []        []         [56, 57]
(1174662, 1357589, 432134)    [39, 44]   [26, 29]   []        [68]       [61, 66]
(3970510, 4588827, 1460669)   []         [28]       []        []         []

slabbe.diophantine_approx.mult_cont_frac_vs_dirichlet_dict(v, dirichlet, algos)

INPUT:

• v – list of real numbers
• dirichlet – list, first dirichlet approximations
• algos – list, list of mult. cont. frac. algorithms

OUTPUT:

dict

EXAMPLES:

   sage: from slabbe.diophantine_approx import best_simultaneous_convergents
sage: from slabbe.diophantine_approx import mult_cont_frac_vs_dirichlet_dict
sage: from slabbe.mult_cont_frac import ARP, Brun
sage: v = [e, pi]
sage: it = best_simultaneous_convergents(v)
sage: dirichlet = [next(it) for _ in range(3)]
sage: mult_cont_frac_vs_dirichlet_dict([e,pi], dirichlet, [Brun(), ARP()])
{Arnoux-Rauzy-Poincar\'e 3-dimensional continued fraction algorithm: defaultdict(<class 'list'>, {(3, 3, 1): [3], (19, 22, 7): [6, 7, 8, 9, 10, 11], (1843, 2130, 678): [16, 17]}),
Brun 3-dimensional continued fraction algorithm: defaultdict(<class 'list'>, {(3, 3, 1): [4, 5], (19, 22, 7): [9, 10, 11, 12, 13, 14, 15, 16], (1843, 2130, 678): [22, 23, 24, 25, 26, 27]})}

Or from precomputed Dirichlet approximations::

sage: dirichlet = [(3, 3, 1), (19, 22, 7), (1843, 2130, 678), (51892, 59973, 19090)]
sage: mult_cont_frac_vs_dirichlet_dict([e,pi], dirichlet, [Brun(), ARP()])
{Arnoux-Rauzy-Poincar\'e 3-dimensional continued fraction algorithm: defaultdict(<class 'list'>, {(3, 3, 1): [3], (19, 22, 7): [6, 7, 8, 9, 10, 11], (1843, 2130, 678): [16, 17]}),
Brun 3-dimensional continued fraction algorithm: defaultdict(<class 'list'>, {(3, 3, 1): [4, 5], (19, 22, 7): [9, 10, 11, 12, 13, 14, 15, 16], (1843, 2130, 678): [22, 23, 24, 25, 26, 27]})}


Simultaneous diophantine approximation

EXAMPLES:

The code gets between 1000x faster and 50000x faster compared to the same function in Python code (diophantine_approximation.py):

sage: from slabbe.diophantine_approx_pyx import good_simultaneous_convergents_upto
sage: from slabbe.diophantine_approx import _best_simultaneous_convergents_upto
sage: good_simultaneous_convergents_upto([e,pi], 203)     # 493 µs
[678, 19090, 19768]
sage: _best_simultaneous_convergents_upto([e,pi], 203)         # 905 ms
((1843, 2130, 678), 203.239442934072)
sage: 905/493. * 1000
1835.69979716024

sage: good_simultaneous_convergents_upto([e,pi], 204)    # 493 µs
[19090, 19768, 41577]
sage: _best_simultaneous_convergents_upto([e,pi], 204)         # 25s (not tested)
((51892, 59973, 19090), 266.17)
sage: 25 / 493. * 1000000
50709.9391480730


AUTHORS:

• Sébastien Labbé, October 19, 2016
slabbe.diophantine_approx_pyx.good_simultaneous_convergents_upto()

Return a list of potential best simultaneous diophantine approximation of vector v at distance 1/Q.

It searches for all possibilities of denominators in the interval [1, Q^d] starting at start by step. Argument step is used for parallel computations purposes.

INPUT:

• v – list of real numbers to approximate
• Q – real number, Q>1
• start – integer (default: 1), starting value to check
• step – integer (default: 1), step

OUTPUT:

• list of integers q such that entries of q*v are at most 1/Q from the integer lattice

EXAMPLES:

sage: from slabbe.diophantine_approx_pyx import good_simultaneous_convergents_upto
sage: good_simultaneous_convergents_upto([e,pi], 2)
[1, 2, 3, 4]
sage: good_simultaneous_convergents_upto([e,pi], 4)
[7, 14, 15]
sage: good_simultaneous_convergents_upto([e,pi], 35)
[7, 678, 685]
sage: good_simultaneous_convergents_upto([e,pi], 36)
[678, 685]
sage: good_simultaneous_convergents_upto([e,pi], 203)
[678, 19090, 19768]
sage: good_simultaneous_convergents_upto([e,pi], 204)
[19090, 19768, 41577]


The previous results mean that these integers are good approximations:

sage: v = [e,pi]
sage: [(7 * a).n() for a in v]
[19.0279727992133, 21.9911485751286]
sage: [(678 * a).n() for a in v]
[1842.99507969523, 2129.99981913388]
sage: [(19090 * a).n() for a in v]
[51892.0001052832, 59973.0037570292]


Knowing the error of previous results, we can start the next computation at step 678:

sage: min(1/abs(678*a - p) for a,p in zip([e,pi], (1843, 2130))).n()
203.239442934072
sage: good_simultaneous_convergents_upto([e,pi], 204, start=678)
[19090, 19768, 41577]


Use start and step for parallel computation purposes:

sage: good_simultaneous_convergents_upto([e,pi], 2935, start=1)
[432134, 1460669, 7961091, 8393225]
sage: good_simultaneous_convergents_upto([e,pi], 2935, start=1, step=3)
[]
sage: good_simultaneous_convergents_upto([e,pi], 2935, start=2, step=3)
[432134, 1460669, 8393225]
sage: good_simultaneous_convergents_upto([e,pi], 2935, start=3, step=3)
[7961091]


TESTS:

sage: good_simultaneous_convergents_upto([1,e,pi], 1)
Traceback (most recent call last):
...
ValueError: argument Q(=1.0) must be > 1


If the interval is too short, then noting is found:

sage: good_simultaneous_convergents_upto([e,pi], 204, start=42000, step=2)
[]


No solution was found when the q was $$cdef int$$ instead of $$cdef long$$:

sage: good_simultaneous_convergents_upto([e,pi], 136500, start=1)  # not tested (2min)
[3111494861]


This used to give a Value Error:

sage: good_simultaneous_convergents_upto([e,pi], 102300, start=10^9, step=2) # not tested (1min)
[]
sage: good_simultaneous_convergents_upto([e,pi], 102300, start=10^9+1, step=2) # not tested (1min)
[3111494861, 3398281569]