Coverage for /usr/lib/python3/dist-packages/sympy/ntheory/egyptian_fraction.py: 13%
87 statements
« prev ^ index » next coverage.py v7.9.1, created at 2025-06-14 15:55 +0200
1from sympy.core.containers import Tuple
2from sympy.core.numbers import (Integer, Rational)
3from sympy.core.singleton import S
4import sympy.polys
6from math import gcd
9def egyptian_fraction(r, algorithm="Greedy"):
10 """
11 Return the list of denominators of an Egyptian fraction
12 expansion [1]_ of the said rational `r`.
14 Parameters
15 ==========
17 r : Rational or (p, q)
18 a positive rational number, ``p/q``.
19 algorithm : { "Greedy", "Graham Jewett", "Takenouchi", "Golomb" }, optional
20 Denotes the algorithm to be used (the default is "Greedy").
22 Examples
23 ========
25 >>> from sympy import Rational
26 >>> from sympy.ntheory.egyptian_fraction import egyptian_fraction
27 >>> egyptian_fraction(Rational(3, 7))
28 [3, 11, 231]
29 >>> egyptian_fraction((3, 7), "Graham Jewett")
30 [7, 8, 9, 56, 57, 72, 3192]
31 >>> egyptian_fraction((3, 7), "Takenouchi")
32 [4, 7, 28]
33 >>> egyptian_fraction((3, 7), "Golomb")
34 [3, 15, 35]
35 >>> egyptian_fraction((11, 5), "Golomb")
36 [1, 2, 3, 4, 9, 234, 1118, 2580]
38 See Also
39 ========
41 sympy.core.numbers.Rational
43 Notes
44 =====
46 Currently the following algorithms are supported:
48 1) Greedy Algorithm
50 Also called the Fibonacci-Sylvester algorithm [2]_.
51 At each step, extract the largest unit fraction less
52 than the target and replace the target with the remainder.
54 It has some distinct properties:
56 a) Given `p/q` in lowest terms, generates an expansion of maximum
57 length `p`. Even as the numerators get large, the number of
58 terms is seldom more than a handful.
60 b) Uses minimal memory.
62 c) The terms can blow up (standard examples of this are 5/121 and
63 31/311). The denominator is at most squared at each step
64 (doubly-exponential growth) and typically exhibits
65 singly-exponential growth.
67 2) Graham Jewett Algorithm
69 The algorithm suggested by the result of Graham and Jewett.
70 Note that this has a tendency to blow up: the length of the
71 resulting expansion is always ``2**(x/gcd(x, y)) - 1``. See [3]_.
73 3) Takenouchi Algorithm
75 The algorithm suggested by Takenouchi (1921).
76 Differs from the Graham-Jewett algorithm only in the handling
77 of duplicates. See [3]_.
79 4) Golomb's Algorithm
81 A method given by Golumb (1962), using modular arithmetic and
82 inverses. It yields the same results as a method using continued
83 fractions proposed by Bleicher (1972). See [4]_.
85 If the given rational is greater than or equal to 1, a greedy algorithm
86 of summing the harmonic sequence 1/1 + 1/2 + 1/3 + ... is used, taking
87 all the unit fractions of this sequence until adding one more would be
88 greater than the given number. This list of denominators is prefixed
89 to the result from the requested algorithm used on the remainder. For
90 example, if r is 8/3, using the Greedy algorithm, we get [1, 2, 3, 4,
91 5, 6, 7, 14, 420], where the beginning of the sequence, [1, 2, 3, 4, 5,
92 6, 7] is part of the harmonic sequence summing to 363/140, leaving a
93 remainder of 31/420, which yields [14, 420] by the Greedy algorithm.
94 The result of egyptian_fraction(Rational(8, 3), "Golomb") is [1, 2, 3,
95 4, 5, 6, 7, 14, 574, 2788, 6460, 11590, 33062, 113820], and so on.
97 References
98 ==========
100 .. [1] https://en.wikipedia.org/wiki/Egyptian_fraction
101 .. [2] https://en.wikipedia.org/wiki/Greedy_algorithm_for_Egyptian_fractions
102 .. [3] https://www.ics.uci.edu/~eppstein/numth/egypt/conflict.html
103 .. [4] https://web.archive.org/web/20180413004012/https://ami.ektf.hu/uploads/papers/finalpdf/AMI_42_from129to134.pdf
105 """
107 if not isinstance(r, Rational):
108 if isinstance(r, (Tuple, tuple)) and len(r) == 2:
109 r = Rational(*r)
110 else:
111 raise ValueError("Value must be a Rational or tuple of ints")
112 if r <= 0:
113 raise ValueError("Value must be positive")
115 # common cases that all methods agree on
116 x, y = r.as_numer_denom()
117 if y == 1 and x == 2:
118 return [Integer(i) for i in [1, 2, 3, 6]]
119 if x == y + 1:
120 return [S.One, y]
122 prefix, rem = egypt_harmonic(r)
123 if rem == 0:
124 return prefix
125 # work in Python ints
126 x, y = rem.p, rem.q
127 # assert x < y and gcd(x, y) = 1
129 if algorithm == "Greedy":
130 postfix = egypt_greedy(x, y)
131 elif algorithm == "Graham Jewett":
132 postfix = egypt_graham_jewett(x, y)
133 elif algorithm == "Takenouchi":
134 postfix = egypt_takenouchi(x, y)
135 elif algorithm == "Golomb":
136 postfix = egypt_golomb(x, y)
137 else:
138 raise ValueError("Entered invalid algorithm")
139 return prefix + [Integer(i) for i in postfix]
142def egypt_greedy(x, y):
143 # assumes gcd(x, y) == 1
144 if x == 1:
145 return [y]
146 else:
147 a = (-y) % x
148 b = y*(y//x + 1)
149 c = gcd(a, b)
150 if c > 1:
151 num, denom = a//c, b//c
152 else:
153 num, denom = a, b
154 return [y//x + 1] + egypt_greedy(num, denom)
157def egypt_graham_jewett(x, y):
158 # assumes gcd(x, y) == 1
159 l = [y] * x
161 # l is now a list of integers whose reciprocals sum to x/y.
162 # we shall now proceed to manipulate the elements of l without
163 # changing the reciprocated sum until all elements are unique.
165 while len(l) != len(set(l)):
166 l.sort() # so the list has duplicates. find a smallest pair
167 for i in range(len(l) - 1):
168 if l[i] == l[i + 1]:
169 break
170 # we have now identified a pair of identical
171 # elements: l[i] and l[i + 1].
172 # now comes the application of the result of graham and jewett:
173 l[i + 1] = l[i] + 1
174 # and we just iterate that until the list has no duplicates.
175 l.append(l[i]*(l[i] + 1))
176 return sorted(l)
179def egypt_takenouchi(x, y):
180 # assumes gcd(x, y) == 1
181 # special cases for 3/y
182 if x == 3:
183 if y % 2 == 0:
184 return [y//2, y]
185 i = (y - 1)//2
186 j = i + 1
187 k = j + i
188 return [j, k, j*k]
189 l = [y] * x
190 while len(l) != len(set(l)):
191 l.sort()
192 for i in range(len(l) - 1):
193 if l[i] == l[i + 1]:
194 break
195 k = l[i]
196 if k % 2 == 0:
197 l[i] = l[i] // 2
198 del l[i + 1]
199 else:
200 l[i], l[i + 1] = (k + 1)//2, k*(k + 1)//2
201 return sorted(l)
204def egypt_golomb(x, y):
205 # assumes x < y and gcd(x, y) == 1
206 if x == 1:
207 return [y]
208 xp = sympy.polys.ZZ.invert(int(x), int(y))
209 rv = [xp*y]
210 rv.extend(egypt_golomb((x*xp - 1)//y, xp))
211 return sorted(rv)
214def egypt_harmonic(r):
215 # assumes r is Rational
216 rv = []
217 d = S.One
218 acc = S.Zero
219 while acc + 1/d <= r:
220 acc += 1/d
221 rv.append(d)
222 d += 1
223 return (rv, r - acc)