如何求解多项式的所有实数根呢？用什么样的算法呢？

郭先抢

任意给一个多项式，如何求解这个多项式的所有实数根呢？其中系数都是实数！

郭先抢

http://blog.163.com/shikang999@1 ... 896201131702358523/
这儿有个答案，但是我不知道原理，原理没有说

郭先抢

函数格式:Math_Matrix_PolyRoots1(A,A_n,Ret,Erro)

A:1*A_n的多项式系数矩阵

A_n:A的列数或大小

Ret:返回的实数根

Erro:误差控制

附加说明:

Math_Matrix_PolySturm:求多项式在某一点的变号数

Math_Matrix_PolyRoots1_Call:获取求解区间,气源代码在下面

Math_Matrix_PolyDiv:多项式除法

特别说明:执行本函数前,最好先处理重根。然后执行本函数后,根据函数返回值的不同以确定是否继续求解复数根。如果返回0表示多项式的根都为实数根,这是Ret里面就是全部的根,否则(不管返回1还是-1)都得再求复数根,这时A里面的多项式就只存在复数解,则可以调用【求多项式复数根贝尔斯托(Bairstow)算法Math_Matrix_PolyRoots2】函数求解复数根。

源代码:

Public Function Math_Matrix_PolyRoots1(ByRef A(,) As Double, ByVal A_n As Integer, ByRef Ret() As Double, ByVal Erro As Integer) As Int16

        '本函数求解多项式,求解其实根到Ret里面

        '如果A里面存在虚根,返回False并且返回求虚根多项式A;A里面全是实根,返回True

        '使用本函数前,请确定多项式没有重根

        '返回1表示求得有根,但还有虚根没有求解

        '返回0表示求得全部根,且无虚根

        '返回-1,没有求得根,即多项式只存在虚根

        Dim min As Double = Math.Abs(A(0, 1))

        Dim Max As Double = Math.Abs(A(0, 0))

        Dim i As Integer

        Dim j As Integer

        A_n -= 1

        For i = 2 To A_n

            If Math.Abs(A(0, i)) > min Then

                min = Math.Abs(A(0, i))

            End If

        Next

        min /= Max

        Max = 1 + min

        min = -Max

        '在区间[min,max]求解

        i = Math_Matrix_PolySturm(A, min)

        j = Math_Matrix_PolySturm(A, Max)

        i = Math.Abs(i - j)

        If i > 0 Then

            Dim N As Integer = i

            Dim temp(2) As Double

            Dim mtemp As Double

            Dim erro1 As Double = Math.Pow(0.1, Erro)

            N -= 1

            ReDim Ret(N)

            Dim Qujian(N, 1) As Double

            Dim Loopnumber As Int16

            If N > 0 Then

                Math_Matrix_PolyRoots1_Call(A, min, Max, N + 1, Qujian)

            Else

                Qujian(0, 0) = min

                Qujian(0, 1) = Max

            End If

            For i = 0 To N

                temp(0) = A(0, A_n)

                For j = 0 To A_n - 1

                    temp(0) += A(0, j) * Math.Pow(Qujian(i, 0), A_n - j)

                Next

                If Math.Abs(temp(0)) <= erro1 Then

                    Ret(i) = Qujian(i, 0)

                Else

                    temp(2) = A(0, A_n)

                    For j = 0 To A_n - 1

                        temp(2) += A(0, j) * Math.Pow(Qujian(i, 1), A_n - j)

                    Next

                    If Math.Abs(temp(2)) <= erro1 Then

                        Ret(i) = Qujian(i, 1)

                    Else

                        temp(1) = 1

                        Loopnumber = 1000

                        While temp(1) > erro1

                            mtemp = (Qujian(i, 0) + Qujian(i, 1)) / 2

                            temp(1) = A(0, A_n)

                            For j = 0 To A_n - 1

                                temp(1) += A(0, j) * Math.Pow(mtemp, A_n - j)

                            Next

                            If Math.Sign(temp(1)) = Math.Sign(temp(0)) Then

                                Qujian(i, 0) = mtemp

                            Else

                                Qujian(i, 1) = mtemp

                            End If

                            temp(1) = Math.Abs(temp(1))

                            Loopnumber -= 1

                            If Loopnumber = 0 Then

                                Exit While

                            End If

                        End While

                        Ret(i) = mtemp

                    End If

                End If

            Next

            If N = A.Length - 2 Then

                Return 0

            Else

                Dim mod1(0, 0) As Double

                Dim Atemp(0, 1) As Double

                For i = 0 To N

                    Atemp(0, 0) = 1

                    Atemp(0, 1) = -Ret(i)

                    Math_Matrix_PolyDiv(A, Atemp, Nothing, mod1, 11)

                    j = mod1.Length - 1

                    ReDim A(0, j)

                    While j >= 0

                        A(0, j) = mod1(0, j)

                        j -= 1

                    End While

                Next

                Return 1

            End If

        Else

            Return -1

        End If

    End Function

    Public Sub Math_Matrix_PolyRoots1_Call(ByVal A(,) As Double, ByVal min As Double, ByVal max As Double, ByVal N As Integer, ByRef Ret(,) As Double)

        '返回求解区域

        Dim i As Integer = Math_Matrix_PolySturm(A, min)

        Dim j As Integer = Math_Matrix_PolySturm(A, max)

        Dim temp As Double = (max + min) / 2

        Dim k As Integer = Math_Matrix_PolySturm(A, temp)

        If N > 1 Then

            i -= k

            i = Math.Abs(i)

            j -= k

            j = Math.Abs(j)

            If i = 0 Then

                Math_Matrix_PolyRoots1_Call(A, temp, max, N, Ret)

            ElseIf j = 0 Then

                Math_Matrix_PolyRoots1_Call(A, min, temp, N, Ret)

            Else

                If i = 1 Then

                    Ret(0, 0) = min

                    Ret(0, 1) = temp

                    N -= 1

                    If N = 1 Then

                        Ret(1, 0) = temp

                        Ret(1, 1) = max

                    Else

                        Dim ret1(N - 1, 1) As Double

                        Math_Matrix_PolyRoots1_Call(A, temp, max, N, ret1)

                        For i = 0 To N - 1

                            Ret(i + 1, 0) = ret1(i, 0)

                            Ret(i + 1, 1) = ret1(i, 1)

                        Next

                    End If

                ElseIf j = 1 Then

                    Ret(0, 0) = temp

                    Ret(0, 1) = max

                    N -= 1

                    Dim ret1(N - 1, 1) As Double

                    Math_Matrix_PolyRoots1_Call(A, min, temp, N, ret1)

                    For i = 0 To N - 1

                        Ret(i + 1, 0) = ret1(i, 0)

                        Ret(i + 1, 1) = ret1(i, 1)

                    Next

                Else

                    Dim ret1(i - 1, 1) As Double

                    N = i - 1

                    Math_Matrix_PolyRoots1_Call(A, min, temp, i, ret1)

                    For i = 0 To N

                        Ret(i, 0) = ret1(i, 0)

                        Ret(i, 1) = ret1(i, 1)

                    Next

                    ReDim ret1(j - 1, 1)

                    i = j - 1

                    Math_Matrix_PolyRoots1_Call(A, temp, max, j, ret1)

                    N += 1

                    For j = 0 To i

                        Ret(j + N, 0) = ret1(j, 0)

                        Ret(j + N, 1) = ret1(j, 1)

                    Next

                End If

            End If

        End If

    End Sub

郭先抢

Polynomial Root-finder (Real Coefficients)
http://www.akiti.ca/PolyRootRe.html

郭先抢

http://www.akiti.ca/rpoly_ak1_cpp.html

郭先抢

http://www.netlib.org/toms/493
读那么长的代码，也是一种折磨呀！！！！！！！！！！！！！！！！！！！！1

郭先抢

General Complex Polynomial Root Solver and Its Further Optimization for Binary Microlenses

http://arxiv.org/abs/1203.1034