Roots of unity

Hans Peter Würmli wurmli at
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 (!=3) return lst();

  ex a0=poly.coeff(z,0), a1=poly.coeff(z,1),

  symbol t("t");

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

  ex save=a2/(a3*3), lead=a3;
  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,

  if (p.is_zero()) 
    return  lst((-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);

  return lst(th*(lrho+lcrho)-save, th*(crho*lrho+rho*lcrho)-save, th*(rho*lrho+crho*lcrho)-save, lead);

//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";
  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

More information about the GiNaC-list mailing list