# Implicitisation et parametrisation des courbes rationnelles

## Calcul de l'indice:

Calcul d'indice, l'algorithme:

TraInd(P):=block(
       [P11,P12,P21,P22,G1,G2,G,index],
       P11:num(P[1]), 
       P12:denom(P[1]), 
       P21:num(P[2]), 
       P22:denom(P[2]),
       G1:ev(P11,t=s)*P12-ev(P12,t=s)*P11, 
       G2:ev(P21,t=s)*P22-ev(P22,t=s)*P21,
       G:gcd(G1,G2), 
       index:hipow(G,t),
       if index=1 then print("la parametrisation est propre")
                  else index);

TraInd([t,t**2]);

TraInd([t**2,t**8]);

P:[(t**2-1)*t/(t**4-t**2+1),(t**2-1)*t**2/(t**6-3*t**4+3*t**2-1-2*t**3)];

TraInd(P);

L'indice est 2, la parametrisation n'est pas propre

## Implicitisation

Implicit(P):=block(
[P11,P12,P21,P22,H1,H2,R,dR],
P11:num(P[1]), P12:denom(P[1]), P21:num(P[2]), P22:denom(P[2]),
H1:P12*x-P11, H2:P22*y-P21,
R:resultant(H1,H2,t), 
if hipow(R,x)#0 then 
    (dR:diff(R,x),
     return(rat(rat(R/gcd(R,dR,x)))))
else
    (dR:diff(R,y),
    return(rat(rat(R/gcd(R,dR,y)))))
);

Implicit([(t**2-1)*t/(t**4-t**2+1),(t**2-1)*t**2/(t**6-3*t**4+3*t**2-1-2*t**3)]);

Implicit([t,1]); Implicit([1,t]);

## Reparametrisation propre

Reparam(P):=block(
       [P11,P12,P21,P22,G1,G2,G,index,alpha,beta,R,g1,g2],
       P11:num(P[1]), 
       P12:denom(P[1]), 
       P21:num(P[2]), 
       P22:denom(P[2]),
       G1:ev(P11,t=s)*P12-ev(P12,t=s)*P11, 
       G2:ev(P21,t=s)*P22-ev(P22,t=s)*P21,
       G:gcd(G1,G2), 
print(G),
       index:hipow(G,t), 
       if index=1 then 
           return(P),
       alpha:0, beta:0,
       while ev(G,s=alpha)=0 do alpha:alpha+1,
       while ev(G,s=alpha,t=beta)=0 do beta:beta+1,
       R:ev(G,s=alpha)/ev(G,s=beta),
print(R),
       g1:expand(Implicit([P11/P12,R])), g2:expand(Implicit([P21/P22,R])),
print([g1,g2]),
       return(rat(ev([-coeff(g1,x,0)/coeff(g1,x,1),-coeff(g2,x,0)/coeff(g2,x,1)],y=t)))
       );

Tests

Rep:Reparam([t**2,t**4]);

I:Implicit(Rep);

PP:Reparam(P);

is(expand(Implicit(PP))=expand(Implicit(P)));

## Comparaisons de traces de courbes

On reprend l'exemple P precedent (precedant ?)

