[R] Solving 100th order equation

Sun May 25 03:10:59 CEST 2008

Oops!  Actually spectacularly bad.  I didn't see the positive exponent!

Curiously, there appears to be just one bad apple:

> sort(Mod(p(z)))
  [1] 1.062855e-10 1.062855e-10 1.328999e-10 1.328999e-10 2.579625e-10
  [6] 2.579625e-10 3.834721e-10 3.834721e-10 3.875288e-10 3.875288e-10
 [11] 5.287459e-10 5.287459e-10 5.306241e-10 5.306241e-10 6.678424e-10
 [16] 6.678424e-10 6.876300e-10 6.876300e-10 8.519089e-10 8.519089e-10
 [21] 9.345531e-10 9.345531e-10 9.369015e-10 9.369015e-10 9.789510e-10
 [26] 9.789510e-10 1.041601e-09 1.041601e-09 1.046996e-09 1.046996e-09
 [31] 1.149073e-09 1.149073e-09 1.209875e-09 1.209875e-09 1.232793e-09
 [36] 1.232793e-09 1.244167e-09 1.244167e-09 1.388979e-09 1.388979e-09
 [41] 1.556354e-09 1.556354e-09 1.596424e-09 1.596424e-09 1.610661e-09
 [46] 1.610661e-09 1.722577e-09 1.722577e-09 1.727842e-09 1.727842e-09
 [51] 1.728728e-09 1.728728e-09 1.769644e-09 1.769644e-09 1.863983e-09
 [56] 1.863983e-09 1.895296e-09 1.895296e-09 1.907509e-09 1.907509e-09
 [61] 1.948662e-09 1.948662e-09 2.027021e-09 2.027021e-09 2.061063e-09
 [66] 2.061063e-09 2.116735e-09 2.116735e-09 2.185764e-09 2.185764e-09
 [71] 2.209158e-09 2.209158e-09 2.398479e-09 2.398479e-09 2.404217e-09
 [76] 2.404217e-09 2.503279e-09 2.503279e-09 2.643436e-09 2.643436e-09
 [81] 2.654788e-09 2.654788e-09 2.695735e-09 2.695735e-09 2.921933e-09
 [86] 2.921933e-09 2.948185e-09 2.948185e-09 2.953596e-09 2.953596e-09
 [91] 3.097433e-09 3.097433e-09 3.420593e-09 3.420593e-09 3.735880e-09
 [96] 3.735880e-09 4.190042e-09 4.554964e-09 4.554964e-09 1.548112e+15
>  

Bill Venables
CSIRO Laboratories
PO Box 120, Cleveland, 4163
AUSTRALIA
Office Phone (email preferred): +61 7 3826 7251
Fax (if absolutely necessary):  +61 7 3826 7304
Mobile:                         +61 4 8819 4402
Home Phone:                     +61 7 3286 7700
mailto:Bill.Venables at csiro.au
http://www.cmis.csiro.au/bill.venables/ 

-----Original Message-----
From: r-help-bounces at r-project.org [mailto:r-help-bounces at r-project.org] On Behalf Of Venables, Bill (CMIS, Cleveland)
Sent: Sunday, 25 May 2008 10:58 AM
To: ggrothendieck at gmail.com; shubhak at ambaresearch.com
Cc: r-help at stat.math.ethz.ch; murdoch at stats.uwo.ca; p.dalgaard at biostat.ku.dk
Subject: [ExternalEmail] Re: [R] Solving 100th order equation

> library(PolynomF)
> x <- polynom()
> p <- x^100 - 2*x^99 + 10*x^50 + 6*x - 4000
> z <- solve(p)
> z
  [1] -1.0741267+0.0000000i -1.0739999-0.0680356i -1.0739999+0.0680356i -1.0655699-0.1354644i
  [5] -1.0655699+0.1354644i -1.0568677-0.2030274i -1.0568677+0.2030274i -1.0400346-0.2687815i
  ...
 [93]  1.0595174+0.2439885i  1.0746575-0.1721335i  1.0746575+0.1721335i  1.0828132-0.1065591i
 [97]  1.0828132+0.1065591i  1.0879363-0.0330308i  1.0879363+0.0330308i  2.0000000+0.0000000i
> 

Now to check how good they are:

