//
// 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 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,
    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(), 
		lead);

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

  }