Coding Challenge for Vieta's Formulas

[Namir]

I am interested in coding general cases of the Vieta's formulas (check the Wikipedia article to see these formulas in a high-level language like VBA, Python, Matlab, and C++. The target is to have a function that sums up the squares of the LHS - RHS of each Vieta formula:

F(X,a) = (LHS1 - RHS1)^2 + (LHS2 - RHS2)^2 + ... + (LHSn - RHSn)^2

Where:

LHS1 = sum(X(i)) for i = 1 to n
RHS1 = - a(n-1)/a(n))

LHS2 = ((X(1)*X(2) + X(1)*X(3) + ... + X(1)*X(n)) +
((X(2)*X(3) + X(2)*X(4) + ... + X(2)*X(n)) +
... +
X(n-1)*X(n)
RHS2 = a(n-2)/a(n))
....
LHSn = product(X(i)) for i = 1 to n
RHSn = (-1)^n * a(0)/a(n)

The parameters are the array X() which represents the guesses for the roots, and array a() that represents the polynomial coefficients. The Wikipedia article designates a(n) as the coefficient for x^n and a(0) as the constant term.

Good luck to us all!!

Namir

[Valentin Albillo]

.
Hi, Namir|

5:30 am here so just a couple quick observations:

(07-12-2018 11:15 PM)Namir Wrote:  F(X,a) = (LHS1 - RHS1)^2 + (LHS2 - RHS2)^2 + ... + (LHSn - RHSn)^2

Where:

LHS1 = sum(X(i)) for i = 1 to n
RHS1 = - a(n-1)/a(n)
....
LHSn = product(X(i)) for i = 1 to n
RHSn = (-1)^2 * a(0)/a(n)

1) without loss of generality you can consider a(n) to be 1. If it's not, simply divide all other coefficients by it before doing anything else.

2) there's a typo in your RHSn, (-1)^2 is obviously wrong.

Regards.
V.
.

[Namir]

1) Yes a(n) can be 1.

2) Error corrected.

Thanks!

Namir

[Namir]

The Vieta formulas are related to the Elementary symmetric polynomial. Check the Wikipedia link.

The solution to my challenge lies in the ability to generate high-level programming language code for the Elementary Symmetric Polynomials.

[Valentin Albillo]

.
G'day, Namir:

Just a question:

(07-12-2018 11:15 PM)Namir Wrote:  The parameters are the array X() which represents the guesses for the roots, and array a() that represents the polynomial coefficients

As you say X() represents the guesses for the roots, which might be complex, then we must assume X() is a complex array, right ?

Also, you say the array a() represents the coefficients. Are they assumed to be complex in the general case or are you considering just polynomials with real coefficients ?

Finally, pretty please, correct that "Foummulas" (sic) in the thread title, it's a real eyesore. XD

Regards.
V.
.

[Namir]

(07-13-2018 02:46 PM)Valentin Albillo Wrote:  .
G'day, Namir:

Just a question:

(07-12-2018 11:15 PM)Namir Wrote:  The parameters are the array X() which represents the guesses for the roots, and array a() that represents the polynomial coefficients

As you say X() represents the guesses for the roots, which might be complex, then we must assume X() is a complex array, right ?

Also, you say the array a() represents the coefficients. Are they assumed to be complex in the general case or are you considering just polynomials with real coefficients ?

Finally, pretty please, correct that "Foummulas" (sic) in the thread title, it's a real eyesore. XD

Regards.
V.
.

Yes both arrays X() and a() can be complex and the mathematical operations work with complex numbers. For programming languages like Matlab, the operations and functions automatically handle real and complex numbers. For other languages (not sure about Python) one may have to use complex operators and functions just to be sure that the math expressions can handle real and complex numbers.

[Valentin Albillo]

.
Hi, Namir:

