1

I'm fairly new to Mathematica. In the past, my usage has mostly been limited to solving the occasional equation, making some plots, and working with small scaled statistics. None of these have been requiring regarding performance so I typically did not give much thought to optimizing my code.

Right now, I am working on a project that requires a large number of polynomial roots. I have written the following code to find these roots for polynomials with coefficients (1) and (-1).

fileloc = FileNameJoin[{NotebookDirectory[], "data"}];
out = OpenWrite [fileloc];

maxDeg = 12;
f[c_, x_, m_] = c*x^m;
coefficients = {-1, 1};

(* Making List of Functions *)
polyList = coefficients;
n = 1;
While[n <= maxDeg, 
  Do[Do[polyList = Append[polyList, j + f[i, x, n]], {i, 
     coefficients}], {j, polyList}]; n++];

(* Solving and Printing *)
Do[Write[out, N[Solve[i == 0, x]]], {i, polyList}];
Close[out];

I have run this script up to a maximum degree of 12. Doing that took roughly 7 hours to complete but it did not create enough roots for my needs. Could someone give me some suggestions as to where I should start to improve my code? Moreover, are there any good resources that could tell me some of the differences between similar built-in-functions (ie Solve, NSolve, FindRoot)?

As an aside, I have noticed in Parallel Kernel Status, the four cores assigned to my Local kernel are all idle while my script is running. I am not sure how to interpret this but it seemed somewhat odd to me.

Frank Wang
  • 43
  • 4
  • Do you want imaginary and real roots? – 2012rcampion Apr 18 '15 at 17:00
  • Yes, I want both imaginary and real roots. – Frank Wang Apr 18 '15 at 17:02
  • (1) Have you checked the source of the bottleneck? I ask because the quadratic complexity of the mode of generating them leads me to suspect that that's where the main problem lies. At maxDeg=9 it took 28 seconds on my machine. For maxDeg=10 it took around 5 minutes. I expect roughly a factor of 9 for each increment (factor of 3 growth, squared due to the AppendTo complexity). – Daniel Lichtblau Apr 19 '15 at 21:06
  • (2) You are also doing too much work by a factor of 2. You can insist that the constant term be 1 in p[x] because case of -1 will show up as the negative version of same (that is, -p[x]`). (This was also noted in a response by @2012rcampion.) – Daniel Lichtblau Apr 19 '15 at 21:08
  • (3) As for faster generation, could use this: Timing[maxDeg = 12; polys = 1 + Flatten[Outer[List, Apply[Sequence, Table[{1, 0, -1}, {maxDeg}]]], maxDeg - 1].x^Range[maxDeg];] – Daniel Lichtblau Apr 19 '15 at 21:09

1 Answers1

6

This function prints out the solutions to all polynomials of degree d with coefficients of unit magnitude.

printAllSolutions[d_] := 
 Do[
  With[{p = FromDigits[c~Prepend~1, x]}, 
   Print[Expand[p] -> (x /. NSolve[p == 0, x])]
  ],
  {c, Tuples[{-1, +1}, d]}
 ]

I made one optimization, which is that the leading coefficient is always 1. If it is allowed to be -1, then for each positive polynomial we will have a corresponding negative polynomial with the same roots. This cuts our work in half.

Instead of generating the list iteratively like you do, I generate the lists of coefficcents directly with Tuple, then use FromDigits to convert to a polynomial. I use NSolve instead of N@Solve, which affords a significant speedup.

2012rcampion
  • 7,851
  • 25
  • 44