▼
Posts: 2,247
Threads: 200
Joined: Jun 2005
Hi All folks who are NOT attending HHC 2010,
Since you are not busy at the HHC, I thought you can handle the following worthwhile challenge.
The Viete's formulas allow you to calculate the coefficients of a polynomial given it's roots. For example for quadratic polynomials:
x1 + x2 =  b / a
x1 * x2 = c / a
Giving the polynomial a x^2 + b X + c. Choosing a = 1 simplifies things.
In the case of a cubic polynomial you have:
x1 + x2 + x3 =  b / a
x1* x2 + x1 * x3 + x2 * x3 = c /a
x1 * x2 * x3 = d / a
Giving the polynomial a x^3 + b x^2 + c x + d. Choosing a = 1 simplifies things.
Wikipedia shows the general equations for the Viete's formulas here.
All polynomials (order 2 ad higher) have the following in common:
1) The coefficient of the n1 term (n being the polynomial order) is calculated as:
sum(x(i)) =  b / a, for i = 1 to n
2) The constant term is calculated as:
product(x(i) = (1)^n * constant_term / a for i = 1 to n
Your challenge is to write a program in RPN, RPL, BASIC or any other popular language (for calculator or for PC) that can calculate the coefficients of any polynomial given the array of roots. Using coefficient a as 1 simplifies things a bit.
The solution requires that you handle a number of n equations to calculate the coefficients of a polynomial of order n.
Namir
AWOL from HHC 2010 and hiding in the south of France!! Merci Sarkozi!!!
Edited: 25 Sept 2010, 12:31 p.m.
▼
Posts: 735
Threads: 34
Joined: May 2007
The HP48GX has this function built in:
[1 2 3] PCOEF > [1 6 11 6]
You may consider it the inverse operation of PROOT.
Best regards
Thomas
▼
Posts: 2,247
Threads: 200
Joined: Jun 2005
You are not implementing a solution using Viete's formulas. I am aware of function PCOEF and PROOT.
Namir
Posts: 709
Threads: 104
Joined: Nov 2005
Are we shooting for program size, speed, or simply an implementation of the formula? =)
Han
Posts: 709
Threads: 104
Joined: Nov 2005
I didn't implement Viete's formula directly, but the algorithm does (in the end) do the exact multiplication and addition used in Viete's formula.
Given a list of roots: { x1 x2 x3 ... xn } it computes the list of coefficients.
<<
DUP SIZE > r n
<<
{ 1. } 0. n 1. 
FOR i
0. OVER + r n i  GET NEG *
SWAP 0. + SWAP ADD
NEXT
>>
>>
If anyone would like to see the math behind it, I can post that as well.
▼
Posts: 735
Threads: 34
Joined: May 2007
Same idea but using STREAM:
\<<
{{1}} SWAP +
\<< \> p x
\<<
p 0 +
0 p +
x * 
\>>
\>>
STREAM
\>>
Best regards
Thomas
▼
Posts: 2,247
Threads: 200
Joined: Jun 2005
Thomas,
Your solution works great too (input has to be a list of roots).
Cool!!
Namir
Posts: 2,247
Threads: 200
Joined: Jun 2005
Sure! Please do! I am very impressed with your solution.
You solution works great! Thanks!
Edited: 25 Sept 2010, 8:48 p.m.
▼
Posts: 709
Threads: 104
Joined: Nov 2005
The idea can be easily seen with the following example. Suppose our roots are 3, 5, and 7. Then the polynomial can be written in factored form is:
(x3)(x5)(x7)
First multiply the first two factors (x3)(x5).
(x3)*x + (x3)*(5) = (x^2  3x) + (5x + 15)
If we use vector (or list) notation, what we are essentially doing is taking { 1 3 } (which corresponds to the polynomial x3) and "shifting" it by a power of x (to get x^2  3x). So in list form, this is { 1 3 0 }. Then, we compute (x3)*(5) which gives us the list { 0 5 15 }. We need a 3element list though, since ADD requires same list sizes. Lastly, we add the two lists via elementwise addition. This is essentially the sum of two polynomials: one whose coefficient is { 1 3 0 } and the other is { 0 5 15 }. So the sum is { 1 8 15 }, or the polynomial x^2  8x + 15.
Now, just repeat this process with the final factor: (x7). It's really just distribution of x and the root x_i (negated).
This is best seen when using the DBUG feature of the HP48/49/50.
Han
Edited: 25 Sept 2010, 9:22 p.m.
▼
Posts: 709
Threads: 104
Joined: Nov 2005
This reminds of another challenge of computing f(n) for some given number n, where f(x) is a polynomial using an RPN calculator that only has 2 levels of stack.
The idea is to rewrite the polynomial in a "recursive" manner. For example, if we have 3x^25x+6, we can view this as:
(((x * 3)  5 ) * x) + 6
So to compute f(4), for example, we would do:
4
3
*
5

