Programming Challenge



#17

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.


#18

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


#19

You are not implementing a solution using Viete's formulas. I am aware of function PCOEF and PROOT.

Namir

#20

Are we shooting for program size, speed, or simply an implementation of the formula? =)

Han

#21

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.


#22

Same idea but using STREAM:

\<<
{{1}} SWAP +
\<< \-> p x
\<<
p 0 +
0 p +
x * -
\>>
\>>
STREAM
\>>

Best regards

Thomas


#23

Thomas,

Your solution works great too (input has to be a list of roots).

Cool!!

Namir

#24

Sure! Please do! I am very impressed with your solution.

You solution works great! Thanks!


Edited: 25 Sept 2010, 8:48 p.m.


#25

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.


#26

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

+


#27

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

#28

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.


#29

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!!


#30

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.

#31

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.

#32

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


#33

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.


Possibly Related Threads...
Thread Author Replies Views Last Post
  HHC 2013 RPN Programming Challenge - Final Thought Jeff O. 7 398 10-17-2013, 11:02 AM
Last Post: Jeff O.
  HHC 2012 RPN Programming Challenge Conundrum Jeff O. 15 564 10-08-2012, 03:34 PM
Last Post: Gerson W. Barbosa
  Mini-challenge: HHC2012 RPL programming contest with larger input David Hayden 14 641 10-05-2012, 10:36 PM
Last Post: David Hayden
  Weekend programming challenge: Euler's Totient function Allen 36 1,191 06-03-2012, 10:39 PM
Last Post: Paul Dale
  A Sunday Programming Challenge Allen 61 1,220 03-17-2012, 06:48 PM
Last Post: Allen
  HHC2012 programming challenge?? Pal G. 28 930 09-25-2011, 05:02 PM
Last Post: Paul Dale
  HHC MMC Programming challenge inC code David Hayden 8 378 12-31-2010, 09:40 AM
Last Post: David Hayden
  RE: 35s sorting routine challenge - Gene's Challenge Miguel Toro 4 272 08-01-2007, 08:36 AM
Last Post: Miguel Toro
  HP-16C programming challenge - Hamming distance Eric Smith 9 422 01-27-2006, 12:53 AM
Last Post: Eric Smith
  a simple programming challenge Don Shepherd 14 527 05-31-2005, 01:09 PM
Last Post: Valentin Albillo

Forum Jump: