如何用SymPy扩展计算(m/n)π的正弦余弦精确值范围？

问题描述 #

我希望通过计算大量有理数倍π的sin和cos精确值来学习数学与Python，编写了如下代码，在特定条件下运行效果较好，但n的适用范围有限：可处理素数n={2,3,5,7,11}，以及形如2^k*3l的n（如n=144），但遇到n=25时会生成无法求解的五次方程，请问如何修改代码扩展n的适用范围？同时希望了解相关的数学公式、方法或理论作为学习参考。

原代码：

from sympy import *
import sympy
from IPython.display import display
from sympy.abc import k, m, n, x
import math
from sympy import init_printing
init_printing() 
print("You are solving for sin and cos of (m/n)pi")
m = int(input("Insert m: "))
n = int(input("Insert n: "))
gcd = sympy.gcd(m,n) #Getting the GCD and dividing by it will make the algorithm have less steps
m = m/gcd
n = n/gcd
prime = primefactors(n) #this outputs an array of all of the factors of n
cosp = sympify(-1) #this is necessary, we set the value of cos(pi) = -1, and format it for sympy
r = int(1)
for p in prime:#the nested loops will go through every factor. The for loop each of the factors,
    while n%p == 0:#the while loop for the repetition of the factors
        r = int(r*p)#the r variable serves as a comparison for the values of the cosine
        n = int(n/p)#dividing by the prime we're checking for will complete the loop
        if p == 2:
            cosp = sqrt((1+cosp)/2)#simple bisection of the cosine formula
        else:
            cosp = (Sum((-1)**k*factorial(p)/(factorial(p-2*k)*factorial(2*k))*x**(p-2*k)*(1-x**2)**k , (k, 0, int((p-1)/2)))).doit()-cosp #check out this explanation: https://brilliant.org/wiki/expansions-of-certain-trigonometric-functions/. This is a version of the cosine formula, where x symbolizes cosine pi, and we bring the cosine value to the other side to get "0 = F(x)"
            expandeq = expand(cosp)
            res = solve(expandeq) #we solve for the roots of the equation
            for root in res:
                if N(sympy.cos(sympy.pi/r), 15) == re(N(root, 15)): #to remove all extra solutions, we check the approximation of each root and verify it is equal to cos(pi/r)
                    if str(root)[0]!= "C":
                        cosp = root
sinp = expand(sqrt(1-(cosp**2))) 

cosm = expand(Sum((-1)**k*factorial(m)/(factorial(m-(2*k))*factorial(2*k))*cosp**(m-2*k)*sinp**(2*k) , (k, 0, int(m/2))).doit(), trig = False) #finally, we use the formula again to multiply the cosine of pi/n by m
sinm = expand(sympy.sqrt(1-(cosm**2)), trig=False)
print("Approximation of sine of pi*",m, "/",r,": " , round(math.sin(math.pi*m/r), 15))
display(sinm)
print("Approximation of of cosine pi*",m, "/",r,": " , round(math.cos(math.pi*m/r), 15))
display(cosm)

代码修改方案 #

原代码的核心问题是手动递推构造高次方程求解，当遇到五次及以上不可约方程时，solve无法直接给出简洁的根式解（即使某些方程理论上有根式解，手动构造的方程也会引入冗余复杂度）。更可靠的方式是直接利用SymPy内置的符号三角函数化简能力，它已经封装了分圆多项式、伽罗瓦理论相关的逻辑，能处理绝大多数有理数倍π的精确值计算：

简化后的代码 #

from sympy import *
init_printing()

print("计算 (m/n)π 的正弦和余弦精确值")
m = int(input("输入 m: "))
n = int(input("输入 n: "))

# 先对m/n约分，减少计算量
g = gcd(m, n)
m_reduced = m // g
n_reduced = n // g

# 直接构造有理数倍π的符号表达式，SymPy会自动化简为精确形式
theta = Rational(m_reduced, n_reduced) * pi
sin_exact = sin(theta).expand(trig=True)
cos_exact = cos(theta).expand(trig=True)

# 输出近似值与精确表达式
print(f"(m/n)π 的正弦近似值: {N(sin_exact, 15)}")
display(sin_exact)
print(f"(m/n)π 的余弦近似值: {N(cos_exact, 15)}")
display(cos_exact)

代码说明 #

1. 约分处理：通过gcd对m/n约分，避免重复计算相同角度的三角函数。
2. 符号计算：SymPy的sin和cos函数接收有理数倍π的符号输入时，会自动利用分圆多项式、三角恒等式化简为根式或其他精确形式（即使无法用简单根式表示，也会输出基于分圆多项式根的精确表达式）。
3. 扩展适用范围：该代码支持所有正整数n，包括n=25这类原代码无法处理的情况。

相关数学理论参考 #

1. 伽罗瓦理论 #

这是决定有理数倍π的三角函数能否用根式表示的核心理论：

• 当且仅当n的素因子仅为2和费马素数（形如2⁽²k)+1的素数，已知的有3,5,17,257,65537）时，cos(πm/n)可以用有限次加减乘除和开方表示为简单根式。
• 对于n包含其他素因子（如5^2=25），对应的分圆域伽罗瓦群虽然可解，但解的形式会是嵌套根式，手动构造方程难以处理，而SymPy已内置相关逻辑。

2. 分圆多项式 #

分圆多项式Φ_n(x)是首一整系数多项式，其根为n次本原单位根e^(2πik/n)（k与n互质）。利用恒等式：
$$\cos\left(\frac{2\pi k}{n}\right) = \frac{e^{2\pi ik/n} + e^{-2\pi ik/n}}{2}$$
可以通过分圆多项式构造cos的精确表达式，SymPy的三角函数化简正是基于这一原理。

3. 多倍角与三角恒等式 #

• 倍角公式：如cos(2θ)=2cos²θ-1，cos(3θ)=4cos³θ-3cosθ，可用于递推低角度的余弦值，但高次倍角公式会生成高次方程，手动求解效率低且易出错。

内容的提问来源于stack exchange，提问作者Nicolas Campailla