# Roots of unity

Hans Peter Würmli wurmli at freesurf.ch
Sun Jan 20 08:40:12 CET 2002

```Dear Bob

I have written a

lst RootOf(const lst & l) returning a list of the three zeros plus the lead coefficient

program for degree 3 polynomials, where the input list contains the polynomial and the indeterminate. In the sample output of the main program you will see the effect of precision with Digits=yyy that "zero" is represented by some long digit sequence times E-xxx that you felt some time ago was a bug of GiNaC, but isn't.

I don't know how one can force GiNaC to simplify expressions. Say, I used the programm to find symbolically the three roots x1, x2, x3 and reassembled the polynomonial by expanding leadcoefficient*(x-x1)*(x-x2)*(x-x3). To take an example:

p(x) = 1 + x + x^3

Reassembling p from the roots and expanding produces:

p(x) = 1/2+x+x^3+(-27/2+3/2*sqrt(93))^(-1)-1/18*sqrt(93)

which is correct, if one simplifies 1/2+(-27/2+3/2*sqrt(93))^(-1)-1/18*sqrt(93) to 1

Best regards

Hans Peter

//
// Roots of third degree polynomial
// (from Kochendörfer "Einführung in die Algebra"
//

#include <iostream>
#include <ginac/ginac.h>

using namespace GiNaC;

//
// lst contains the polynomial and the indeterminate
//

lst RootOf(const lst& l){

//
// some checks - throwing an error would be appropriate
//

if (l.nops()!=2) return lst();
ex poly=l.op(0);

//
// Cint cannot resolve the next line
//

if (!is_a<symbol>(l.op(1))) return lst();

//
// my Cint complains if I use the name "x" again, therefore "z"
//

ex z=l.op(1);
if (poly.degree(z)!=3) return lst();

ex a0=poly.coeff(z,0), a1=poly.coeff(z,1),
a2=poly.coeff(z,2),a3=poly.coeff(z,3);

symbol t("t");

//
// normalisation to get the form x^3+a1*x+a0
//

ex pol=  (poly/a3).subs(z==t-a2/(a3*3)).expand();

//
// polynomial is now p(t)=t^3+p*t+q
//

ex  q=pol.coeff(t,0), p=pol.coeff(t,1);

//
// rho^3=1, crho^3=1, the two non-trivial roots of x^3-1
//

ex rho=(-1+pow(3,numeric(1,2))*I) /2,
crho=(-1-pow(3,numeric(1,2))*I)/2;

if (p.is_zero())
return  lst((-pow(q,1/numeric(3))-save).expand(),
(-rho*pow(q,1/numeric(3))-save).expand(),
(-crho*pow(q,1/numeric(3))-save).expand(),

//
// the discriminant
//

ex D= -numeric(4)*pow(p,3)-numeric(27)*pow(q,2);

// Lagrange resolvents

ex lrho=pow(-numeric(27,2)*q+numeric(3,2)*pow(-3*D,numeric(1,2)),1/numeric(3));
ex lcrho=-3*p/lrho;
ex th=numeric(1,3);

}

//
//Some examples
//

symbol y("y"), x("x"), a0("a0"),a1("a1"),a2("a2"),a3("a3");
ex poly=1+x+pow(x,3);
ex poly1=a0+a1*x+a2*pow(x,2)+a3*pow(x,3);
ex poly2=27+pow(x,3);
lst l(poly,x), l2(poly2,x),l1(poly1,x);

void main(void){

cout << RootOf(l) << "\n\n\n";
cout << RootOf(l1) << "\n\n\n";
cout << RootOf(l2) << "\n\n\n";

ex test=RootOf(l2).op(3)*(x-RootOf(l2).op(0))*(x-RootOf(l2).op(1))*(x-RootOf(l2).op(2));
cout << l2.op(0) << " = " << test.expand() << "\n\n\n";

ex test1=RootOf(l).op(3)*(x-RootOf(l).op(0))*(x-RootOf(l).op(1))*(x-RootOf(l).op(2));
cout << l.op(0) << " = " << test1.expand() << "\n\n\n";

Digits=50;
cout << evalf(poly.subs(x==RootOf(l).op(0))) <<  "\n\n\n";
cout << evalf(poly.subs(x==RootOf(l).op(1))) <<  "\n\n\n";
cout << evalf(poly.subs(x==RootOf(l).op(2))) <<  "\n\n\n";
cout << evalf(poly2.subs(x==RootOf(l2).op(0))) <<  "\n\n\n";
cout << evalf(poly2.subs(x==RootOf(l2).op(1))) <<  "\n\n\n";
cout << evalf(poly2.subs(x==RootOf(l2).op(2))) <<  "\n\n\n";

cout << "This list should be empty = " << RootOf(lst(pow(x,3),x+y)) <<  "\n\n\n";

}

-------------- next part --------------
An embedded and charset-unspecified text was scrubbed...
Name: thirdroots.cpp
Url: http://www.cebix.net/pipermail/ginac-list/attachments/20020120/539e838d/thirdroots.cpp
```