_trigonometric_special.py 7.1 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261
  1. r"""A module for special angle formulas for trigonometric functions
  2. TODO
  3. ====
  4. This module should be developed in the future to contain direct square root
  5. representation of
  6. .. math
  7. F(\frac{n}{m} \pi)
  8. for every
  9. - $m \in \{ 3, 5, 17, 257, 65537 \}$
  10. - $n \in \mathbb{N}$, $0 \le n < m$
  11. - $F \in \{\sin, \cos, \tan, \csc, \sec, \cot\}$
  12. Without multi-step rewrites
  13. (e.g. $\tan \to \cos/\sin \to \cos/\sqrt \to \ sqrt$)
  14. or using chebyshev identities
  15. (e.g. $\cos \to \cos + \cos^2 + \cdots \to \sqrt{} + \sqrt{}^2 + \cdots $),
  16. which are trivial to implement in sympy,
  17. and had used to give overly complicated expressions.
  18. The reference can be found below, if anyone may need help implementing them.
  19. References
  20. ==========
  21. .. [*] Gottlieb, Christian. (1999). The Simple and straightforward construction
  22. of the regular 257-gon. The Mathematical Intelligencer. 21. 31-37.
  23. 10.1007/BF03024829.
  24. .. [*] https://resources.wolframcloud.com/FunctionRepository/resources/Cos2PiOverFermatPrime
  25. """
  26. from __future__ import annotations
  27. from typing import Callable
  28. from functools import reduce
  29. from sympy.core.expr import Expr
  30. from sympy.core.singleton import S
  31. from sympy.core.intfunc import igcdex
  32. from sympy.core.numbers import Integer
  33. from sympy.functions.elementary.miscellaneous import sqrt
  34. from sympy.core.cache import cacheit
  35. def migcdex(*x: int) -> tuple[tuple[int, ...], int]:
  36. r"""Compute extended gcd for multiple integers.
  37. Explanation
  38. ===========
  39. Given the integers $x_1, \cdots, x_n$ and
  40. an extended gcd for multiple arguments are defined as a solution
  41. $(y_1, \cdots, y_n), g$ for the diophantine equation
  42. $x_1 y_1 + \cdots + x_n y_n = g$ such that
  43. $g = \gcd(x_1, \cdots, x_n)$.
  44. Examples
  45. ========
  46. >>> from sympy.functions.elementary._trigonometric_special import migcdex
  47. >>> migcdex()
  48. ((), 0)
  49. >>> migcdex(4)
  50. ((1,), 4)
  51. >>> migcdex(4, 6)
  52. ((-1, 1), 2)
  53. >>> migcdex(6, 10, 15)
  54. ((1, 1, -1), 1)
  55. """
  56. if not x:
  57. return (), 0
  58. if len(x) == 1:
  59. return (1,), x[0]
  60. if len(x) == 2:
  61. u, v, h = igcdex(x[0], x[1])
  62. return (u, v), h
  63. y, g = migcdex(*x[1:])
  64. u, v, h = igcdex(x[0], g)
  65. return (u, *(v * i for i in y)), h
  66. def ipartfrac(*denoms: int) -> tuple[int, ...]:
  67. r"""Compute the partial fraction decomposition.
  68. Explanation
  69. ===========
  70. Given a rational number $\frac{1}{q_1 \cdots q_n}$ where all
  71. $q_1, \cdots, q_n$ are pairwise coprime,
  72. A partial fraction decomposition is defined as
  73. .. math::
  74. \frac{1}{q_1 \cdots q_n} = \frac{p_1}{q_1} + \cdots + \frac{p_n}{q_n}
  75. And it can be derived from solving the following diophantine equation for
  76. the $p_1, \cdots, p_n$
  77. .. math::
  78. 1 = p_1 \prod_{i \ne 1}q_i + \cdots + p_n \prod_{i \ne n}q_i
  79. Where $q_1, \cdots, q_n$ being pairwise coprime implies
  80. $\gcd(\prod_{i \ne 1}q_i, \cdots, \prod_{i \ne n}q_i) = 1$,
  81. which guarantees the existence of the solution.
  82. It is sufficient to compute partial fraction decomposition only
  83. for numerator $1$ because partial fraction decomposition for any
  84. $\frac{n}{q_1 \cdots q_n}$ can be easily computed by multiplying
  85. the result by $n$ afterwards.
  86. Parameters
  87. ==========
  88. denoms : int
  89. The pairwise coprime integer denominators $q_i$ which defines the
  90. rational number $\frac{1}{q_1 \cdots q_n}$
  91. Returns
  92. =======
  93. tuple[int, ...]
  94. The list of numerators which semantically corresponds to $p_i$ of the
  95. partial fraction decomposition
  96. $\frac{1}{q_1 \cdots q_n} = \frac{p_1}{q_1} + \cdots + \frac{p_n}{q_n}$
  97. Examples
  98. ========
  99. >>> from sympy import Rational, Mul
  100. >>> from sympy.functions.elementary._trigonometric_special import ipartfrac
  101. >>> denoms = 2, 3, 5
  102. >>> numers = ipartfrac(2, 3, 5)
  103. >>> numers
  104. (1, 7, -14)
  105. >>> Rational(1, Mul(*denoms))
  106. 1/30
  107. >>> out = 0
  108. >>> for n, d in zip(numers, denoms):
  109. ... out += Rational(n, d)
  110. >>> out
  111. 1/30
  112. """
  113. if not denoms:
  114. return ()
  115. def mul(x: int, y: int) -> int:
  116. return x * y
  117. denom = reduce(mul, denoms)
  118. a = [denom // x for x in denoms]
  119. h, _ = migcdex(*a)
  120. return h
  121. def fermat_coords(n: int) -> list[int] | None:
  122. """If n can be factored in terms of Fermat primes with
  123. multiplicity of each being 1, return those primes, else
  124. None
  125. """
  126. primes = []
  127. for p in [3, 5, 17, 257, 65537]:
  128. quotient, remainder = divmod(n, p)
  129. if remainder == 0:
  130. n = quotient
  131. primes.append(p)
  132. if n == 1:
  133. return primes
  134. return None
  135. @cacheit
  136. def cos_3() -> Expr:
  137. r"""Computes $\cos \frac{\pi}{3}$ in square roots"""
  138. return S.Half
  139. @cacheit
  140. def cos_5() -> Expr:
  141. r"""Computes $\cos \frac{\pi}{5}$ in square roots"""
  142. return (sqrt(5) + 1) / 4
  143. @cacheit
  144. def cos_17() -> Expr:
  145. r"""Computes $\cos \frac{\pi}{17}$ in square roots"""
  146. return sqrt(
  147. (15 + sqrt(17)) / 32 + sqrt(2) * (sqrt(17 - sqrt(17)) +
  148. sqrt(sqrt(2) * (-8 * sqrt(17 + sqrt(17)) - (1 - sqrt(17))
  149. * sqrt(17 - sqrt(17))) + 6 * sqrt(17) + 34)) / 32)
  150. @cacheit
  151. def cos_257() -> Expr:
  152. r"""Computes $\cos \frac{\pi}{257}$ in square roots
  153. References
  154. ==========
  155. .. [*] https://math.stackexchange.com/questions/516142/how-does-cos2-pi-257-look-like-in-real-radicals
  156. .. [*] https://r-knott.surrey.ac.uk/Fibonacci/simpleTrig.html
  157. """
  158. def f1(a: Expr, b: Expr) -> tuple[Expr, Expr]:
  159. return (a + sqrt(a**2 + b)) / 2, (a - sqrt(a**2 + b)) / 2
  160. def f2(a: Expr, b: Expr) -> Expr:
  161. return (a - sqrt(a**2 + b))/2
  162. t1, t2 = f1(S.NegativeOne, Integer(256))
  163. z1, z3 = f1(t1, Integer(64))
  164. z2, z4 = f1(t2, Integer(64))
  165. y1, y5 = f1(z1, 4*(5 + t1 + 2*z1))
  166. y6, y2 = f1(z2, 4*(5 + t2 + 2*z2))
  167. y3, y7 = f1(z3, 4*(5 + t1 + 2*z3))
  168. y8, y4 = f1(z4, 4*(5 + t2 + 2*z4))
  169. x1, x9 = f1(y1, -4*(t1 + y1 + y3 + 2*y6))
  170. x2, x10 = f1(y2, -4*(t2 + y2 + y4 + 2*y7))
  171. x3, x11 = f1(y3, -4*(t1 + y3 + y5 + 2*y8))
  172. x4, x12 = f1(y4, -4*(t2 + y4 + y6 + 2*y1))
  173. x5, x13 = f1(y5, -4*(t1 + y5 + y7 + 2*y2))
  174. x6, x14 = f1(y6, -4*(t2 + y6 + y8 + 2*y3))
  175. x15, x7 = f1(y7, -4*(t1 + y7 + y1 + 2*y4))
  176. x8, x16 = f1(y8, -4*(t2 + y8 + y2 + 2*y5))
  177. v1 = f2(x1, -4*(x1 + x2 + x3 + x6))
  178. v2 = f2(x2, -4*(x2 + x3 + x4 + x7))
  179. v3 = f2(x8, -4*(x8 + x9 + x10 + x13))
  180. v4 = f2(x9, -4*(x9 + x10 + x11 + x14))
  181. v5 = f2(x10, -4*(x10 + x11 + x12 + x15))
  182. v6 = f2(x16, -4*(x16 + x1 + x2 + x5))
  183. u1 = -f2(-v1, -4*(v2 + v3))
  184. u2 = -f2(-v4, -4*(v5 + v6))
  185. w1 = -2*f2(-u1, -4*u2)
  186. return sqrt(sqrt(2)*sqrt(w1 + 4)/8 + S.Half)
  187. def cos_table() -> dict[int, Callable[[], Expr]]:
  188. r"""Lazily evaluated table for $\cos \frac{\pi}{n}$ in square roots for
  189. $n \in \{3, 5, 17, 257, 65537\}$.
  190. Notes
  191. =====
  192. 65537 is the only other known Fermat prime and it is nearly impossible to
  193. build in the current SymPy due to performance issues.
  194. References
  195. ==========
  196. https://r-knott.surrey.ac.uk/Fibonacci/simpleTrig.html
  197. """
  198. return {
  199. 3: cos_3,
  200. 5: cos_5,
  201. 17: cos_17,
  202. 257: cos_257
  203. }