This is my entry for your challenge, which fully meets your requirements. I'm not sure if you consider HP-71B BASIC a "high-level" language but this 4-liner (199 bytes) delivers the goods in the most efficient, economical way and further, if you don't like 71B BASIC you can convert it to your preferred dialect (VBA, Excel, Maple, Python, HPPPL, even RPL) in 5 minutes flat as it's eminently understandable. Let's see:

>LIST

      1  DESTROY ALL @ OPTION BASE 0 @ INPUT "Deg: ";N @ COMPLEX A(N),B(N),X(N-1)
      2  MAT INPUT X,A @ MAT B=ZER @ B(N)=1 @ FOR I=1 TO N @ FOR J=N-I TO N-1
      3  B(J)=B(J+1)-B(J)*X(I-1) @ NEXT J @ B(N)=-B(N)*X(I-1) @ NEXT I @ MAT DISP B
      4  MAT A=(1/A(0))*A @ MAT B=B-A @ DISP "F(X,a) = ";FNORM(B)^2

Notes:

• Both X() and a() are declared as COMPLEX arrays so the program works for both real/complex roots and/or real/complex coefficients
• Complex values are represented as an ordered pair (a, b) = a + i*b. The pair (a, 0) is the real value a
• The coefficients of a() should ideally begin with 1 (i.e.: a(n)=1) for efficiency but the program works equally well with any non-zero a(n)
• If your preferred language doesn't have FNORM (Frobenius Norm), then use this line 4 instead:
  
        4  MAT A=(1/A(0))*A @ S=0 @ FOR I=1 TO N @ S=S+ABS(B(I)-A(I))^2 @ NEXT I @ DISP "F(X,a) =";S
  
  which gives the same result but takes more time and more RAM than FNORM.
• For complex arrays FNORM^2 is the sum of the squares of the absolute values of each (possibly complex) element.

1) Let's see some examples. First a 3rd degree polynomial with 3 real roots (2, 3 and 4):

      P(x) = 5*x^3 -45*x^2 +130*x -120

>RUN

Deg: 3
X(0)? 2.1,3.2,3.9    { approximations to the exact roots }
A(0)? 5,-45,130,-120 { coefficients, its roots are 2, 3 and 4 }

      (1,0)          { always 1 }
      (-9.2,-0)      { x1 + x2 + x3 = 2.1 + 3.2 + 3.9 = 9.2 }
      (27.39,0)      { x1*x2 + x1*x3 + x2*x3 = 2.1*3.2+2.1*3.9+3.2*3.9 = 27.39 }
      (-26.208,-0)   { x1 * x2 * x3 = 2.1*3.2*3.9 = 26.208 }

            F(X,a) = 6.847364 { (-9.2-(-9))^2+(27.39-26)^2+(-26.208-(-24))^2 = 6.847364 }

When the approximate roots finally converge to the exact roots, you get F(X,a) = 0 as it should:

>RUN

Deg: 3
X(0)? 2,3,4
A(0)? 5,-45,130,-120

      (1,0)
      (-9,-0)
      (26,0)
      (-24,-0)

             F(X,a) = 0

2) Let's see a more complicated example, this time a 6th degree polynomial with 2 real roots and 4 complex roots:

      P(x) = 5*x^6 -45*x^5 +225*x^4 -425*x^3 +170*x^2 +370*x -500

>RUN

Deg: 6
X(0)? (0.99,1.01),(0.99,-1.01),-1.02,1.98,(2.99,3.99),(2.99,-3.99)  { approximations }
A(0)? 5,-45,225,-425,170,370,-500                                   { its roots are (1,+-1), -1, 2, (3,+-4) }

      (1,0)                { always 1 }
      (-8.92,0)            { (0.99,1.01)+(0.99,-1.01)-1.02+1.98+(2.99,3.99)+(2.99,-3.99) = 8.92 }
      (44.3228,0)
      (-82.261144,0)
      (30.30225268,0)
      (75.8316409248,0)
      (-100.425361372,0)   { (0.99,1.01)*(0.99,-1.01)*(-1.02)*1.98*(2.99,3.99)*(2.99,-3.99) = -100.425361372 }

            F(X,a) = 25.1755080455

