continued_fraction.py 10 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369
  1. from __future__ import annotations
  2. import itertools
  3. from sympy.core.exprtools import factor_terms
  4. from sympy.core.numbers import Integer, Rational
  5. from sympy.core.singleton import S
  6. from sympy.core.symbol import Dummy
  7. from sympy.core.sympify import _sympify
  8. from sympy.utilities.misc import as_int
  9. def continued_fraction(a) -> list:
  10. """Return the continued fraction representation of a Rational or
  11. quadratic irrational.
  12. Examples
  13. ========
  14. >>> from sympy.ntheory.continued_fraction import continued_fraction
  15. >>> from sympy import sqrt
  16. >>> continued_fraction((1 + 2*sqrt(3))/5)
  17. [0, 1, [8, 3, 34, 3]]
  18. See Also
  19. ========
  20. continued_fraction_periodic, continued_fraction_reduce, continued_fraction_convergents
  21. """
  22. e = _sympify(a)
  23. if all(i.is_Rational for i in e.atoms()):
  24. if e.is_Integer:
  25. return continued_fraction_periodic(e, 1, 0)
  26. elif e.is_Rational:
  27. return continued_fraction_periodic(e.p, e.q, 0)
  28. elif e.is_Pow and e.exp is S.Half and e.base.is_Integer:
  29. return continued_fraction_periodic(0, 1, e.base)
  30. elif e.is_Mul and len(e.args) == 2 and (
  31. e.args[0].is_Rational and
  32. e.args[1].is_Pow and
  33. e.args[1].base.is_Integer and
  34. e.args[1].exp is S.Half):
  35. a, b = e.args
  36. return continued_fraction_periodic(0, a.q, b.base, a.p)
  37. else:
  38. # this should not have to work very hard- no
  39. # simplification, cancel, etc... which should be
  40. # done by the user. e.g. This is a fancy 1 but
  41. # the user should simplify it first:
  42. # sqrt(2)*(1 + sqrt(2))/(sqrt(2) + 2)
  43. p, d = e.expand().as_numer_denom()
  44. if d.is_Integer:
  45. if p.is_Rational:
  46. return continued_fraction_periodic(p, d)
  47. # look for a + b*c
  48. # with c = sqrt(s)
  49. if p.is_Add and len(p.args) == 2:
  50. a, bc = p.args
  51. else:
  52. a = S.Zero
  53. bc = p
  54. if a.is_Integer:
  55. b = S.NaN
  56. if bc.is_Mul and len(bc.args) == 2:
  57. b, c = bc.args
  58. elif bc.is_Pow:
  59. b = Integer(1)
  60. c = bc
  61. if b.is_Integer and (
  62. c.is_Pow and c.exp is S.Half and
  63. c.base.is_Integer):
  64. # (a + b*sqrt(c))/d
  65. c = c.base
  66. return continued_fraction_periodic(a, d, c, b)
  67. raise ValueError(
  68. 'expecting a rational or quadratic irrational, not %s' % e)
  69. def continued_fraction_periodic(p, q, d=0, s=1) -> list:
  70. r"""
  71. Find the periodic continued fraction expansion of a quadratic irrational.
  72. Compute the continued fraction expansion of a rational or a
  73. quadratic irrational number, i.e. `\frac{p + s\sqrt{d}}{q}`, where
  74. `p`, `q \ne 0` and `d \ge 0` are integers.
  75. Returns the continued fraction representation (canonical form) as
  76. a list of integers, optionally ending (for quadratic irrationals)
  77. with list of integers representing the repeating digits.
  78. Parameters
  79. ==========
  80. p : int
  81. the rational part of the number's numerator
  82. q : int
  83. the denominator of the number
  84. d : int, optional
  85. the irrational part (discriminator) of the number's numerator
  86. s : int, optional
  87. the coefficient of the irrational part
  88. Examples
  89. ========
  90. >>> from sympy.ntheory.continued_fraction import continued_fraction_periodic
  91. >>> continued_fraction_periodic(3, 2, 7)
  92. [2, [1, 4, 1, 1]]
  93. Golden ratio has the simplest continued fraction expansion:
  94. >>> continued_fraction_periodic(1, 2, 5)
  95. [[1]]
  96. If the discriminator is zero or a perfect square then the number will be a
  97. rational number:
  98. >>> continued_fraction_periodic(4, 3, 0)
  99. [1, 3]
  100. >>> continued_fraction_periodic(4, 3, 49)
  101. [3, 1, 2]
  102. See Also
  103. ========
  104. continued_fraction_iterator, continued_fraction_reduce
  105. References
  106. ==========
  107. .. [1] https://en.wikipedia.org/wiki/Periodic_continued_fraction
  108. .. [2] K. Rosen. Elementary Number theory and its applications.
  109. Addison-Wesley, 3 Sub edition, pages 379-381, January 1992.
  110. """
  111. from sympy.functions import sqrt, floor
  112. p, q, d, s = list(map(as_int, [p, q, d, s]))
  113. if d < 0:
  114. raise ValueError("expected non-negative for `d` but got %s" % d)
  115. if q == 0:
  116. raise ValueError("The denominator cannot be 0.")
  117. if not s:
  118. d = 0
  119. # check for rational case
  120. sd = sqrt(d)
  121. if sd.is_Integer:
  122. return list(continued_fraction_iterator(Rational(p + s*sd, q)))
  123. # irrational case with sd != Integer
  124. if q < 0:
  125. p, q, s = -p, -q, -s
  126. n = (p + s*sd)/q
  127. if n < 0:
  128. w = floor(-n)
  129. f = -n - w
  130. one_f = continued_fraction(1 - f) # 1-f < 1 so cf is [0 ... [...]]
  131. one_f[0] -= w + 1
  132. return one_f
  133. d *= s**2
  134. sd *= s
  135. if (d - p**2)%q:
  136. d *= q**2
  137. sd *= q
  138. p *= q
  139. q *= q
  140. terms: list[int] = []
  141. pq = {}
  142. while (p, q) not in pq:
  143. pq[(p, q)] = len(terms)
  144. terms.append((p + sd)//q)
  145. p = terms[-1]*q - p
  146. q = (d - p**2)//q
  147. i = pq[(p, q)]
  148. return terms[:i] + [terms[i:]] # type: ignore
  149. def continued_fraction_reduce(cf):
  150. """
  151. Reduce a continued fraction to a rational or quadratic irrational.
  152. Compute the rational or quadratic irrational number from its
  153. terminating or periodic continued fraction expansion. The
  154. continued fraction expansion (cf) should be supplied as a
  155. terminating iterator supplying the terms of the expansion. For
  156. terminating continued fractions, this is equivalent to
  157. ``list(continued_fraction_convergents(cf))[-1]``, only a little more
  158. efficient. If the expansion has a repeating part, a list of the
  159. repeating terms should be returned as the last element from the
  160. iterator. This is the format returned by
  161. continued_fraction_periodic.
  162. For quadratic irrationals, returns the largest solution found,
  163. which is generally the one sought, if the fraction is in canonical
  164. form (all terms positive except possibly the first).
  165. Examples
  166. ========
  167. >>> from sympy.ntheory.continued_fraction import continued_fraction_reduce
  168. >>> continued_fraction_reduce([1, 2, 3, 4, 5])
  169. 225/157
  170. >>> continued_fraction_reduce([-2, 1, 9, 7, 1, 2])
  171. -256/233
  172. >>> continued_fraction_reduce([2, 1, 2, 1, 1, 4, 1, 1, 6, 1, 1, 8]).n(10)
  173. 2.718281835
  174. >>> continued_fraction_reduce([1, 4, 2, [3, 1]])
  175. (sqrt(21) + 287)/238
  176. >>> continued_fraction_reduce([[1]])
  177. (1 + sqrt(5))/2
  178. >>> from sympy.ntheory.continued_fraction import continued_fraction_periodic
  179. >>> continued_fraction_reduce(continued_fraction_periodic(8, 5, 13))
  180. (sqrt(13) + 8)/5
  181. See Also
  182. ========
  183. continued_fraction_periodic
  184. """
  185. from sympy.solvers import solve
  186. period = []
  187. x = Dummy('x')
  188. def untillist(cf):
  189. for nxt in cf:
  190. if isinstance(nxt, list):
  191. period.extend(nxt)
  192. yield x
  193. break
  194. yield nxt
  195. a = S.Zero
  196. for a in continued_fraction_convergents(untillist(cf)):
  197. pass
  198. if period:
  199. y = Dummy('y')
  200. solns = solve(continued_fraction_reduce(period + [y]) - y, y)
  201. solns.sort()
  202. pure = solns[-1]
  203. rv = a.subs(x, pure).radsimp()
  204. else:
  205. rv = a
  206. if rv.is_Add:
  207. rv = factor_terms(rv)
  208. if rv.is_Mul and rv.args[0] == -1:
  209. rv = rv.func(*rv.args)
  210. return rv
  211. def continued_fraction_iterator(x):
  212. """
  213. Return continued fraction expansion of x as iterator.
  214. Examples
  215. ========
  216. >>> from sympy import Rational, pi
  217. >>> from sympy.ntheory.continued_fraction import continued_fraction_iterator
  218. >>> list(continued_fraction_iterator(Rational(3, 8)))
  219. [0, 2, 1, 2]
  220. >>> list(continued_fraction_iterator(Rational(-3, 8)))
  221. [-1, 1, 1, 1, 2]
  222. >>> for i, v in enumerate(continued_fraction_iterator(pi)):
  223. ... if i > 7:
  224. ... break
  225. ... print(v)
  226. 3
  227. 7
  228. 15
  229. 1
  230. 292
  231. 1
  232. 1
  233. 1
  234. References
  235. ==========
  236. .. [1] https://en.wikipedia.org/wiki/Continued_fraction
  237. """
  238. from sympy.functions import floor
  239. while True:
  240. i = floor(x)
  241. yield i
  242. x -= i
  243. if not x:
  244. break
  245. x = 1/x
  246. def continued_fraction_convergents(cf):
  247. """
  248. Return an iterator over the convergents of a continued fraction (cf).
  249. The parameter should be in either of the following to forms:
  250. - A list of partial quotients, possibly with the last element being a list
  251. of repeating partial quotients, such as might be returned by
  252. continued_fraction and continued_fraction_periodic.
  253. - An iterable returning successive partial quotients of the continued
  254. fraction, such as might be returned by continued_fraction_iterator.
  255. In computing the convergents, the continued fraction need not be strictly
  256. in canonical form (all integers, all but the first positive).
  257. Rational and negative elements may be present in the expansion.
  258. Examples
  259. ========
  260. >>> from sympy.core import pi
  261. >>> from sympy import S
  262. >>> from sympy.ntheory.continued_fraction import \
  263. continued_fraction_convergents, continued_fraction_iterator
  264. >>> list(continued_fraction_convergents([0, 2, 1, 2]))
  265. [0, 1/2, 1/3, 3/8]
  266. >>> list(continued_fraction_convergents([1, S('1/2'), -7, S('1/4')]))
  267. [1, 3, 19/5, 7]
  268. >>> it = continued_fraction_convergents(continued_fraction_iterator(pi))
  269. >>> for n in range(7):
  270. ... print(next(it))
  271. 3
  272. 22/7
  273. 333/106
  274. 355/113
  275. 103993/33102
  276. 104348/33215
  277. 208341/66317
  278. >>> it = continued_fraction_convergents([1, [1, 2]]) # sqrt(3)
  279. >>> for n in range(7):
  280. ... print(next(it))
  281. 1
  282. 2
  283. 5/3
  284. 7/4
  285. 19/11
  286. 26/15
  287. 71/41
  288. See Also
  289. ========
  290. continued_fraction_iterator, continued_fraction, continued_fraction_periodic
  291. """
  292. if isinstance(cf, list) and isinstance(cf[-1], list):
  293. cf = itertools.chain(cf[:-1], itertools.cycle(cf[-1]))
  294. p_2, q_2 = S.Zero, S.One
  295. p_1, q_1 = S.One, S.Zero
  296. for a in cf:
  297. p, q = a*p_1 + p_2, a*q_1 + q_2
  298. p_2, q_2 = p_1, q_1
  299. p_1, q_1 = p, q
  300. yield p/q