> 
> range(Mod(p(z)))
[1] 1.062855e-10 1.548112e+15
> 

Not brilliant, but not too bad. 

Bill Venables
CSIRO Laboratories
PO Box 120, Cleveland, 4163
AUSTRALIA
Office Phone (email preferred): +61 7 3826 7251
Fax (if absolutely necessary):  +61 7 3826 7304
Mobile:                         +61 4 8819 4402
Home Phone:                     +61 7 3286 7700
mailto:Bill.Venables at csiro.au
http://www.cmis.csiro.au/bill.venables/ 

-----Original Message-----
From: r-help-bounces at r-project.org [mailto:r-help-bounces at r-project.org] On Behalf Of Gabor Grothendieck
Sent: Sunday, 25 May 2008 10:09 AM
To: Shubha Vishwanath Karanth
Cc: r-help at stat.math.ethz.ch; Duncan Murdoch; Peter Dalgaard
Subject: Re: [R] Solving 100th order equation

Actually maybe I was premature. It does not handle the polynomial
I tried it on in the example earlier in this thread but it does seem to work
with the following very simple polynomials of order 100.  At any rate it
would not take long to try it on the real problem and see.

> Solve(x^100 - 1, x)
[1] "Starting Yacas!"
expression(list(x == complex_cartesian(cos(pi/50), sin(pi/50)),
    x == complex_cartesian(cos(pi/25), sin(pi/25)), x ==
complex_cartesian(cos(3 *
        pi/50), sin(3 * pi/50)), x == complex_cartesian(cos(2 *
        pi/25), sin(2 * pi/25)), x == complex_cartesian(cos(pi/10),
        (root(5, 2) - 1)/4), x == complex_cartesian(cos(3 * pi/25),
        sin(3 * pi/25)), x == complex_cartesian(cos(7 * pi/50),
        sin(7 * pi/50)), x == complex_cartesian(cos(4 * pi/25),
        sin(4 * pi/25)), x == complex_cartesian(cos(9 * pi/50),
        sin(9 * pi/50)), x == complex_cartesian(cos(pi/5), sin(pi/5)),
    x == complex_cartesian(cos(11 * pi/50), sin(11 * pi/50)),
    x == complex_cartesian(cos(6 * pi/25), sin(6 * pi/25)), x ==
        complex_cartesian(cos(13 * pi/50), sin(13 * pi/50)),
    x == complex_cartesian(cos(7 * pi/25), sin(7 * pi/25)), x ==
        complex_cartesian(cos(3 * pi/10), sin(3 * pi/10)), x ==
        complex_cartesian(cos(8 * pi/25), sin(8 * pi/25)), x ==
        complex_cartesian(cos(17 * pi/50), sin(17 * pi/50)),
    x == complex_cartesian(cos(9 * pi/25), sin(9 * pi/25)), x ==
        complex_cartesian(cos(19 * pi/50), sin(19 * pi/50)),
    x == complex_cartesian((root(5, 2) - 1)/4, sin(2 * pi/5)),
    x == complex_cartesian(cos(21 * pi/50), sin(21 * pi/50)),
    x == complex_cartesian(cos(11 * pi/25), sin(11 * pi/25)),
    x == complex_cartesian(cos(23 * pi/50), sin(23 * pi/50)),
    x == complex_cartesian(cos(12 * pi/25), sin(12 * pi/25)),
    x == complex_cartesian(0, 1), x == complex_cartesian(-cos(12 *
        pi/25), sin(12 * pi/25)), x == complex_cartesian(-cos(23 *
        pi/50), sin(23 * pi/50)), x == complex_cartesian(-cos(11 *
        pi/25), sin(11 * pi/25)), x == complex_cartesian(-cos(21 *
        pi/50), sin(21 * pi/50)), x == complex_cartesian(-((root(5,
        2) - 1)/4), sin(2 * pi/5)), x == complex_cartesian(-cos(19 *
        pi/50), sin(19 * pi/50)), x == complex_cartesian(-cos(9 *
        pi/25), sin(9 * pi/25)), x == complex_cartesian(-cos(17 *
        pi/50), sin(17 * pi/50)), x == complex_cartesian(-cos(8 *
        pi/25), sin(8 * pi/25)), x == complex_cartesian(-cos(3 *
        pi/10), sin(3 * pi/10)), x == complex_cartesian(-cos(7 *
        pi/25), sin(7 * pi/25)), x == complex_cartesian(-cos(13 *
        pi/50), sin(13 * pi/50)), x == complex_cartesian(-cos(6 *
        pi/25), sin(6 * pi/25)), x == complex_cartesian(-cos(11 *
        pi/50), sin(11 * pi/50)), x == complex_cartesian(-cos(pi/5),
        sin(pi/5)), x == complex_cartesian(-cos(9 * pi/50), sin(9 *
        pi/50)), x == complex_cartesian(-cos(4 * pi/25), sin(4 *
        pi/25)), x == complex_cartesian(-cos(7 * pi/50), sin(7 *
        pi/50)), x == complex_cartesian(-cos(3 * pi/25), sin(3 *
        pi/25)), x == complex_cartesian(-cos(pi/10), (root(5,
        2) - 1)/4), x == complex_cartesian(-cos(2 * pi/25), sin(2 *
        pi/25)), x == complex_cartesian(-cos(3 * pi/50), sin(3 *
        pi/50)), x == complex_cartesian(-cos(pi/25), sin(pi/25)),
    x == complex_cartesian(-cos(pi/50), sin(pi/50)), x == -1,
    x == complex_cartesian(-cos(pi/50), -sin(pi/50)), x ==
complex_cartesian(-cos(pi/25),
        -sin(pi/25)), x == complex_cartesian(-cos(3 * pi/50),
        -sin(3 * pi/50)), x == complex_cartesian(-cos(2 * pi/25),
        -sin(2 * pi/25)), x == complex_cartesian(-cos(pi/10),
        -((root(5, 2) - 1)/4)), x == complex_cartesian(-cos(3 *
        pi/25), -sin(3 * pi/25)), x == complex_cartesian(-cos(7 *
        pi/50), -sin(7 * pi/50)), x == complex_cartesian(-cos(4 *
        pi/25), -sin(4 * pi/25)), x == complex_cartesian(-cos(9 *
        pi/50), -sin(9 * pi/50)), x == complex_cartesian(-cos(pi/5),
        -sin(pi/5)), x == complex_cartesian(-cos(11 * pi/50),
        -sin(11 * pi/50)), x == complex_cartesian(-cos(6 * pi/25),
        -sin(6 * pi/25)), x == complex_cartesian(-cos(13 * pi/50),
        -sin(13 * pi/50)), x == complex_cartesian(-cos(7 * pi/25),
        -sin(7 * pi/25)), x == complex_cartesian(-cos(3 * pi/10),
        -sin(3 * pi/10)), x == complex_cartesian(-cos(8 * pi/25),
        -sin(8 * pi/25)), x == complex_cartesian(-cos(17 * pi/50),
        -sin(17 * pi/50)), x == complex_cartesian(-cos(9 * pi/25),
        -sin(9 * pi/25)), x == complex_cartesian(-cos(19 * pi/50),
        -sin(19 * pi/50)), x == complex_cartesian(-((root(5,
        2) - 1)/4), -sin(2 * pi/5)), x == complex_cartesian(-cos(21 *
        pi/50), -sin(21 * pi/50)), x == complex_cartesian(-cos(11 *
        pi/25), -sin(11 * pi/25)), x == complex_cartesian(-cos(23 *
        pi/50), -sin(23 * pi/50)), x == complex_cartesian(-cos(12 *
        pi/25), -sin(12 * pi/25)), x == complex_cartesian(0,
        -1), x == complex_cartesian(cos(12 * pi/25), -sin(12 *
        pi/25)), x == complex_cartesian(cos(23 * pi/50), -sin(23 *
        pi/50)), x == complex_cartesian(cos(11 * pi/25), -sin(11 *
        pi/25)), x == complex_cartesian(cos(21 * pi/50), -sin(21 *
        pi/50)), x == complex_cartesian((root(5, 2) - 1)/4, -sin(2 *
        pi/5)), x == complex_cartesian(cos(19 * pi/50), -sin(19 *
        pi/50)), x == complex_cartesian(cos(9 * pi/25), -sin(9 *
        pi/25)), x == complex_cartesian(cos(17 * pi/50), -sin(17 *
        pi/50)), x == complex_cartesian(cos(8 * pi/25), -sin(8 *
        pi/25)), x == complex_cartesian(cos(3 * pi/10), -sin(3 *
        pi/10)), x == complex_cartesian(cos(7 * pi/25), -sin(7 *
        pi/25)), x == complex_cartesian(cos(13 * pi/50), -sin(13 *
        pi/50)), x == complex_cartesian(cos(6 * pi/25), -sin(6 *
        pi/25)), x == complex_cartesian(cos(11 * pi/50), -sin(11 *
        pi/50)), x == complex_cartesian(cos(pi/5), -sin(pi/5)),
    x == complex_cartesian(cos(9 * pi/50), -sin(9 * pi/50)),
    x == complex_cartesian(cos(4 * pi/25), -sin(4 * pi/25)),
    x == complex_cartesian(cos(7 * pi/50), -sin(7 * pi/50)),
    x == complex_cartesian(cos(3 * pi/25), -sin(3 * pi/25)),
    x == complex_cartesian(cos(pi/10), -((root(5, 2) - 1)/4)),
    x == complex_cartesian(cos(2 * pi/25), -sin(2 * pi/25)),
    x == complex_cartesian(cos(3 * pi/50), -sin(3 * pi/50)),
    x == complex_cartesian(cos(pi/25), -sin(pi/25)), x ==
complex_cartesian(cos(pi/50),
        -sin(pi/50)), x == 1))

