The Museum of HP Calculators

HP Forum Archive 19

 Programming ChallengeMessage #1 Posted by Namir on 25 Sept 2010, 12:18 p.m. 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 n-1 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.

 Re: Programming ChallengeMessage #2 Posted by Thomas Klemm on 25 Sept 2010, 5:03 p.m.,in response to message #1 by Namir The HP-48GX 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

 Re: Programming ChallengeMessage #3 Posted by Namir on 25 Sept 2010, 5:38 p.m.,in response to message #2 by Thomas Klemm You are not implementing a solution using Viete's formulas. I am aware of function PCOEF and PROOT. Namir

 Re: Programming ChallengeMessage #4 Posted by Han on 25 Sept 2010, 5:59 p.m.,in response to message #1 by Namir Are we shooting for program size, speed, or simply an implementation of the formula? =) Han

 Re: Programming ChallengeMessage #5 Posted by Han on 25 Sept 2010, 6:24 p.m.,in response to message #1 by Namir 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.

 Re: Programming ChallengeMessage #6 Posted by Thomas Klemm on 25 Sept 2010, 7:32 p.m.,in response to message #5 by Han Same idea but using STREAM: ```\<< {{1}} SWAP + \<< \-> p x \<< p 0 + 0 p + x * - \>> \>> STREAM \>> ``` Best regards Thomas

 Re: Programming ChallengeMessage #7 Posted by Namir on 25 Sept 2010, 8:42 p.m.,in response to message #6 by Thomas Klemm Thomas, Your solution works great too (input has to be a list of roots). Cool!! Namir

 Re: Programming ChallengeMessage #8 Posted by Namir on 25 Sept 2010, 7:55 p.m.,in response to message #5 by Han Sure! Please do! I am very impressed with your solution. You solution works great! Thanks! Edited: 25 Sept 2010, 8:48 p.m.

 Re: Programming ChallengeMessage #9 Posted by Han on 25 Sept 2010, 9:16 p.m.,in response to message #8 by Namir 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: (x-3)(x-5)(x-7) First multiply the first two factors (x-3)(x-5). (x-3)*x + (x-3)*(-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 x-3) 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 (x-3)*(-5) which gives us the list { 0 -5 15 }. We need a 3-element list though, since ADD requires same list sizes. Lastly, we add the two lists via element-wise 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: (x-7). 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.

 Re: Programming ChallengeMessage #10 Posted by Han on 25 Sept 2010, 9:33 p.m.,in response to message #9 by Han 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^2-5x+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 +

 Re: Programming ChallengeMessage #11 Posted by Namir on 26 Sept 2010, 6:14 a.m.,in response to message #10 by Han 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

 Re: Programming ChallengeMessage #12 Posted by Olivier De Smet on 26 Sept 2010, 11:20 a.m.,in response to message #1 by Namir 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";n-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)<= n-k+l THEN DISP "+";@ GOTO 50 70 l=l-1 @ IF l=0 THEN 200 80 IF ik(l)=n-k+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(n-k)=s*((-1)^k) @ DISP "" @ DISP "a(";n-k;")=";a(n-k) @ NEXT k 999 END 1000 DATA 5,1,2,3,4,5 ``` Edited: 26 Sept 2010, 11:47 a.m.

 Re: Programming ChallengeMessage #13 Posted by Namir on 26 Sept 2010, 1:26 p.m.,in response to message #12 by Olivier De Smet 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!!

 Re: Programming ChallengeMessage #14 Posted by Olivier De Smet on 26 Sept 2010, 2:13 p.m.,in response to message #13 by Namir 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.

 Re: Programming ChallengeMessage #15 Posted by Namir on 26 Sept 2010, 5:50 p.m.,in response to message #12 by Olivier De Smet 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 two-level 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.

 Re: Programming ChallengeMessage #16 Posted by Olivier De Smet on 27 Sept 2010, 8:03 a.m.,in response to message #12 by Olivier De Smet Optimizing a bit (remove an unnecessary for-next) ```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";n-k;"="; 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)<= n-k+l THEN DISP "+";@ GOTO 50 70 l=l-1 @ IF l=0 THEN 90 80 IF ik(l)=n-k+l THEN 70 ELSE DISP "+";@ GOTO 40 90 a(n-k)=s*((-1)^k) @ DISP "" @ DISP "a(";n-k;")=";a(n-k) @ NEXT k 100 END 110 DATA 5,1,2,3,4,5 ``` Olivier

 Re: Programming ChallengeMessage #17 Posted by Namir on 27 Sept 2010, 1:13 p.m.,in response to message #16 by Olivier De Smet 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^(n-1), 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.

Go back to the main exhibit hall