nolabels:true;
display2d:false;linel:12500;fpprec:16;algexact:true;
declare(const,constant);
logexpand:false;
logabs:true;
/* mettre les valeurs absolues quand on intègre en ln*/
lognegint:true;
/* pour que les ln de négatifs soient repérés comme non réels..*/
/*logarc:true;*/
/* pour convertir les arcsh..s'applique aussi aux trigos inverses!!*/
argsinh(xx):=logarc(asinh(xx));
argtanh(xx):=logarc(atanh(xx));
argcosh(xx):=logarc(acosh(xx));
/*pour éviter les solutions implicites*/
solveexplicit:true;
solveradcan : true;


maplist(lambda([tt],texput(tt,concat("\\text{",tt,"}"))),[A,B,C,F,G,H,I,J,K,L,M,N,P,Q,R,S,T,U,V,W,Y,Z,a,b,c,e,f,g,h,i,j,k,l,m,n,p,q,r,s,t,u,v,w,x,y,z]);

maplist(lambda([tt],texput(concat("t",tt),concat("t_",tt))),[A,B,C,F,G,H,I,J,K,L,M,N,P,Q,R,S,T,U,V,W,Y,Z,a,b,c,e,f,g,h,i,j,k,l,m,n,p,q,r,s,t,u,v,w,x,y,z]);

maplist(lambda([tt],texput(concat("x",tt),concat("x_",tt))),[A,B,C,F,G,H,I,J,K,L,M,N,P,Q,R,S,T,U,V,W,Y,Z,a,b,c,e,f,g,h,i,j,k,l,m,n,p,q,r,s,t,u,v,w,y,z]);/*pas xx*/

maplist(lambda([tt],texput(concat("y",tt),concat("y_",tt))),[A,B,C,F,G,H,I,J,K,L,M,N,P,Q,R,S,T,U,V,W,Y,Z,a,b,c,e,f,g,h,i,j,k,l,m,n,p,q,r,s,t,u,v,w,x,y,z]);


texput(pi,"{\\it pi}");
degre(ee):=hipow(expand(ee),xx);
realonly:true;
algebraic:true;

Realpart(pp):=realpart(ev(pp,infeval));

Imagpart(pp):=imagpart(ev(pp,infeval));
sen(tt):=sin(tt);
asen(tt):=asin(tt);

conjuguate(pp):=Realpart(pp)-%i*Imagpart(pp);

solution(expr,inc):=block([ss],ss:algsys([expr],[inc]),map(lambda([tt],rhs(tt[1])),ss));

solutions(expr,inc):=block([ss], ratvars(), ss:solve(expr, inc), ratvars(xx), ss:sublist(ss, lambda([tt], lhs(tt)=inc)), map(lambda([tt], ratsimp(rhs(tt))),ss));
/*or constantp(tt) retire car pb avec ln de nb negat*/
part_rec(tt)::=block([simp:false],if atom(tt)  then [] elseif member(op(tt),[sqrt,log,"/",tan,"^"]) then cons(tt,flatten(map(part_rec,args(tt)))) else flatten(map(part_rec,args(tt))));

conditions(tt,BorneInf,BorneSup):= block([LL:(xx>BorneInf),MM:(xx<BorneSup),ee],  assume(LL,MM), ee:
if operatorp(tt,sqrt)  then is(part(tt,1)>=0)or is(QuadraSignum(part(tt,1))=1) 
elseif operatorp(tt,"/")then is(abs(denom(tt))>0) or not is(QuadraSignumStrict(denom(tt))=0)  
elseif operatorp(tt,log) then is(part(tt,1)>0)or is (QuadraSignumStrict(part(tt,1))=1)  
elseif operatorp(tt,tan) then is(abs(cos(part(tt,1)))>0) 
elseif operatorp(tt,"^") then ((integerp(part(tt,2)) or is(part(tt,1)>=0)))or is(QuadraSignum(part(tt,1))=1)  
else true, forget(LL,MM), ee);