4
*
6
+
▼
Posts: 2,247
Threads: 200
Joined: Jun 2005
I see .. so you are doing the polynomial multiplication approach. This is the easy route. I was able to write a version of PCOEF in Matlab (similar to Matlab's ROOTS function).
The Viete's formulas are a bit tough to code. Hence the challenge!!
:)
Namir
Posts: 275
Threads: 41
Joined: Mar 2010
Hi,
Here is a brutal solution of Viète formula (non optimized, lot of room for optimizing I think) for HP86/7 with nice display of sums.
Olivier
10 DISP "Viete's formulas" @ OPTION BASE 0
20 DIM x(99),a(99),ik(99)
30 RESTORE 1000 @ READ n@ FOR k=1 TO n @ READ x(k)@ NEXT k @ a(n)=1
40 FOR k=1 TO n @ s=0 @ FOR i=1 TO k @ ik(i)=i @ NEXT i @ l=k @ DISP "S";nk;" =";
50 p=1 @ FOR i=1 TO k @ p=p*x(ik(i)) @ DISP ik(i);@ NEXT i @ s=s+p
60 ik(l)=ik(l)+1 @ IF ik(l)<= nk+l THEN DISP "+";@ GOTO 50
70 l=l1 @ IF l=0 THEN 200
80 IF ik(l)=nk+l THEN 70
90 j=ik(l)+1 @ FOR i=l TO k @ ik(i)=j @ j=j+1 @ NEXT i @ l=k @ DISP "+";@ GOTO 50
200 a(nk)=s*((1)^k) @ DISP "" @ DISP "a(";nk;")=";a(nk) @ NEXT k
999 END
1000 DATA 5,1,2,3,4,5
Edited: 26 Sept 2010, 11:47 a.m.
▼
Posts: 2,247
Threads: 200
Joined: Jun 2005
Olivier,
Quoting the famous French mathematician Bernoulli, let me say "I recognize the lion by its tale!"
The quote comes from when Bernoulli challenged the mathematicians of Europe to solve a problem. At that time, Newton was sick, so Bernoulli did NOT send him a copy of the problem. After waiting several months, Bernoulli got frustrated and upon hearing that Newton's health was better, he mailed Newton a copy of the problem. Upon receipt Newton looked at the problem and mailed THE SOLUTION THE NEXT DAY!!! When Bernoulli got an envelope from England that quickly (relative to speed of mail at the time), he knew it was from Newton and said, before opening the envelope, "I recognize the lion by its tale!"
I am not surprised that you solved it, since it takes a first class programmer like yourself to solve a tough problem like implementing the Viete's formulas!
Thanks!!
▼
Posts: 275
Threads: 41
Joined: Mar 2010
Thanks a lot for your nice comments Namir.
Just a little remark ... could it be "paw" in place of "tale" in the quote of Bernouilli ? (see cite)
Paw (as "patte" in french) seems better I think.
Anyway thanks again for your kind comments.
Olivier
Edited: 27 Sept 2010, 2:07 a.m.
Posts: 2,247
Threads: 200
Joined: Jun 2005
I typed in your solution in the EMU71 (I have no physical 71B with me while I am on vacation). I changed the name of the array ik() to i1() (a letter and a digit) to make the code work. Other than that the code is a beauty!
I am impressed on how you were able to use a twolevel nested FOR loops to manage calculating the various Viete's summations (and consequently the related polynomial coefficients). I kept thinking that the solution needed more nested loops.
Cool!
Namir
Edited: 26 Sept 2010, 5:56 p.m.
Posts: 275
Threads: 41
Joined: Mar 2010
Optimizing a bit (remove an unnecessary fornext)
10 DISP "Viete's formulas" @ OPTION BASE 0@ DIM x(99),a(99),ik(99)
20 RESTORE @ READ n@ FOR k=1 TO n @ READ x(k)@ NEXT k @ a(n)=1
30 FOR k=1 TO n @ s=0 @ l=1 @ ik(1)=0 @ DISP "S";nk;"=";
40 j=ik(l)+1 @ FOR i=l TO k @ ik(i)=j @ j=j+1 @ NEXT i @ l=k
50 p=1 @ FOR i=1 TO k @ p=p*x(ik(i)) @ DISP ik(i);@ NEXT i @ s=s+p
60 ik(l)=ik(l)+1 @ IF ik(l)<= nk+l THEN DISP "+";@ GOTO 50
70 l=l1 @ IF l=0 THEN 90
80 IF ik(l)=nk+l THEN 70 ELSE DISP "+";@ GOTO 40
90 a(nk)=s*((1)^k) @ DISP "" @ DISP "a(";nk;")=";a(nk) @ NEXT k
100 END
110 DATA 5,1,2,3,4,5
Olivier
▼
Posts: 2,247
Threads: 200
Joined: Jun 2005
Thanks!!! I was able to translate the 71B code into Excel VBA. Had to watch for those zigzagging GOTOs!!
Here is the VBA code:
Option Explicit
Option Base 0
' Version 1.1
Sub CalcCoeff()
Dim N As Integer, S As Double, P As Double
Dim I As Integer, J As Integer, K As Integer, M As Integer
Dim Root(99) As Double, Coeff(99) As Double, IK(99) As Integer
Dim bJump As Boolean, Row As Integer, Col As Integer
Dim aStr As String
N = 2
Do While Cells(N, 1) <> ""
Root(N  1) = Cells(N, 1)
N = N + 1
Loop
N = N  2 ' number of roots
Coeff(N) = 1
' display Coeff(N)
Cells(2, 2) = Coeff(N)
Range("B3:Z1000").Value = ""
Row = 3
' main loop
For K = 1 To N
Col = 3
S = 0
M = 1
IK(1) = 0
' Line 40
Do
J = IK(M) + 1
For I = M To K
IK(I) = J
J = J + 1
Next I
M = K
' Line 50
Do
P = 1
For I = 1 To K
P = P * Root(IK(I))
aStr = Cells(Row, Col)
If aStr = "" Then
Cells(Row, Col) = "r" & CStr(IK(I))
Else
Cells(Row, Col) = aStr & "*r" & CStr(IK(I))
End If
Next I
S = S + P
' L60
IK(M) = IK(M) + 1
' move column for Viete's terms to the right
If IK(M) <= N  K + M Then Col = Col + 1
Loop While IK(M) <= N  K + M
bJump = False
Do
M = M  1
If M = 0 Then
bJump = True
Exit Do
End If
' move column for Viete's terms to the right
If IK(M) <> (N  K + M) Then Col = Col + 1
Loop While IK(M) = N  K + M
' jump out of outer loop?
If bJump Then Exit Do
Loop
' Line 200
Coeff(N  K) = S * ((1) ^ K)
' show the coefficient value
Cells(Row, 2) = Coeff(N  K)
Row = Row + 1
Next K
End Sub
I replaced variables l, x() and a() with M, Root() and Coeff(), respectively. I also added the variable Row to display the polynomial coefficients in a descending power order. The program expects the roots to appear in column 1 starting with Row 2. The output polynomial coefficients appear in column 2, starting with row 2. Row 1 for columns 1 and 2 are reserved for labeling the columns. The output places the coefficients of the polynomial in a descending power order. Row 2 has 1 (the term for X^n), row 3 has the coefficient for X^(n1), and so on. Row n+2 has the constant term.
The above code also shows the different terms in Viete's formula. These symbolic terms appear in columns C, D, and up. Each cell contains a set of root products. The variable Col controls which column each root product group appears.
Namir
Edited: 27 Sept 2010, 11:12 p.m.
