Project Lovecraft Problem 5

[danbaron]

Problem 5:

Find the 5 positive integer roots (the 5 values of X which cause the equation to be true) for the following polynomial.

X^5 - 973447643102 * X^4 + 1132119607924548 * X^3 - 340211189835457006 * X^2 + 6379933987025808250 * X - 6040853943350632690 = 0

DTB.

[danbaron]

Here is one method of solution.

Dan



'~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
'file = Lovecraft5.tbasicc
'~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
Uses "console"
Uses "file"
'~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

%maxdegree = 5
%eps = 10^-12 As Ext

Global root(%maxdegree) As Quad
Global polycoeff(%maxdegree + 1) As Quad
Global derivcoeff(%maxdegree + 1) As Quad
Global reducedpolycoeff(%maxdegree + 1) As Quad
Global lowerlimit As Quad
Global currentdegree As DWord
Global currentroot As DWord
Global ti, tf As DWord
Global tt As Ext
Global filestring As String = "Lovecraft5.txt"
Global outfile As DWord

'~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

Function TBMain()
Local i As DWord

currentroot = 0
lowerlimit = 0

ti = GetTickCount
setinitialpolycoeff()

Do
currentroot += 1
currentdegree = %maxdegree - currentroot + 1
calcderivcoeff()
calcnextroot()
If currentdegree = 1 Then Exit Do
calcreducedpolycoeff()
setnewpolycoeff()
Loop

tf = GetTickCount
calctotaltime()

outfile = FILE_Open(filestring, "output")
writeoutput()
FILE_Close(outfile)

WaitKey
End Function

'~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

Function calcnextroot()

'#########################

'As the script is currently constructed, the polynomial must have all integer
'roots.
'
'Using the method of this script, we always find the leftmost root (the lowest in
'value).
'
'At any point, the value of the derivative of a function, is a straight line that
'is tangent to the curve of the function.
'
'Here, we follow the derivative line to the x axis. The intersection gives the
'new value of x.
'
'For any right triangle, with legs "x" and "y"; and "t" defined as the acute angle
'between the hypotenuse and x,
'
'y = x * tan(t).
'
'Therefore,
'
'x = y / tan(t).
'
'Below, "calcpolyval(oldx)" corresponds to y, and "calcderivval(oldx)" corresponds to
'tan(t).
'
'When the absolute value of, the new value of x minus the old value of x, is very
'small, we round it to the nearest integer, and that integer is the root.

'#########################

Local oldx, newx As Ext
Local count As DWord

oldx = lowerlimit

Do
newx = oldx - calcpolyval(oldx) / calcderivval(oldx)

If Abs(newx - oldx) < %eps Then
root(currentroot) = newx
lowerlimit = root(currentroot)
Exit Do
End If

oldx = newx
Loop

End Function

'~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

Function setinitialpolycoeff()

polycoeff(1) = -06040853943350632690, +06379933987025808250, -00340211189835457006
polycoeff(4) = +00001132119607924548, -00000000973447643102, +00000000000000000001

currentdegree = %maxdegree
End Function

'~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

Function zeropolycoeff()
Local i As DWord
For i = 0 To %maxdegree
setpolycoeff(i, 0)
Next
End Function

'--------------------------------------------------------------

Function setpolycoeff(degree As DWord, coeff As Quad)
polycoeff(degree + 1) = coeff
End Function

'--------------------------------------------------------------

Function getpolycoeff(degree As DWord) As Quad
Return polycoeff(degree + 1)
End Function

'--------------------------------------------------------------

Function calcpolyval(x As Ext) As Ext
Local total As Ext = 0
Local degree As Ext
Local i As DWord
For i = 0 To currentdegree
degree = i
total += getpolycoeff(i) * x ^ degree
Next
Return total
End Function

'~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

Function zeroderivcoeff()
Local i As DWord
For i = 0 To %maxdegree
setderivcoeff(i, 0)
Next
End Function

'--------------------------------------------------------------

Function calcderivcoeff()

'#########################

'For the polynomial form,
'
'c * x ^ n,
'
'where c and n are integers,
'
'the derivative is,
'
'n * c * x ^ (n - 1).

'#########################

Local i As DWord
zeroderivcoeff()
For i = 1 To currentdegree
setderivcoeff(i - 1, getpolycoeff(i) * i)
Next
End Function

'--------------------------------------------------------------

Function setderivcoeff(degree As DWord, V As Quad)
derivcoeff(degree + 1) = V
End Function

'--------------------------------------------------------------

Function getderivcoeff(degree As DWord) As Quad
Return derivcoeff(degree + 1)
End Function

'--------------------------------------------------------------

