4
$\begingroup$

How does Mathematica's RootOfUnityQ function work?

That is, how does Mathematica know if a number is a root of unity?

Example:

Let $x = \frac{1-i \sqrt{2+\sqrt{5}}}{1+i \sqrt{2+\sqrt{5}}}$.

Then

In[1]:= RootOfUnityQ[x]
Out[1]= False

Wrong Answers:

This bahavior is not the same as

In[2]:= Reduce[Abs[x] == 1 && Simplify[Element[Arg[x]/(2 Pi), Rationals]]]
Out[2]= ArcTan[(2 Sqrt[2 + Sqrt[5]])/(1 + Sqrt[5])]/\[Pi] \[Element] Rationals
Improve this question
$\endgroup$
1

2 Answers 2

Reset to default
6
$\begingroup$

As long as we're giving methods that mimic the behavior, here is a quick implementation based on Kronecker's theorem. 

RoUQ[u_] :=
 If[! (Abs[N[u]] == 1), False,
  If[! AlgebraicIntegerQ[u], False,
   (f = MinimalPolynomial[u, x];
    n = Exponent[f, x];
    cf = CoefficientList[f, x]/Coefficient[f, x^n]; 
    M = Table[
      If[j == i + 1, 1,
       If[i == n, -cf[[j]], 0]], {i, 1, n}, {j, 1, n}];
    lambda = First[Eigenvalues[M, 1]];
    Abs[N[lambda]] == 1)]]

The sole reason for building the matrix $M$ is that Eigenvalues[M,1] is guaranteed to return the largest eigenvalue of $M$ (in absolute value), while I don't know how to tell FindRoot that it must give me the largest root.

J.M. below suggests just finding all the roots and taking the largest one. Based on a tiny amount of testing, I think that this is slower. Here is what I did:

(* Both methods will be tested on the same set of polynomials *)
data=Table[Sum[RandomInteger[100] x^j, {j,0,10}], {i,1,10^4}];

(* Find all roots and extract the largest *)
BigRoot1[f_]:=Max[Abs[x] /. NSolve[f==0,x]]

(* Find the largest eigenvalue of the companion matrix. *)
BigRoot2[f_]:=
(n = Exponent[f, x];
    cf = CoefficientList[f, x]/Coefficient[f, x^n]; 
    M = Table[
      If[j == i + 1, 1,
       If[i == n, -cf[[j]], 0.0]], {i, 1, n}, {j, 1, n}];
    lambda = First[Eigenvalues[M, 1]];
    N[Abs[lambda]])                                                               

Timing[Map[BigRoot1, data]][[1]]

10.9893

Timing[Map[BigRoot2, data]][[1]]

3.00254

I wouldn't take these results too seriously, because I'm sure both implementations do the list processing inefficiently, but it suggests that I'm not crazy to use Eigenvalues[].

Improve this answer
$\endgroup$
4
  • $\begingroup$ Well, you could use x /. First @ Solve[p[x] == 0, x] to generate the roots, and then Pick[roots, #, Max[#]] &[Abs[roots]] to get the largest of them all in magnitude. But using the companion matrix is cool anyway. :) $\endgroup$
    – J. M.'s missing motivation ♦
    Commented Aug 15, 2012 at 23:37
  • $\begingroup$ @J.M. It strikes me as wasteful to compute all those roots which will be discarded. I've added some tests above to support this claim. $\endgroup$
    – David E Speyer
    Commented Aug 16, 2012 at 0:14
  • $\begingroup$ Hmm, right. Then, I suppose if one wants to use FindRoot[], there's always the option of trying to find the smallest root of $x^n p\left(\frac1{x}\right)$ and reciprocating that... the only reason I was somewhat lukewarm to using the companion matrix is that it looks uncomfortable to use $O(n^2)$ storage and $O(n^3)$ effort for just a single root... $\endgroup$
    – J. M.'s missing motivation ♦
    Commented Aug 16, 2012 at 1:45
  • $\begingroup$ @J.M. Yeah, memory usage is a good point. I guess the next question then is how often FindRoot[], starting at 0, fails to find the smallest root. A good first test would be p[x]=x^10-x^9+x^7-x^6+x^5-x^4+x^3-x+1, or perhaps this after the change of variable x /. E^(2 Pi I RandomReal[]) x. This has two roots at 0.85 and 1.176 and the other 8 are all on the unit circle with roughly equal spacing, so a random walk from 0 is very likely to hit one of the roots on the circle before the small root. I imagine FindRoot[] does something smarter than a random walk, though. $\endgroup$
    – David E Speyer
    Commented Aug 16, 2012 at 12:13
4
$\begingroup$

If Abs[x]==1 and Element[Arg[x]/(2 Pi), Rationals], then the number is a root of unity

$$ \arg[x]= 2 \pi \frac{p}{q}$$ $p,q\in \mathbb{Z}$

and $$ |x|=1 $$

$$ \to x^q=1 $$

because
$ |x^q|=|x|^q=1^q=1 $ and
$ \arg[x^q] = \arg[x] q = 2 \pi p $

So perhaps this mimics the behaviour

ru = TrueQ[Abs[#] == 1 && Simplify@Element[Arg[#]/(2 Pi), Rationals] ]&;
Improve this answer
$\endgroup$
4
  • $\begingroup$ It does not mimic the behavior. I updated the question to discuss your answer. $\endgroup$
    – Tyson Williams
    Commented Aug 1, 2012 at 19:48
  • $\begingroup$ @TysonWilliams in your counterexample, it's just a matter of adding TrueQ. There are probably other counterexamples where this doesn't work, but please post them. I'm editing my answer $\endgroup$
    – Rojo
    Commented Aug 1, 2012 at 20:08
  • $\begingroup$ @R.M I am not looking for (numerically) approximate answers. I want to know with certainty. $\endgroup$
    – Tyson Williams
    Commented Aug 1, 2012 at 20:10
  • $\begingroup$ Ok, suppose that this is the definition of RootOfUnityQ. The purpose of my question was to be able to verify Mathematica's answer for myself. So my question now becoems this: how does Mathematica compute ArcTan[Sqrt[1/2 (1+Sqrt[5])]]/\[Pi]\[Element]Rationals? $\endgroup$
    – Tyson Williams
    Commented Aug 1, 2012 at 21:34

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.