On Sat, May 24, 2008 at 8:56 AM, Shubha Vishwanath Karanth
<shubhak at ambaresearch.com> wrote:
> Was also wondering which theoretical method is used to solve this problem?
>
> Thanks,
> Shubha Karanth | Amba Research
> Ph +91 80 3980 8031 | Mob +91 94 4886 4510
> Bangalore * Colombo * London * New York * San José * Singapore * www.ambaresearch.com
>
> -----Original Message-----
> From: Gabor Grothendieck [mailto:ggrothendieck at gmail.com]
> Sent: Saturday, May 24, 2008 6:13 PM
> To: Peter Dalgaard
> Cc: Shubha Vishwanath Karanth; r-help at stat.math.ethz.ch; Duncan Murdoch
> Subject: Re: [R] Solving 100th order equation
>
> On Sat, May 24, 2008 at 8:31 AM, Peter Dalgaard
> <p.dalgaard at biostat.ku.dk> wrote:
>> Shubha Vishwanath Karanth wrote:
>>>
>>> To apply uniroot I don't even know the interval values... Does numerical
>>> methods help me? Or any other method?
>>>
>>> Thanks and Regards,
>>> Shubha
>>>
>>> -----Original Message-----
>>> From: Duncan Murdoch [mailto:murdoch at stats.uwo.ca] Sent: Saturday, May 24,
>>> 2008 5:08 PM
>>> To: Shubha Vishwanath Karanth
>>> Subject: Re: [R] Solving 100th order equation
>>>
>>> Shubha Vishwanath Karanth wrote:
>>>
>>>>
>>>> Hi R,
>>>>
>>>>
>>>> I have a 100th order equation for which I need to solve the value for x.
>>>> Is there a package to do this?
>>>>
>>>>
>>>> For example my equation is:
>>>>
>>>>
>>>> (x^100 )- (2*x^99) +(10*x^50)+.............. +(6*x ) = 4000
>>>>
>>>>
>>>> I have only one unknown value and that is x. How do I solve for this?
>>>>
>>>>
>>>
>>> uniroot() will find one root.  If you want all of them, I don't know what
>>> is available.
>>>
>>> Duncan Murdoch
>>>
>>
>> polyroot() is built for this, but it stops at 48th degree polynomials, at
>> least as currently implemented. Not sure that it (or anything else) would be
>> stable beyond that limit. YACAS perhaps?
>>
>
> Unfortunately yacas does not seem to be able to handle it:
>
>> library(Ryacas)
>> x <- Sym("x")
>> Solve((x^100 )- (2*x^99) +(10*x^50)+(6*x ) - 4000 == 0, x)
> [1] "Starting Yacas!"
> expression(list())
>
> Simpler one works ok:
>
>> Solve(x^2 - 1, x)
> expression(list(x == 1, x == -1))
> This e-mail may contain confidential and/or privileged...{{dropped:12}}

______________________________________________
R-help at r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.