Function calcderivval(x As Ext) As Ext
Local total As Ext = 0
Local i As DWord
For i = 0 To currentdegree - 1
total += getderivcoeff(i) * x ^ i
Next
Return total
End Function

'~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

Function zeroreducedpolycoeff()
Local i As DWord
For i = 0 To %maxdegree
setreducedpolycoeff(i, 0)
Next
End Function

'--------------------------------------------------------------

Function calcreducedpolycoeff()

'#########################

'A polynomial's degree, is the value of its highest exponent.
'
'An nth degree polynomial, has n roots.
'
'Here we are reducing the polynomial (lowering its degree by one).
'
'Say, we begin with the polynomial,
'
'x^3 - 6*x^2 + 11*x - 6.
'
'I know its roots are x = 1, 2, 3.
'
'How do I know?
'
'Because, I constructed it that way.
'
'I wanted the roots to be x = 1, 2, 3, so I began with,
'
'(x - 1) * (x - 2) * (x - 3).
'
'Notice, that if x = 1, 2, or 3, then one of the three factors is 0, and
'therefore, the result is 0.
'
'When the three factors are multiplied together, we get,
'
'x^3 - 6*x^2 + 11*x - 6.
'
'Using the method of this script, we always find the leftmost root, so, the first
'root found, would be 1, using,
'
'x^3 - 6*x^2 + 11*x - 6 = 0.
'
'Then, we reduce the polynomial from degree 3, to degree 2, using what is called
'synthetic division.
'
'We perform the division,
'
'(x^3 - 6*x^2 + 11*x - 6) / (x - 1),
'
'and we get,
'
'x^2 -5*x + 6,
'
'which is equal to,
'
'(x - 2) * (x - 3).
'
'Then, we find the left root of,
'
'x^2 -5*x + 6 = 0,
'
'which is 2.
'
'Then, we reduce the polynomial from degree 2, to degree 1, by performing the
'division,
'
'(x^2 -5*x + 6) / (x - 2),
'
'and we get,
'
'x - 3.
'
'Then, we solve,
'
'x - 3 = 0,
'
'and we are done, we have found all three roots.
'
'So, the method works by alternately finding the leftmost root, and then reducing
'the polynomial by one degree, until all of the roots have been found.

'#########################

Local i, cd, cr As DWord
Local pc, rpc As Quad

cd = currentdegree
cr = currentroot

zeroreducedpolycoeff()

For i = 1 To cd
pc = getpolycoeff(cd - i + 1)
rpc = getreducedpolycoeff(cd - i + 1)
setreducedpolycoeff(cd - i, pc + rpc * root(cr))
Next

End Function

'--------------------------------------------------------------

Function setreducedpolycoeff(degree As DWord, V As Quad)
reducedpolycoeff(degree + 1) = V
End Function

'--------------------------------------------------------------

Function getreducedpolycoeff(degree As DWord) As Quad
Return reducedpolycoeff(degree + 1)
End Function

'--------------------------------------------------------------

Function setnewpolycoeff()
Local i As DWord

zeropolycoeff()
For i = 0 To currentdegree - 1
setpolycoeff(i, getreducedpolycoeff(i))
Next
End Function

'~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

Function calctotaltime()
tt = (tf - ti)/1000.0
End Function

'~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

Function writeoutput()
Local s1, s2, s3 As String
Local i As Integer

s1 = "Root Value"
writeoutputline(s1)
writeoutputline("")

For i = 1 To %maxdegree
s1 = Format$(i, "00")
s2 = s1 & " "
s1 = Format$(root(i), "000000000000")
s3 = s2 & s1
writeoutputline(s3)
Next

writeoutputline("")

s1 = "elapsed time = "
s2 = Format$(tt, "0000.000")
s3 = s1 & s2
s1 = " seconds"
s2 = s3 & s1
writeoutputline(s2)

End Function

'--------------------------------------------------------------

Function writeoutputline(s4 As String)
Console_WriteLine(s4)
FILE_LinePrint(outfile, s4)
End Function

'~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

[danbaron]

You can solve a degree 16 polynomial, if in the program you change,

%maxdegree = 5,

to

%maxdegree = 16.

And in function "setinitialpolycoeff()", you use these coefficients, -->

polycoeff(01) = +020922789888000, -070734282393600, +102992244837120, -087077748875904
polycoeff(05) = +048366009233424, -018861567058880, +005374523477960, -001146901283528
polycoeff(09) = +000185953177553, -000023057159840, +000002185031420, -000000156952432
polycoeff(13) = +000000008394022, -000000000323680, +000000000008500, -000000000000136
polycoeff(17) = +000000000000001

I tried solving a degree 20 polynomial. It didn't work. My guess is that the
reason is, we are losing too much precision. I think we would need more
significant digits, than Quad and Ext provide.

Dan