share/sym/resolvante.mac

   1 /* CALCUL DES RESOLVANTES PRODUITS DE 1 AU DEGRE DU POLYNOMES */
   2
   3 /* AVEC LES RESULTANTS */
   4 resolvante_produit_res(p,x):=
   5    block([n,rmax,puissances,elementaires,rh,k,rf],
   6           n:hipow(p,x),
   7           rmax : binomial(n,quotient(n,2)),
   8           rh: makelist(resultant(p,y-x^i,x),i,0,rmax),
   9           puissances: makelist(maplist(lambda([pol],
  10                                               (-1)^k*coeff(pol,y,n-k)),
  11                                        expand(rh)),
  12                                k,0,n),
  13           elementaires :maplist(lambda([puissances] ,
  14                                     pui2ele(first(puissances),puissances)),
  15                                 puissances),
  16           k : 0,
  17           maplist(lambda([elem],
  18                          (deg:first(elem),
  19                           rf[k] : y^deg,
  20                           for i:1 thru deg do
  21                               (elem:rest(elem),
  22                                rf[k]:rf[k] +
  23                                        (-1)^(i)*first(elem)*y^(deg-i)),
  24                           k:k+1,
  25                           rf[k-1])),
  26                   elementaires));
  27
  28 /* AVEC LES FONCTIONS SYMETRIQUES */
  29 /* Les A[k] sont les fonctions puissances des racines de la re'solvante */
  30
  31 resolvante_produit_sym (p,x):=
  32     block([n,rmax,krmax,a,aa,resol,pui_depart],
  33        (n:hipow(p,x),
  34         krmax : if oddp(n) then 1+quotient(n,2) else quotient(n,2),
  35         a[0]:makelist(binomial(n,r),r,0,n),
  36         a[1] :cons(n,makelist(coeff(p,x,n-i)*(-1)^i/coeff(p,x,n),i,1,n)),
  37         pui_depart:ele2pui(binomial(n,krmax)*krmax,a[1]),
  38         for i:1 thru quotient(n,2)
  39            do for k:binomial(n,i-1)+1 thru binomial(n,i)
  40               do a[k] : pui2ele(n-i,makelist(part(pui_depart,r*k+1),
  41                                           r,0,n-i)),
  42         makelist((bin : binomial(n,r),
  43                   aa : pui2ele(bin,makelist(part(a[k],r+1),k,0,bin)),
  44                   resol : y^bin,
  45                   for j:bin-1 step -1 thru 0 do
  46                       (aa:rest(aa),
  47                        resol : resol + (-1)^(bin-j)*first(aa)*y^j),
  48                  resol), r,1,n)));
  49
  50 resolvante_unitaire(p,q,x):=
  51     block([aa,ele,pp,ppui,n,m,alt,resol],
  52     (p:expand(p),q:expand(q),
  53      n:hipow(p,x),
  54      ele:cons(n,makelist(coeff(p,x,n-j)*(-1)^j/coeff(p,x,n),j,1,n)),
  55      m:hipow(q,x),
  56      pp:ele2pui(n*m,ele), print(pp),
  57      ppui : expand(makelist(q^j,j,1,n)), print(ppui),
  58      ppui : makelist((aa:part(ppui,j),
  59                       aa:aa + (part(pp,1)-1)*coeff(aa,x,0),
  60                 for k:m*j step -1 thru 0 do
  61                     aa:ratsubst(part(pp,k+1),x^k,aa),aa),
  62                 j,1,n), print(ppui),
  63      pp:pui2ele(n,cons(n,ppui)),
  64      resol:y^n,alt:1,
  65      for i:n-1 step -1 thru 0 do
  66         ( pp:rest(pp),
  67           alt:-1*alt,
  68           resol:alt*first(pp)*y^i + resol),
  69      resol));
  70
  71 /*   Pour calculer prod_{1<= i<j<= 7} (x_i-x_j)   */
  72
  73 resolvante_alternee1(f,x) :=
  74     block([r,c,delta,n],
  75       n:hipow(f,x),
  76       r:resultant(f,y-diff(f,x)^2,x),
  77       c : ev(r,y=0),
  78       r:expand(subst(y=1/z,r/y^8)),
  79       delta: poly_discriminant(f,x),
  80       r:expand(delta^8*subst(z=y/delta,r)/c),
  81       ev(r,y=z^2));
  82
  83
  84 /* Calcul de $x_1x_2+x_3$ */
  85
  86
  87 resolvante_klein(polynome,x) :=
  88        block([degre,p,pui,elem,e,n],
  89        e:polynome2ele(polynome,x),
  90        n:hipow(polynome,x),
  91        degre : 3*binomial(n,3),
  92        p : rest(ele2pui(2*degre,e)),
  93        print(fait),
  94        pui: cons(degre,
  95                cons((n-2)*e[3]+(n-2)*(n-1)/2*e[2],
  96                     makelist(monterme(p,r)
  97                              + (n-2)*(p[r]^2/2-p[2*r]/2)
  98                              + (n-2)*(n-1)/2*p[r]
  99                              +3*h(p,r), r,2,degre))),
 100        print(fait),
 101        elem : rest(pui2ele(degre,pui)),
 102        print(fait),
 103        x^degre + sum((-1)^i*elem[i]*x^(degre-i),i,1,degre))$
 104
 105 monterme(p,r) := block([somme],
 106                somme:0,
 107                for i:quotient(r,2)+1 thru r-1 do
 108                  (somme : somme +
 109                       binomial(r,i)*(p[r+i]-p[i]*p[r]-p[2*i]*p[r-i]/2+
 110                           p[i]^2*p[r-i]/2 + p[2*r-i]-p[r-i]*p[r]-
 111                           p[2*r-2*i]*p[i]/2+ p[r-i]^2*p[i]/2)),
 112                somme)$
 113
 114 h(p,r) := if oddp(r) then 0
 115              else binomial(r,r/2)*(p[3*r/2]/3-p[r/2]*p[r]/2+p[r/2]^3/6)$
 116
 117 /* Calcul de $x_1x_2x_3+x_4$ */
 118
 119 resolvante_klein3(polynome,x) :=
 120        block([degre,p,lim,e,n],
 121        e:polynome2ele(polynome,x),
 122        n:hipow(polynome,x),
 123        degre : 4*binomial(n,4),
 124        p : rest(ele2pui(3*degre,e)), print(fait),
 125        p: cons(degre,cons((n-3)*e[4]+(n-3)*(n-2)*(n-1)/3*e[2],
 126              makelist(monterme3(p,r)+
 127                        + (n-3)*(p[3*r]/3-p[r]*p[2*r]/2 + p[r]^3/6) +
 128                        (n-3)*(n-2)*(n-1)/3*p[r]
 129                        +3*h3(p,r),r,2,degre))), print(fait),
 130        p : rest(pui2ele(degre,p)), print(fait),
 131        x^degre + sum((-1)^i*p[i]*x^(degre-i),i,1,degre))$
 132
 133 monterme3(p,r) :=
 134        block([somme],
 135                somme:0,
 136                for i:quotient(r,2)+1 thru r-1 do
 137                  (somme : somme +
 138                              binomial(r,i)*
 139                  (-p[2*i+r]+p[i]*p[i+r] + p[2*i]*p[r]/2 - p[i]^2*p[r]
 140          +p[3*i]*p[r-i]/3 -p[i]*p[2*i]*p[r-i]/2 + p[i]^3*p[r-i]/6
 141          -p[3*r-2*i]+p[r-i]*p[2*r-i] + p[2*r-2*i]*p[r]/2 - p[r-i]^2*p[r]
 142          +p[3*r-3*i]*p[i]/3 -p[r-i]*p[2*r-2*i]*p[i]/2 + p[r-i]^3*p[i]/6)),
 143                  somme)$
 144
 145
 146 h3(p,r) := if oddp(r) then 0
 147              else binomial(r,r/2)*
 148           (-p[2*r]/4+p[3*r/2]*p[r/2]/3+p[r]^2/8 +p[r]*p[r/2]^2/4+p[r/2]^4/24)$
 149
 150 /*  CALCUL DE LA RESOLVANTE x1x2-x3x4 */
 151 /* e EST LISTE DES FONCTIONS SYMETRIQUES ELEMENTAIRES [n,e1,...,en]
 152    DU POLYNOME QUE L'ON TRANSFORME DEGRE 420 pour degre 8 */
 153   /* (NON TESTEE) */
 154
 155 resolvante_vierer (polynome,x) :=
 156       block( [degre,p,e,n],
 157         e:polynome2ele(polynome,x),
 158         n:hipow(polynome,x),
 159         degre : 6*binomial(n,4),
 160         p : rest(ele2pui(2*degre,e)),
 161         p: cons(degre,
 162                 makelist(if oddp(r) then 0
 163                             else 2*termevierer(p,r)
 164                                  + (-1)^r*(n-3)*(n-2)*(p[r]^2-p[2*r])/2
 165                                  + 6*(-1)^(r/2)*binomial(r,r/2)*
 166                (-p[2*r]/4+p[3*r/2]*p[r/2]/3+p[r]^2/8 - p[r]*p[r/2]^2/4
 167                    +p[r/2]^4/24),r,1,degre)),
 168        p:rest(pui2ele(degre,p)),
 169        x^degre + sum((-1)^i*p[i]*x^(degre-i),i,1,degre))$
 170
 171 termevierer(p,r) := block([somme],
 172             somme:0,
 173            for i:r/2+1 thru r-1 do
 174                somme : somme +
 175                       (-1)^i*binomial(r,i)*(-3*p[2*r]/2+p[r-i]*p[r+i]
 176                         +p[2*i]*p[2*(r-i)]/4 -p[r-i]^2*p[2*i]/4
 177                         -p[i]^2*p[2*(r-i)]/4
 178                          +p[i]*p[2*r-i]+p[r]^2/2
 179                      -p[i]*p[r-i]*p[r] +p[i]^2*p[r-i]^2/4),
 180            somme)$
 181
 182 /*  CALCUL DE LA RESOLVANTE x1x2 + x3x4  de D_8 */
 183 /* e EST LISTE DES FONCTIONS SYMETRIQUES ELEMENTAIRES [n,e1,...,en]
 184    DU POLYNOME QUE L'ON TRANSFORME DEGRE 210 pour degre le degre 8 */
 185
 186 resolvante_diedrale (polynome,x) :=
 187       block( [degre,p,e,n],
 188         e:polynome2ele(polynome,x),
 189         n:hipow(polynome,x),
 190         degre : 3*binomial(n,4),
 191         p : rest(ele2pui(2*degre,e)),
 192         p: cons(degre,
 193                 makelist(termediedral(p,r)
 194                          + (n-3)*(n-2)*(p[r]^2-p[2*r])/4
 195                                  + if oddp(r) then 0
 196                                       else 3*binomial(r,r/2)*
 197                (-p[2*r]/4+p[3*r/2]*p[r/2]/3+p[r]^2/8 - p[r]*p[r/2]^2/4
 198                    +p[r/2]^4/24),r,1,degre)),
 199        p:rest(pui2ele(degre,p)),
 200        x^degre + sum((-1)^i*p[i]*x^(degre-i),i,1,degre))$
 201
 202 termediedral(p,r) := block([somme],
 203             somme:0,
 204            for i:quotient(r,2)+1 thru r-1 do
 205                somme : somme +
 206                       binomial(r,i)*(-3*p[2*r]/2+p[r-i]*p[r+i]
 207                         +p[2*i]*p[2*(r-i)]/4 -p[r-i]^2*p[2*i]/4
 208                         -p[i]^2*p[2*(r-i)]/4
 209                          +p[i]*p[2*r-i]+p[r]^2/2
 210                      -p[i]*p[r-i]*p[r] +p[i]^2*p[r-i]^2/4),
 211            somme)$
 212
 213
 214 /* RESOLVANTE BIPARTITE : x1x2x3..x(n/2) + x(n/2+1)....xn
 215    NE FONCTIONNE QUE SI n EST PAIR
 216    SE CALCULE EN 3 OU 5 SECONDES POUR LE DEGRE 6
 217    SE CALCULE EN 4 mn POUR LE DEGRE 8 */
 218
 219 resolvante_bipartite(polynome,x) :=
 220       block( [degre,pui_pol,elem_pol,n,pui_resol,elem_resol],
 221         elem_pol:polynome2ele(polynome,x),
 222         n:hipow(polynome,x),
 223         degre : binomial(n-1,n/2-1),
 224         pui_pol : ele2pui(4*degre,elem_pol),
 225         pui_resol : [],
 226         for r:degre step -1 thru 1 do
 227              (print(r),
 228              pui_resol : cons(pui(pui_pol,
 229                                    polynome_bipartite(r),
 230                                    makelist(concat(a,j),j,1,n)),
 231                                pui_resol)),
 232         elem_resol : pui2ele(degre,cons(degre,pui_resol)),
 233         ele2polynome(elem_resol,y))$
 234
 235 polynome_bipartite(r) :=
 236            block([s,borne],
 237            s : 0,
 238            borne:quotient(r,2)+remainder(r,2),
 239            for i:borne thru r do
 240                 s : s +
 241                     binomial(r,i)*prod(concat(a,j)^i,j,1,n/2)
 242                                  *prod(concat(a,j)^(r-i),j,n/2+1,n),
 243            if oddp(r) then  s
 244               else s +(degre-1)*binomial(r,r/2)*prod(concat(a,j)^(r/2),j,1,n))$
 245
 246 /* resolvante_bipartite(x^6+2*x^2+2,x);
 247
 248    resolvante_bipartite(x^6+2*x^3-2,x);  */
 249
 250