# makes a polynom from the listof simples permutations Lpoly: [[size,number of permutation of this size],...] makeNP:=proc(Lpoly) convert([seq(c[2]*T^c[1],c=Lpoly)],`+`) end; # makes a system of derivatives corresponding to the system of generating functions made from Lpoly oracle_diff:=proc(Lpoly) #option trace; local f_dnp,deq; f_dnp:=convert([seq(c[2]*c[1]*T^(c[1]-1),c=Lpoly)],`+`); deq := [DT = 1/(-f_dnp-2*Cp+Cp^2+1-2*Cm+Cm^2), DCm = (1-2*Cm+Cm^2)/(-f_dnp-2*Cp+Cp^2+1-2*Cm+Cm^2), DTp = -Cp*(-2+Cp)/(-f_dnp-2*Cp+Cp^2+1-2*Cm+Cm^2), DTm = -Cm*(-2+Cm)/(-f_dnp-2*Cp+Cp^2+1-2*Cm+Cm^2), DCp = (1-2*Cp+Cp^2)/(-f_dnp-2*Cp+Cp^2+1-2*Cm+Cm^2), DNP = f_dnp/(-f_dnp-2*Cp+Cp^2+1-2*Cm+Cm^2)]; deq end: # makes a Boltzmann oracle (a function which evaluates the generating at a given point, with a given precision) # for the specification made from Lpoly gen_oracle:=proc(Lpoly,ifprint:=false) #option trace; local i,c, gfs,fnames,dim,sys,jac,jacm1,newtit,newton; gfs := [T=Z+Tp+Tm+NP, Tp=Cp^2/(1-Cp), Tm=Cm^2/(1-Cm), Cp=Z+Tm+NP, Cm=Z+Tp+NP, NP=makeNP(Lpoly)]; if ifprint then print(gfs) end; fnames:= [T, Tp, Tm, Cp, Cm, NP]; dim:=nops(fnames); sys:=[seq(op(1,eq)-op(2,eq),eq=gfs)]; jac:=convert(linalg['jacobian'](sys,fnames),Matrix); jacm1:=linalg['inverse'](jac); newtit:=convert(evalm(fnames-jacm1&*sys),list); newton:=proc(val,prec) option remember; local i,st,t,res,oldres,nb_ite,diff_res,err,oldDigits; st:=time(); res:=[0\$dim]; oldDigits:=Digits; Digits:=prec+3; for nb_ite to 1000 do oldres:=res; res:=eval(newtit,[seq(fnames[i]=res[i],i=1..dim),Z=val]); for t in res do if not type(t,numeric) then return [seq(fnames[i]=-1,i=1..dim)] end if od; diff_res:=seq(res[i]-oldres[i],i=1..dim); if convert(map(x->x<-Float(1,-prec),[diff_res]),`or`) then Digits:=oldDigits; return [seq(fnames[i]=-1,i=1..dim)] fi; err:=max(diff_res); if nb_ite>5 and err-1 and tmpmax then error("trop grand\n") end if; if bern(x/O_T(x))=1 then _size:=_size+1; return Z elif bern(O_Tp(x)/(O_Tp(x)+O_Tm(x)+O_NP(x)))=1 then return genTp(x,max) elif bern(O_Tm(x)/(O_Tm(x)+O_NP(x)))=1 then return genTm(x,max) else return genNP(x,max) fi; end; genTp:=proc(x,max) global _size; local g; if _size>max then error("trop grand\n") end if; g:=fast_geom(O_Cp(x))+2; return p(seq(genCp(x,max),i=1..g)); end; genTm:=proc(x,max) #option trace; global _size; local g; if _size>max then error("trop grand\n") end if; g:=fast_geom(O_Cm(x))+2; return m(seq(genCm(x,max),i=1..g)); end; genCp:=proc(x,max) global _size; if _size>max then error("trop grand\n") end if; if bern(x/O_Cp(x))=1 then _size:=_size+1; return Z elif bern(O_Tm(x)/(O_Tm(x)+O_NP(x)))=1 then return genTm(x,max) else return genNP(x,max) fi; end; genCm:=proc(x,max) global _size; if _size>max then error("trop grand\n") end if; if bern(x/O_Cm(x))=1 then _size:=_size+1; return Z elif bern(O_Tp(x)/(O_Tp(x)+O_NP(x)))=1 then return genTp(x,max) else return genNP(x,max) fi; end; genNP:=proc(x,max) global _size; local tmp,n,nb_arg,i; if _size>max then error("trop grand\n") end if; nb_arg:=nops(L); tmp:=O_NP(x); for i from 1 to nb_arg do n:=nops(L[i]); if bern(O_T(x)^n/tmp)=1 then return NP[L[i]](seq(genT(x,max),k=1..n)) end if; tmp:=tmp-O_T(x)^n; od; return TOTO end; proc(x,n) global _size; _size:=0; genT(x,n) end; end; # makes an approximate size Boltzmann sampler from a list of simple permutations genPermApproxSize:=proc(listeSimples,targetSize,epsilon,nbSamples) #option trace; local LS,efunc,paramBoltz,gen,minSize,maxSize,i,tab; LS:=makeLpoly(listeSimples); efunc:=gen_oracle_esp(LS): paramBoltz:=choose_param(efunc,targetSize,Digits,0,1): gen:=genPerm(listeSimples,Digits): minSize:=targetSize-targetSize*epsilon: maxSize:=targetSize+targetSize*epsilon: for i from 1 to nbSamples do tab[i]:=FAIL; while(tab[i]=FAIL) do try tab[i]:=gen(paramBoltz,maxSize); if _sizemaxSize then error("trop grand") end if; print(i); catch "trop grand": tab[i]:=FAIL; catch "trop petit": tab[i]:=FAIL; end try; od; od: tab end; tree_size:=proc(T) option remember; if T=Z then return 1 else convert(map(tree_size,[op(T)]),`+`) end if; end; tree2perm:=proc(T,L) local aux; aux:=proc(t,d) #option trace; local ts,i,cpt,s,i_size,simple,simple_L,tab,tab_size,n,r; ts:=tree_size(T); if t=Z then return d; else cpt:=d; s:=NULL; if op(0,t)=p then for i from 1 to nops(t) do i_size:=tree_size(op(i,t)); s:=(s,aux(op(i,t),cpt)); cpt:=cpt+i_size; od elif op(0,t)=m then for i from nops(t) by -1 to 1 do i_size:=tree_size(op(i,t)); s:=(aux(op(i,t),cpt),s); cpt:=cpt+i_size od; else simple:=op(op(0,t)); i:=1; for n in simple do tab_size[n]:=tree_size(op(i,t)); i:=i+1 od; tab_size[0]:=0; tab[0]:=d; for i from 1 to nops(t) do tab[i]:=tab[i-1]+tab_size[i-1]; od; i:=1; for n in simple do s:=(s,aux(op(i,t),op(tab[n]))); i:=i+1 od fi; s fi; end; [aux(T,1)] end; tab_cycles:=proc(tab) #option trace; local i,j,tlg; tlg[1]:=0; for i from 1 to nops(tab) do if (tlg[i]<>-1) then j:=tab[i]; tlg[i]:=1; while( j<>i ) do tlg[i]:=tlg[i]+1; tlg[j]:=-1; j:=tab[j]; od; end if; od; [seq(`if`(tlg[i]=-1,NULL,tlg[i]),i=1..nops(tab))]; end: