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 - vat distance- |qv-p|<=1/Qsuch 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 - vat 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]