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(<type 'list'>, {(19, 22, 7): [6, 7, 8, 9, 10, 11],
    (3, 3, 1): [3], (1843, 2130, 678): [16, 17]}),
    Brun 3-dimensional continued fraction algorithm:
    defaultdict(<type 'list'>, {(19, 22, 7): [9, 10, 11, 12, 13,
    14, 15, 16], (3, 3, 1): [4, 5], (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(<type 'list'>, {(19, 22, 7): [6, 7, 8, 9, 10, 11],
    (3, 3, 1): [3], (1843, 2130, 678): [16, 17]}),
    Brun 3-dimensional continued fraction algorithm:
    defaultdict(<type 'list'>, {(19, 22, 7): [9, 10, 11, 12, 13,
    14, 15, 16], (3, 3, 1): [4, 5], (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]