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 numbersQ
– real number, Q>1start
– integer (default:1
), starting value to checkverbose
– 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 numbersn
– integer, number of iterationsverbose
– 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 numbersdirichlet
– list, first dirichlet approximationsalgos
– 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 numbersdirichlet
– list, first dirichlet approximationsalgos
– list, list of mult. cont. frac. algorithms
OUTPUT:
dictEXAMPLES:
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 distance1/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 approximateQ
– real number, Q>1start
– integer (default:1
), starting value to checkstep
– 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]