existe(ee):=block([boo:part_rec(ee)],emptyp(boo) or apply("and",map(lambda([tt], conditions(tt,minf,inf)),boo)));
ratvars(xx);
declare(xx,mainvar );
elimine_reel(ee):=block([aa:if op(ee) = "-" then part(ee,1) else ee],if op(aa) = "*" then map(lambda([tt],if constantp(tt) then 1 else tt),aa) else aa);
zero(expr,liste):=block([aa], for aa in liste do (if errcatch(expand(subst(aa,xx,expr)))=[0] then return(aa)));
islinear(expr):= block ([cc],  cc: bothcoef (expand(expr), xx), is (freeof (xx, cc)  and cc[1]#0));
isquadratic(expr):=block([cc:coeff(expand(expr),xx,0)],is(freeof(xx,cc) and islinear((expr-cc)/xx)));
ishomographic(ee):=(islinear(ratnumer(ee)) or freeof(xx,ratnumer(ee))) and islinear (ratdenom(ee));
isexponentiel(pp):=block([aa,bb],aa:diff(pp,xx),bb:aa/diff(aa,xx),is(not(bb=0)) and freeof(xx,bb));
isradical(pp):=islinear(1/diff(pp,xx)^2);
islog(pp):=islinear(1/diff(pp,xx));
isTrigInv(pp):=islinear(tan(pp)) or islinear(cos(pp)) or islinear(sin(pp));
NumFoncRef(ee):=   if islinear(ee) then 0    elseif isquadratic(ee) then 1     elseif ishomographic(ee)then 2     elseif isexponentiel(ee)then 3     elseif isradical(ee)then 4     elseif islog(ee)then 5    else -1;
canonique(ee):=block([cc,ff,dd],cc:ratcoeff(ee,xx,2),ff:-ratcoeff(ee,xx,1)/cc/2,dd:subst([xx=ff],ee),cc*(xx-ff)^2+dd);

QuadraSignum(ee):=if isquadratic(ee) then block([cc,ff,dd],cc:ratcoeff(ee,xx,2),ff:-ratcoeff(ee,xx,1)/cc/2,dd:subst([xx=ff],ee),if cc*dd>=0 then if cc<0 then -1 else 1 else 0) else 0; 

QuadraSignumStrict(ee):=if isquadratic(ee) then block([cc,ff,dd],cc:ratcoeff(ee,xx,2),ff:-ratcoeff(ee,xx,1)/cc/2,dd:subst([xx=ff],ee),if cc*dd>0 then if cc<0 then -1 else 1 else 0) else 0; 

FormVariFoncRef(ee,ii):=if  ii=1 then canonique(ee) elseif ii=2  then partfrac(ee) else expand(ee) ;
coefHomo(ee):=block([cc:solution(denom(rat(ee)),xx)[1]],limit(ee*(xx-cc),xx,cc));
CoefVariFoncRef(ee,ii):=if ii<2 then coeff(expand(ee),xx,ii+1) elseif ii=2 then coefHomo(ee) else limit(ee,xx,temp[2],minus)-limit(ee,xx,temp[1],plus);
SimpliContext(expr,BorneInf,BorneSup):=block([LL:(xx>=BorneInf),MM:xx<=BorneSup,ee],assume(LL,MM),ee:ev(expr),forget(LL,MM),ee);
FactElem(expr,liste):=block([aa], if (aa:zero(expr,liste)) = done then expr else if atom(expr) then expr else if operatorp(expr,"^") then FactElem(part(expr,1),[aa])^part(expr,2) else if isquadratic(expr) then block([ee:expand(expr)], coeff(ee,xx,2)*(xx-aa)*(xx + aa + coeff(ee, xx, 1)/coeff(ee, xx, 2))) else expr);
/* factorisation des entiers*/
MaxFactor(expr,BorneInf,BorneSup,liste):=block([ee],ee:SimpliContext(expr,BorneInf,BorneSup),ee:factor(ee),if atom(ee) or numberp(expr) then ee else (if operatorp(ee,"-") then -MaxFactor(-ee,BorneInf,BorneSup,liste) else (if operatorp(ee,"/") then MaxFactor(num(ee),BorneInf,BorneSup,liste)/MaxFactor(denom(ee),BorneInf,BorneSup,liste) else (if operatorp(ee,"*") then map(lambda([tt],FactElem(tt,liste)),ee) else FactElem(ee,liste)))));
FactorSansConstante(ee):=
 block(([tt:factor(ee)]),
 if not atom(tt) and operatorp(tt,"*")then
     map(lambda([u],if freeof(xx,u) then 1 else u),tt)
 else tt);
NumAppendCasyop(liste,valeur):= if (Imagpart(valeur)=0)  then NumAppendCasyopExact(bfloat(liste),bfloat(valeur)) else -1; 
NumAppendCasyopExact(liste,valeur):=
  block([indices:sublist_indices (liste, lambda ([xx], xx=valeur))],
    if length(indices)>0 then (-indices[1]-1)
  else
   block([listp:sort (append([valeur],delete(ß,liste)),lambda([uu,vv],is(uu<vv)))],
          indices: sublist_indices (listp, lambda([x],x=valeur)),
     if ((indices[1]=1) and (length(listp)>1) and (unknown=is(valeur<listp[indices[1]+1]))) then
         -1
     else
     if (indices[1]>1)then
      sublist_indices(liste, lambda([x],x=listp[indices[1]-1]))[1]+1
     else 1));
AppendCasyop(liste,valeur,ii):=
  block([list1: makelist(liste[k],k,1,ii),list2:makelist(liste[k],k,ii+1,length(liste))],
  append(list1,[valeur],list2));


PointRepere(x,y):=x+%i*y;
PointLibre(nom):=
   block([paramx,paramy],
    paramx:concat('x,nom),
    paramy:concat('y,nom),texput(paramx,concat("x_{",nom,"}")),texput(paramy,concat("y_{",nom,"}")), paramx+%i*paramy);
PointReporte(Pref,Pvect1,Pvect2,distreport):=
  expand(Pref+(Pvect2-Pvect1)/d(Pvect2,Pvect1) * distreport);
PointTransformeAffine(pt, a11, a12, a21, a22, b1, b2) :=
      block([xx, yy], xx:Realpart(pt), yy:Imagpart(pt),
      resTA:expand(a11*xx + a12*yy + b1) + %i*expand(a21*xx + a22*yy + b2),
      resTA);
PointLibreSurDroitePassant2Points(nom,nn,pp):=
block([param],
       param:concat('t,nom),texput(param,concat("t_{",nom,"}")), 
       expand(nn*(1-param)+pp*param));
PointLibreSurDroitePointDroiteOLD(nom,pt1,dd):=
block([param,aa,bb],
       param:concat('t,nom),
       aa:dd(1),bb:dd(%i),
       expand(pt1+param*(-bb+%i*aa)/sqrt(aa^2+bb^2)));
PointLibreSurDroitePointDroiteOLD2(nom,pt1,aa,bb):=
block([param],
       param:concat('t,nom),
       DISTANCE_FONCTION_VAR:aa^2 + bb^2,
       assume(DISTANCE_FONCTION_VAR>0),
       assume(expand(DISTANCE_FONCTION_VAR)>0),
       expand(pt1+param*(-bb+%i*aa)/sqrt(DISTANCE_FONCTION_VAR)));
PointLibreSurDroitePointDroite(nom,pt1,aa,bb):=
block([param],
       param:concat('t,nom),texput(param,concat("t_{",nom,"}")), 
       expand(pt1+param*(-bb+%i*aa)));
PointLibreSurDemiDroite(nom,nn,pp):=
block([param],
       param:concat('t,nom),texput(param,concat("t_{",nom,"}")), 
       assume(param>0),
       expand(nn*(1-param)+pp*param));
PointLibreSurSegment(nom,nn,pp):=
block([param],
       param:concat('t,nom),texput(param,concat("t_{",nom,"}")), 
       assume(param>0),
       assume(param<1),
       expand(nn*(1-param)+pp*param));
milieu(mm,nn):=(mm+nn)/2;
PointLibreSurCercle(nom,nn,dd):=
block([param],
   	param:concat('t,nom),texput(param,concat("t_{",nom,"}")), 
        assume(param>-%pi),
   	assume(param<%pi),
   	nn +dd*exp(param*%i));
xAbs(mm,nn):=Realpart(nn-mm);
yAbs(mm,nn):=Imagpart(nn-mm);

DroiteTransformeeAffine(a11, a12, a21, a22, u, v, PointA, pp):=
 ([nu, nv, cc],
  nu : -v*a21 + u*a22,
  nv : v*a11 - u*a12,
  cc : -(nu*Realpart(PointA) + nv*Imagpart(PointA)),
  nu*Realpart(pp) + nv*Imagpart(pp) + cc);

d(mm,nn):=block(DISTANCE_FONCTION_VAR:xAbs(mm,nn)^2 + yAbs(mm,nn)^2,
                  assume(DISTANCE_FONCTION_VAR>0),
                  assume(expand(DISTANCE_FONCTION_VAR)>0),
                  factor(sqrt(DISTANCE_FONCTION_VAR)));
InterDroites(d1,d2):=
block([xx,yy,sol],
 sol:solve([d1,d2],[xx,yy]),
TRACE_D1:d1,
TRACE_D2:d2,
TRACE_SOL:sol,
TRACE_SUBST:subst(sol,xx+%i*yy));

InterLignesCercles(dd,cc,ListeSubs,ListeNoms,coordNumPtsInter):=
block([epsilon:0.2,xx,yy,sol,sol1,sol2,NomInstancie,SolInstancie1,NomInstancie2],sol:solve([dd,cc],[xx,yy]),ETAPE:1,sol1:subst(part(sol,1),xx+%i*yy),ETAPE:2,v1:part(ListeNoms,1),if length(sol)=1 then v1::sol1 else block(v2:part(ListeNoms,2),sol2:subst(part(sol,2),xx+%i*yy),ETAPE:3,SolInstancie1:float(subst(ListeSubs,''sol1)),ETAPE:4,if ((abs(part(coordNumPtsInter,1)-Realpart(SolInstancie1))<epsilon) and (abs(part(coordNumPtsInter,2)-Imagpart(SolInstancie1))<epsilon))        then block( CAS:1,v2::sol2,v1::sol1) else block( CAS:2,v2::sol1,v1::sol2),[[v1=sol1],[v2=sol2]]));


InterCercles(ce1,ce2,ListeSubs,ListeNoms,coordNumPtsInter):=InterLignesCercles(ce2-ce1,ce1,ListeSubs,ListeNoms,coordNumPtsInter);

InterCourbeDroiteCasSolve(expr):=
  block(if not(freeof(x,expr)) then SOL_DEPEND_DE_X
        else if listofvars(expr)=[] then if not(Imagpart(float(expr))=0) then SOL_IMAGINAIRE
                                                                  else block(NBR_SOL_POSSIBLE:NBR_SOL_POSSIBLE+1,SOL_POSSIBLE)
                                    else block(NBR_SOL_POSSIBLE:NBR_SOL_POSSIBLE+1,SOL_POSSIBLE));
IntersectionCourbeDroitePointsSolve(Droite,aa,bb,cc,Courbe):=
block([ii],
      DR:Droite,CO:Courbe,
      if bb=0 then block(CAS:DROITE_VERTICALE,sol:solve(Droite,x),NBR_SOL_POSSIBLE:0,
                        xsol:subst(sol,x),ysol:subst(sol,Courbe),PTS:[InterCourbeDroiteCasSolve(xsol),xsol,ysol])
      else block(CAS:DROITE_NON_VERTICALE,NBR_SOL_POSSIBLE:0,
	         EQ:-(aa*x+cc)/bb=Courbe,XSOLU:solve(EQ,x),PTS:[],
                 for ii:1 thru length(XSOLU) do
                      block(XSOLUI:part(XSOLU,ii),PTS:cons(InterCourbeDroiteCasSolve(part(XSOLUI,2)),append(subst(XSOLUI,[x,-(aa*x+cc)/bb]),PTS)))
                ),
      [COMMANDE_SOLVE,NBR_SOL_POSSIBLE,CAS,PTS]);
InterCourbeDroiteCasAlgSys(expr):=
  block(if not(freeof(x,expr)) then SOL_DEPEND_DE_X
        else if listofvars(expr)=[] then if not(Imagpart(float(expr))=0) then SOL_IMAGINAIRE
                                                                  else block(NBR_SOL_POSSIBLE:NBR_SOL_POSSIBLE+1,SOL_POSSIBLE)
                                    else block(NBR_SOL_POSSIBLE2:NBR_SOL_POSSIBLE2+1,SOL_POSSIBLE));
IntersectionCourbeDroitePointsAlgsys(Droite,aa,bb,cc,Courbe):=
block([ee],
      if bb=0 then block(CAS2:DROITE_VERTICALE,sol:solve(Droite,x),NBR_SOL_POSSIBLE2:0,
                        xsol:subst(sol,x),ysol:subst(sol,Courbe),PTS2:[InterCourbeDroiteCasAlgSys(xsol),xsol,ysol])
      else block(CAS2:DROITE_NON_VERTICALE,
                EQ2:-(aa*x+cc)/bb=Courbe,XSOLU2:algsys([EQ2],[x]),PTS2:[],NBR_SOL_POSSIBLE2:0,
                for ii:1 thru length(XSOLU2) do
                      block(XSOLUI2:part(XSOLU2,ii,1),PTS2:cons(InterCourbeDroiteCasAlgSys(part(XSOLUI2,2)),append(subst(XSOLUI2,[x,-(aa*x+cc)/bb]),PTS2)))
                 ),
          [COMMANDE_ALGSYS,NBR_SOL_POSSIBLE2,CAS2,PTS2]
     );
demoivre :true;
Im(x):=imagpart(x);Re(x):=realpart(x);
offset(param,pp,qq):=block([Xcentre:subst(param=%pi/2,Realpart(pp)),Ycentre:subst(param=0,Imagpart(pp))],if qq=Xcentre+%i*Ycentre then 0 else atan2(Imagpart(qq)-Ycentre,Realpart(qq)-Xcentre));
ListCoeffSymetrieDroite(dr):=block([aaa,bbb,ccc],ccc:dr(0),aaa:dr(1)-ccc,bbb:dr(%i)-ccc,map(lambda([tt],rhs(tt)),algsys([aaa*(aa+1)+bbb*bb,aaa*bb -bbb*(aa-1),
   if (bbb=0)then ccc*(aa-1)-aaa*cc
    else
     bbb*cc-ccc*bb,
 if (bbb=0)then ccc*bb-aaa*dd
    else
    bbb*dd+ccc*(aa+1)],
  [aa,bb,cc,dd])[1]));

pente(point,var,val):=block([dd,dx,dy,ptx,pty],dd:trigsimp(diff(point,var)),dx:Realpart(dd),dy:Imagpart(dd),ptx:Realpart(subst(val,var,point)),pty:Imagpart(subst(val,var,point)),if dx=0 then [true,ptx, 0,inf] elseif dy=0 then [true,0, pty,0] else[is(equal(diff(dy/dx,var),0)),0, pty-ptx*dy/dx,subst(val,var,dy/dx)] );

centre(point,var,val):=block([dd1,dd2,dx,dy,ddx,ddy,rr,pt],dd1:trigsimp(diff(point,var)),dd2:trigsimp(diff(dd1,var)),dx:Realpart(dd1),dy:Imagpart(dd1),ddx:Realpart(dd2),ddy:Imagpart(dd2),rr:trigsimp((dy^2+dx^2)^(3/2)/(dx*ddy-dy*ddx)),pt: trigsimp(point-rr*(dy-%i*dx)/sqrt(dy^2+dx^2)), [is(equal(diff(rr,var),0)) ,Realpart(pt), Imagpart(pt),abs(subst(val,var,rr))] );

 PointExiste(pp):=block([points:sublist(values,lambda([tt],is(charat(string(tt),1) = "z")))],points:makelist(expand(ev(points[i],infeval)),i,1,length(points)-1),member(expand(ev(pp,infeval)),points));

/* il faut éviter les identificaterus de paramètres. Donc les droites sont DD(pp) et donc ici il faut le paramètre ppp*/
 DroiteCercleExiste(ppp):=block([droitecercle:sublist(functions,lambda([tt],is(substring(string(tt),1,3) = substring(string(ppp),1,3)))),ee],droitecercle:map(lambda([tt],part(tt,0)),droitecercle),droitecercle:map(fundef,droitecercle),droitecercle:map(lambda([tt],part(tt,2)),droitecercle),droitecercle:subst([pp = x+%i*y],droitecercle),droitecercle:ev(droitecercle,infeval),ee:droitecercle[length(droitecercle)],droitecercle:makelist(constantp(ratsimp(droitecercle[i]/ee)),i,1,length(droitecercle)-1),apply("or",droitecercle));

Sce(ll,fonc):=sum((fonc(Realpart(ll[ii]))-Imagpart(ll[ii]))^2,ii,1,length(ll));

/* mean( [zp1,zp2,zp3])*/
mean(lz) := apply("+", args(lz)) / length(lz);

varx(zz):=block([me:apply("+",args(realpart(zz)))/length(zz)],ratsimp(apply("+",args((realpart(zz)-me)^2))/length(zz)));

/* Covariance */
 covxy(zz):=block([mex:apply("+",args(realpart(zz)))/length(zz),mey:apply("+",args(imagpart(zz)))/length(zz)],ratsimp(apply("+",args((realpart(zz)-mex)*(imagpart(zz)-mey)))/length(zz)));

regress(zz):=covxy(zz)/varx(zz);


boucle():= for i:0 thru 10 do i:i-1;

FINI:OUIVRAI;