wxplot2d([['parametric, (t**2-1)*t/(t**4-t**2+1), (t**2-1)*t**2/(t**6-3*t**4+3*t**2-1-2*t**3), [t, -1000, 1000], [nticks, 30000]]], [x,-0.7,0.7])$

wxplot2d([['parametric, t/(t**2+1), t**2/(-2*t**3+1), [t, -1000, 100], [nticks, 30000]]], [x,-0.7,0.7])$

load(implicit_plot);

R:Implicit(P); wximplicit_plot(R,[x,-0.6,0.6],[y,-2,2]);

## Parametrisation par des droites

### Calcul du degre d'un polynome bivarie

degre(f):=block(
    [d],
    d:0,
    for i:0 thru hipow(f,x) do
        d:max(d,i+hipow(coeff(f,x,i),y)), 
    return(d)
);

degre(y^8+x**5+3+x**2*y**2);

### Calcul de la composante homogene de degre k

homogene(f,k):=block(
    [ff,r],
    ff:expand(f),
    r:sum(coeff(coeff(ff,x,i),y,k-i)*x**i*y**(k-i),i,0,k),
    return(r)
);

### Le cas des coniques

on suppose f absolument irreductible de degre 2

ParaConic(f):=block(
    [a,g,g1,g2],
    if degre(f)#2 then
        (print("la courbe n'est pas une conique"),
        return()),
    if hipow(f,y)=0 or hipow(f,x)=0 then 
        (print("la courbe de f est union de deux droites et n'est pas rationnelle"),
        return()),
    alpha:0,
    while hipow(ev(f,x=alpha),y)=0 do 
        alpha:alpha+1,
    a:solve(ev(f,x=alpha),y),
print(a),
    g:expand(ev(f,x=x+alpha,y=y+args(a[1])[2])),
print(g),
    [g1,g2]:[homogene(g,1),homogene(g,2)],
    [g1,g2]:ev([g1,g2],x=1,y=t),
print([g1,g2]),
    return(rat([alpha-g1/g2,args(a[1])[2]-t*g1/g2]))
);

f:x**2+y**2+3; C:ParaConic(f);

is(Implicit(C)=f);

### Le cas des courbes avec une singularite de multiplicte d-1

On cherche si une telle singularite existe et on la retourne

Multiplicite(f,p,prec):=block([d,g,j],
    d:degre(f),
    g:ev(f,x=x+p[1],y=y+p[2]),
    print(),
    j:0,
    while sqrt(sum(
        abs(coeff(coeff(homogene(g,j),y,i),x,j-i))**2,i,0,j)
        ) <10**(-prec) do j:j+1,
    return(j)
);

Multiplicite(x**7+y**3,[0,0],20);

Singula(f,prec):=block(
    [d,dfx,dfy,S,L,p],
    d:degre(f),
    dfx:diff(f,x),
    dfy:diff(f,y),
    S:solve([f=0,dfx=0,dfy=0],[x,y]),
print(length(S)),
    m:0,
    L:[],
    for i:1 thru length(S) do
        (p:[args(S[i][1])[2],args(S[i][2])[2]],
         L:append(L,[[p,Multiplicite(f,p,prec)]])
        ),
    return(L)
    );

f:expand((x**3+y**3)*((x-1)**2+(y-2)**3));

Singula(f,10);

SingMax(f,prec):=block([S,d,bool,j],
    S:Singula(f,prec),
    d:degre(f),
    bool:false,
    j:0,
    while bool=false and j<length(S) do
        (j:j+1,
         if S[j][2]=d-1 then bool:true),
    if bool=true then 
        return(S[j])
    else return(false)
);

SingMax(x**5+3*x**4*y+y**6-x**6,5);

Calcul des multiplicites via les reultants ? ....

factor(ev(f,y=0));

r1:resultant(f,diff(f,x),x); r2:resultant(f,diff(f,y),x);

factor(gcd(r1,r2));

p1:resultant(f,diff(f,x),y); p2:resultant(f,diff(f,y),y); p3:resultant(diff(f,x),diff(f,y),y);

factor(p3);

factor(resultant(f,x-1,x));

factor(gcd(p1,p3));

Singula(fex,10);

[x=0,y=0];

args(%);

L'algo complet

ParamDroite(f,prec):=block(
    [g,g1,g2,d],
    S:SingMax(f,prec),
    if S=false then return("pas de parametrisation par droite"),
    p:S[1],
    d:S[2]+1,
    g:expand(ev(f,x=x+p[1],y=y+p[2])),
    [g1,g2]:[homogene(g,d-1),homogene(g,d)],
    [g1,g2]:ev([g1,g2],x=1,y=t),
    return(rat([p[1]-g1/g2,p[2]-t*g1/g2]))
);

Exemple 1 (2 tangentes complexes et une reelle)

g:expand(ev(30*y**3-25*x*y**2+x**3*y+x**4-6*y**4+6*x**2*y,x=x+1,y=y-2));

Singula(g,1);

Pg:ParamDroite(g,5);

expand(Implicit(Pg));

wximplicit_plot(g, [x,-1.5,0], [y,1.5,2.5]);

wxplot2d([['parametric, -(6*t^4-30*t^3+25*t^2-7*t-1)/(6*t^4-t-1),(42*t^4-25*t^3+6*t^2-2*t-2)/(6*t^4-t-1), [t, -100, 100], [nticks, 30000]]], [x,-1.5,0.7],[y,1,3])$

Exemple 2 (cf Winkler et al.)

gg:1+x-15*x**2-29*y**2+30*y**3-25*x*y**2+x**3*y+35*x*y+x**4-6*y**4+6*x**2*y;

wximplicit_plot(gg, [x,-100,100], [y,-100,100]);

Pgg:ParamDroite(gg,5);

is(gg=expand(Implicit(Pgg)));

## Parametrisation par les adjoints

### Un exemple

F:y**4+(x-1)*y**2+(x+2)*x**2;

wximplicit_plot(F, [x,-5,5], [y,-5,5]);

dFx:diff(F,x);

dFy:diff(F,y);

H:solve([F=0,dFx=0,dFy=0],[x,y]);

Ici, les singularités sont de multiplicité 2, les courbes adjointes sont donc simplement les courbes passant par les singularités.

kill(P,a,b,c,d,e,f);

P:a*x**2+b*x*y+c*y**2+d*x+e*y+f;

Les polynomes P doivent passer par les points singuliers de F (l'ensemble H) et on choisit de les faire aussi passer par (-2,0) qui est un point de la courbe definie par F.

solve([ev(P,x=0,y=0)=0,ev(P,x=-1,y=1)=0,ev(P,x=-1,y=-1)=0,ev(P,x=-2,y=0)=0],[a,b,c,d,e,f]);

Cela signifie a=c=d/2, b=e et f=0. Par multiplication par une constante, les polynomes definissent les memes courbes, on multiplie donc par 1/a et on pose t=b/a.

P_t:x**2+t*x*y+y**2+2*x+t*y;

R_x:factor(resultant(F,P_t,x));

R_y:factor(resultant(F,P_t,y));

Par les resultants, on voit les coordonnees des points d'intersection entre les courbes de F et de P_t. Le resultant par rapport a x montre la coordonnee en y, et inversement pour R_y.
Les facteurs sans t correspondent aux points que l'on a fixe precedemment (les points de H et (-2,0).

P:[3*t**2-t**4-2,t**3-2*t];

Implicit(P);

### Test pour savoir si la fonction n'a bien que des singularités ordinaires

SinguOrdi(f,prec):= block(
    [L,a,p,m,g,gm,gmy,gmx,dgmy,dgmx],
    L:Singula(f,prec),
    a:true,
    for i:1 thru length(L) do
                (p:L[i][1],
                 m:L[i][2],
                 g:ev(f,x=x+p[1],y=y+p[2]),
                 gm:homogene(g,m),
                 gmy:ev(gm,x=1),
                 dgmy:diff(gmy,y),
                 gmx:ev(gm,y=1),
                 dgmx:diff(gmx,x),
                 if (hipow(gcd(gmy,dgmy),y)#0 or hipow(gcd(gmx,dgmx),x)#0) then a:false),
    return(a)
    );

SinguOrdi(F,10);

### Calcul de l'ensemble des adjoints d'une courbe

Adjoint(f,prec):=block(
    [L,d,G,A,a,eq,p,m,H,S,P],
    L:Singula(f,prec),
    d:degre(f),
    G:0,
    A:[],
    eq:[],
    for i:0 thru d-2 do
           (for j:0 thru i do
               (A:append(A,[concat(a,j,i-j)]),
                G:G+concat(a,j,i-j)*x**j*y**(i-j)
               )
           ),
    for i:1 thru length(L) do
               (p:L[i][1],
                m:L[i][2],
                H:ev(G,x=x+p[1],y=y+p[2]),
                for j:0 thru m-2 do
                        eq:append(eq,[homogene(H,j)=0])
               ),
    S:solve(eq,A),
print(S),
    P:0,
    for i:0 thru d-2 do
           for j:0 thru i do
                    P:P+args(S[1][i+j+1])[2]*x**j*y**(i-j),
    return(P)
    );

A:Adjoint(F,10);