3) Finally, a 3rd degree polynomial with complex coefficients: (3 complex roots, none are conjugates):

      P(x) = x^3 +(-9 -12*i)*x^2 +(-21 +64*i)*x +(85 -20*i)

>RUN
Deg: 3

X(0)? (1.1,1.9),(2.9,4.1),(4.9,6.1)   { complex approximations }
A(0)? 1,(-9,-12),(-21,64),(85,-20)    { coefficients, its roots are (1,2), (3,4), (5,6) }

      (1,0)                { always 1 }
      (-8.9,-12.1)         { (1.1,1.9)+(2.9,4.1)+(4.9,6.1) = (8.9,12.1) }
      (-21.6,63.82)        { (1.1,1.9)*(2.9,4.1)+(1.1,1.9)*(4.9,6.1)+(2.9,4.1)*(4.9,6.1) = (-21.6,63.82) }
      (83.662,-21.038)     { (1.1,1.9)*(2.9,4.1)*(4.9,6.1) = (-83.662,21.038) }

            F(X,a) = 3.280088

Bst regards and have a nice weekend.
V.
.

[Namir]

Most fascinating reading and solution!!! The simplicity of your solution makes me speechless. I have not seen any code, on the Internet, for a general calculations for Vieta's Formulas!

Namir

[Namir]

I am still shaking my head at the simplicity and elegance of your solution. I have seen anything comparable on the Internet!!

My attempt to solve the challenge would take me thought a complicated and difficult route!! :-(

Namir

[Namir]

I found a C++ listing that calculates the polynomial coefficnts, given the roots of that polynomial, using the Vieta's Formulas. I checked the code using Visual Studio 2017 Community Edition. The C++ code had a small error in the range of the nested loop (for j = ...). I corrected it and then translated it into Excel VBA. Here is the version that uses Option Base 1 arrays:

Code:

Option Explicit
Option Base 1

Function Vieta(ByRef A() As Double, ByRef X() As Double, _
               ByRef B() As Double, ByVal N As Integer) As Double
  Dim I As Integer, J As Integer
  Dim Sum As Double, A1 As Double
  
  A1 = A(1)
  For I = 1 To N + 1
    A(I) = A(I) / A1
  Next I

  B(N + 1) = 1
  For I = 1 To N
    For J = N - I To N - 1
      B(J + 1) = B(J + 1) - X(I) * B(J + 2)
    Next J
  Next I
  
  Sum = 0
  For I = 1 To N + 1
    Sum = Sum + (A(I) - B(N + 2 - I)) ^ 2
  Next I
  
  Vieta = Sqr(Sum)
  
End Function

Sub go()
   Dim I As Integer, J As Integer, N As Integer
   Dim A() As Double, B() As Double, X() As Double
   
   Range("B7").Clear
   N = Range("A1").CurrentRegion.Rows.Count - 2
   ReDim A(N + 1), X(N), B(N + 1)

   For I = 1 To N
     A(I) = Cells(I + 1, 1)
     X(I) = Cells(I + 1, 2)
   Next I
   A(N + 1) = Cells(N + 2, 1)
  
   Cells(N + 3, 2) = Vieta(A, X, B, N)
   J = 1
   For I = N + 1 To 1 Step -1
     Cells(J + 1, 3) = B(I)
     J = J + 1
   Next I
End Sub

The above code uses the following spreadsheet columns (the first row has headings and so value appear as of row 2):

1) Column A has the polynomial coefficients with cell A2 being the coefficient of the higher term. There are N + 1 entries in this column.
2 Column B has N roots, starting with cell B2.
3) Column C has the output array B() which should match the values in the first columns if the roots match the polynomial coefficients. The first value appears in cell C2.
4) Cell B7 displays the value of the Vieta's Formula function.

The VBA code is similar to, but not an exact match, to Vallentin's HP-71B code. Note that the values in array B() are in reverse order of the values in array A().