Post Reply 
Coding Challenge for Vieta's Formulas
07-12-2018, 11:15 PM (This post was last modified: 07-13-2018 06:37 PM by Namir.)
Post: #1
Coding Challenge for Vieta's Formulas
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
Find all posts by this user
Quote this message in a reply
07-13-2018, 03:34 AM
Post: #2
RE: Coding Challenge for Vieta's Foummulas
.
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.
.

  
All My Articles & other Materials here:  Valentin Albillo's HP Collection
 
Visit this user's website Find all posts by this user
Quote this message in a reply
07-13-2018, 06:19 AM (This post was last modified: 07-13-2018 06:38 PM by Namir.)
Post: #3
RE: Coding Challenge for Vieta's Formulas
1) Yes a(n) can be 1.

2) Error corrected.

Thanks!

Namir
Find all posts by this user
Quote this message in a reply
07-13-2018, 11:50 AM (This post was last modified: 07-13-2018 06:36 PM by Namir.)
Post: #4
RE: Coding Challenge for Vieta's Foumulas
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.
Find all posts by this user
Quote this message in a reply
07-13-2018, 02:46 PM
Post: #5
RE: Coding Challenge for Vieta's Foummulas
.
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.
.

  
All My Articles & other Materials here:  Valentin Albillo's HP Collection
 
Visit this user's website Find all posts by this user
Quote this message in a reply
07-13-2018, 06:35 PM (This post was last modified: 07-13-2018 07:42 PM by Namir.)
Post: #6
RE: Coding Challenge for Vieta's Formulas
(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.
Find all posts by this user
Quote this message in a reply
07-14-2018, 12:50 AM
Post: #7
RE: Coding Challenge for Vieta's Formulas
.
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.
.

  
All My Articles & other Materials here:  Valentin Albillo's HP Collection
 
Visit this user's website Find all posts by this user
Quote this message in a reply
07-14-2018, 01:36 AM (This post was last modified: 07-14-2018 11:05 AM by Namir.)
Post: #8
RE: Coding Challenge for Vieta's Formulas
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
Find all posts by this user
Quote this message in a reply
07-14-2018, 03:13 PM
Post: #9
RE: Coding Challenge for Vieta's Formulas
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
Find all posts by this user
Quote this message in a reply
07-21-2018, 01:10 PM (This post was last modified: 07-21-2018 01:47 PM by Namir.)
Post: #10
RE: Coding Challenge for Vieta's Formulas
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().
Find all posts by this user
Quote this message in a reply
Post Reply 




User(s) browsing this thread: 1 Guest(s)