// gcc -Wall -g poly.c -std=c99 -lm
#include<stdio.h>   // printf
#include<stdlib.h>  // exit
#include<math.h>    // rintl
#define dmax 100
typedef long double poly[dmax+1];
int deg(poly a)
{ int i=dmax;
  while(i>=0 && !a[i]) i--;    // dmax-d°a < dmax
  return i;
}
void aff(poly a)
{ int nul=1, i;
  for(i=0;i<=dmax;i++)
  { long double x=a[i];
    if(!x) continue;
    if(x>0 && !nul) printf("+");
    if(x<0) printf("-"), x=-x;
    nul=0;
    if(x>1e9 || x<1e-2) printf("%Le",x);
    else if(rintl(x)!=x) printf("%Lf",x);
    else if(x>1 || !i) printf("%.0Lf",x);
    switch(i)
    { case 0: break;
      case 1: printf("X"); break;
      case 2: printf("X²"); break;
      default: printf("X^%d",i);
    }
  }
  if(nul) printf("0");
  printf("\n");
}
void zer(poly a)                 { for(int i=dmax;i>=0;i--) a[i]=0        ; } // a=0     dmax
void cop(poly a, poly b)         { for(int i=dmax;i>=0;i--) b[i]=a[i]     ; } // b=a     dmax
void add(poly a, poly b, poly c) { for(int i=dmax;i>=0;i--) c[i]=a[i]+b[i]; } // c=a+b   dmax
void sub(poly a, poly b, poly c) { for(int i=dmax;i>=0;i--) c[i]=a[i]-b[i]; } // c=a-b   dmax

void mul(poly a, poly b, poly c) // c=a*b
{ int da=deg(a), db=deg(b), i, j;                                    //  dmax
  if(da+db>dmax) printf("mul %d+%d>%d\n",da,db,dmax), exit(1);
  zer(c);                                                            //  dmax
  for(i=0;i<=da;i++) if(a[i]) for(j=0;j<=db;j++) c[i+j]+=a[i]*b[j];  //  d°b x nombre de monômes de A < d°axd°b < dmax²
}
void divi(poly a, poly b, poly q, poly r)
{ int db=deg(b), i, j;                                               //  dmax
  if(db<0) printf("division par 0\n"), exit(1);
  zer(q), cop(a,r);                                                  //  dmax
  for(i=dmax;i>=db;i--) if(r[i])
  { q[i-db]=r[i]/b[db];                                              //  max(d°a-d°b+1,0)
    for(j=db;j--;) r[i-db+j]-=q[i-db]*b[j];                          //  max(d°a-d°b+1,0) d°b
    r[i]=0;
  }
}
void pgcd(poly a, poly b, poly c)   // Le temps de calcul suppose d°a>d°b, ce qui vrai lors des appels récursifs
{ if(deg(b)<0) cop(a,c); else       // dmax <=  d°a dmax
  { poly q, r;
    divi(a,b,q,r),                  // (d°a-d°b+1)d°b + dmax < (d°a)²-(d°b)² + dmax(d°a-d°b)
    pgcd(b,r,c);                    //  (d°b)²+d°b dmax
  }
}                                   // le temps total est donc inférieur à   (d°a)²+d°a dmax <= 2dmax²
void pgcd2(poly a0, poly b0, poly c)  // Version en O(max(d°a,d°b)²+dmax)
{ poly a1, b1;                        // On économise le temps dmax(d°b) car on ne recalcule pas d°b 
  long double q, *x, *a=a1, *b=b1;
  cop(a0,a), cop(b0,b);
  int da=deg(a), db=deg(b), i, j;
  for(;db>=0;) if(!b[db]) db--;
  else if(db>da) x=a, a=b, b=x, i=da, da=db, db=i;
  else if(!a[da]) da--;
  else for(q=a[i=da]/b[j=db],a[da--]=0;j;) a[--i]-=q*b[--j];
  cop(a,c);
}
int main()
{ poly a, b, c, q, r;
  zer(a);   aff(a);
  a[3]=2.5; aff(a);
  a[5]=-4;  aff(a);
  a[0]=1;   aff(a);
  a[7]=-1;  aff(a);
  printf("d°=%d\n",deg(a));
  mul(a,a,b); aff(b);
  add(a,b,c); aff(c);
  divi(c,a,q,r); aff(q); aff(r);
  a[0]-=1; aff(a);
  divi(c,a,q,r); aff(q); aff(r);
  
  zer(a); a[60]=a[30]=a[0]=1; aff(a);
  zer(b); b[66]=b[33]=b[0]=1; aff(b);
  pgcd2(a,b,c); aff(c);
  pgcd (a,b,c); aff(c);

  return 0;
}