######################### # /home/yquem/algo/zimmermann/try/gaia.mpl ### # Convert from 5.2 to Release 3 # ## Title: Gaia ## Created: Tue Nov 17 15:33:26 1992 ## Author: Zimmermann Paul ## ## ## Description: draws a combinatorial object at random, uniformly macro( compile = gaia[compile], draw = gaia[draw], makedraw = `gaia/makedraw`, makeSubst = `gaia/makedraw/Subst`, makeProduct = `gaia/makedraw/Product`, doSubst = `gaia/doSubst`, doSubst_aux = `gaia/doSubst_aux` ): # using macro names rather than with( , ) macro( phi = numtheory[phi] ): # list of placeholders in function templates - something will # always be subs into these names macro( A_NONTERM=`gaia/placeholder/nontermA`, B_NONTERM=`gaia/placeholder/nontermB`, C_NONTERM=`gaia/placeholder/nontermC`, START=`gaia/placeholder/start`, INCR=`gaia/placeholder/incr`, VALC=`gaia/placeholder/valc`, VALB=`gaia/placeholder/valb`, VALA=`gaia/placeholder/vala`, PHIFUN=`gaia/placeholder/phifun`, NONTERM=`gaia/placeholder/nonterm`, INT_K=`gaia/placeholder/k`, UOFK=`gaia/placeholder/uofk`, SUBSFLAG=`gaia/placeholder/subsflag`, LABELFLAG=`gaia/placeholder/labelflag`, COUNT=`gaia/placeholder/count`, COUNT2=`gaia/placeholder/count2`, DRAW=`gaia/placeholder/draw`, DRAW2=`gaia/placeholder/draw2`, TERMLIST=`gaia/placeholder/termlist`, POLY=`gaia/placeholder/poly`, INDS=`gaia/placeholder/inds`, INT_K0=`gaia/placeholder/k0`, INT_K1=`gaia/placeholder/k1`, INT_V0=`gaia/placeholder/v0`, INT_V1=`gaia/placeholder/v1`, BODY=`gaia/placeholder/body` ): gaia['method']:='boustrophedon': # using Int and Theta # sys : { list of type declarations } # type declarations == name = type # type == Epsilon | Atom | Union(type list) | Product(type list) # | MultiConstructor(type compile := proc(Sys,typ) global `gaia/sys`, Draw; # typ = labelled | unlabelled # generates a set of drawing procedures local sys,eq,A,Phi,p,rhs,val,i,conflict,zeroproducers; option remember; # make sure the non-terminals don't clash with template variables conflict := {seq(gArg_.i,i=1..3), seq(gLoc_.i,i=1..9)} intersect indets(Sys,name); if conflict <> {} then ERROR(`the names`, conflict, `are reserved for the gaia system`); fi; # clear out remember tables needed for the valuations forget(`type/atom`); forget(`type/epsilon`); # clean out the remember table of user-given non-terminal names forget(`gaia/user_spec`); `gaia/sys`:=Sys,typ; if not type(Sys,set(name=anything)) or nops(Sys)=0 then ERROR(`invalid specification`,Sys) fi; if not member(typ,['labelled','unlabelled']) then ERROR(`invalid labelled type`,typ) fi; # syntax error checking plus some illegal grammars zeroproducers := has_zero_values(Sys); sys := gaia[standardform](Sys,typ); # check is specification is well founded userinfo(1,gaia,`Checking if the specification is well founded ...`); val := `gaia/is_well_founded`(Sys,typ,zeroproducers); if type(val,InfiniteValuation(set)) then ERROR(op(op(val) intersect indets(Sys,name)), `do not derive any object of finite size`) elif val = false then ERROR(`specification is not well founded`) fi; # generates counting procedures `gaia/count_CNF`(Sys,typ, val); # generates drawing procedures for A in indets(sys,epsilon) do Draw[args][A]:=makedraw(A,val) od; for eq in sys do A := op(1,eq); rhs := op(2,eq); # due to a small glitch caused by using a remember table in # the type functions for atom and epsilon, the test for # non-terminal should be done first # but not to worry - this only affects productions of the # form a=b=Atom/Epsilon if type(rhs, nonterminal) then p:=subs(A_NONTERM=rhs, proc(gArg_1,gArg_2) Draw[A_NONTERM](gArg_1,gArg_2) end) elif type(rhs,atom) then p:=subs(A_NONTERM=A,proc() A_NONTERM end) elif type(rhs,epsilon) then # ack! Sometimes the system introduces a production of the form # T1=Epsilon, which should print as Epsilon, not as T1 # better solution would be to stop the system from creating that # production in the first place? if member(eq,Sys) then # it's the user's idea p:=subs(A_NONTERM=A,proc() A_NONTERM end); else p:=makedraw(rhs,val); fi elif type(rhs,function) then p := makedraw(rhs, val); Phi := op(0,rhs); p := subs(`gaia/print`[Phi]=cat(`gaia/print/`,Phi),op(p)) else ERROR(`unexpected production`,rhs) fi; Draw[args][A] := subs(Draw=Draw[args],Count=Count[args], Count2=Count2[args],Draw2=Draw2[args],op(p)) od; forget(`type/atom`); userinfo(1,gaia,`number of non terminals=`,nops(sys)); end: draw := proc(sys,typ,y:name,n:nonnegint) #draws an object of type y and size n # where y is the first element defined in sys local yn; if not member(y,indets(sys,name)) then ERROR(`undefined non-terminal`,y) fi; compile(sys,typ); yn:=Count[sys,typ][y](n); if yn=0 then ERROR(`there is no structure of this size`) fi; if typ='labelled' then label(Draw[sys,typ][y](n,Rand(yn)()),n) else Draw[sys,typ][y](n,Rand(yn)()) fi end: # draw Rand := proc(n) option remember; rand(n) end: makedraw := proc(rhs, val) local Phi; if type(rhs,epsilon) then op(`gaia/print/Epsilon`); elif type(rhs,function) then Phi := op(0,rhs); if Phi=Union then `makedraw/Union`(val,op(rhs)) elif type(Phi,product) then makeProduct(op(1,rhs),op(2,rhs),Phi, val) elif type(Phi,subst) then makeSubst(op(rhs),Phi) elif type(Phi,delta) then `makedraw/Delta`(op(rhs),op(Phi)) elif Phi=Theta then `makedraw/Theta`(op(rhs)) elif Phi=Int then `makedraw/Int`(op(rhs)) else ERROR(`unknown constructor`,Phi) fi else ERROR(`bad production`,rhs) fi end: `makedraw/Union` := proc(val,A,B) if nargs=3 then if type(A,epsilon) and val[B]>0 then subs(A_NONTERM=A,B_NONTERM=B, proc(gArg_1,gArg_2) if gArg_1=0 then Draw[A_NONTERM](gArg_1,gArg_2) else Draw[B_NONTERM](gArg_1,gArg_2) fi end) else # general case subs(A_NONTERM=A,B_NONTERM=B, proc(gArg_1,gArg_2) if gArg_2 if op(0,gArg_1)=Cycle then op(gArg_1) else gArg_1 fi, gLoc_1); `gaia/print/Cycle`(op(gLoc_1)); end)) elif type(typ,integer) then RETURN(subs(A_NONTERM=A,INT_K=typ, proc(gArg_1,gArg_2) [Draw[A_NONTERM](gArg_1/INT_K,gArg_2) $ INT_K] end)) else ERROR(`unknown Delta type`,typ) fi; subs(UOFK=uofk,A_NONTERM=A,START=start,PHIFUN=p, proc(gArg_1,gArg_2) local gLoc_1,gLoc_2,gLoc_3,gLoc_4; gLoc_2 := gArg_2; for gLoc_1 from START to gArg_1 do if gArg_1 mod gLoc_1 = 0 then gLoc_3 := UOFK * Count[A_NONTERM](gArg_1/gLoc_1); if gLoc_2 < gLoc_3 then break else gLoc_2 := gLoc_2-gLoc_3 fi; fi od; if UOFK=1 then gLoc_4 := gLoc_2 else gLoc_4 := iquo(gLoc_2,UOFK) fi; PHIFUN(Draw[A_NONTERM](gArg_1/gLoc_1,gLoc_4),gLoc_1) end) end: `gaia/power/template` := proc(gArg_1,gArg_2,gArg_3) local gLoc_1,gLoc_2,gLoc_3,gLoc_4,gLoc_5,gLoc_6; option remember; if gArg_2=1 then DRAW[A_NONTERM](gArg_1,gArg_3) else gLoc_1:=gArg_3; if LABELFLAG then gLoc_6 := gLoc_5-1 else gLoc_6 := gArg_1-1 fi; for gLoc_5 to gArg_1 - gArg_2 + 1 do gLoc_4 := binomial(gArg_1-1, eval(gLoc_6)); gLoc_3 := gLoc_4 * COUNT[A_NONTERM](gLoc_5)* COUNT2[A_NONTERM](gArg_1-gLoc_5,gArg_2-1); if gLoc_1 < gLoc_3 then break else gLoc_1 := gLoc_1 - gLoc_3 fi od; gLoc_2 := iquo(iquo(gLoc_1,gLoc_4), COUNT[A_NONTERM](gLoc_5),'gLoc_1'); DRAW[A_NONTERM](gLoc_5, gLoc_1), DRAW2[A_NONTERM](gArg_1 - gLoc_5, gArg_2 - 1, gLoc_2) fi end: `gaia/Subst/template` := proc(gArg_1,gArg_2) # B(Z <- A) local gLoc_1,gLoc_2,gLoc_3,gLoc_4,gLoc_5; option remember; gLoc_2 := gArg_2; if SUBSFLAG then gLoc_4 := gArg_1-1 else gLoc_4 := gArg_1 fi; for gLoc_1 to gArg_1 + gLoc_4 do # 2 gLoc_5 := Count[B_NONTERM](gLoc_1) * Count2[A_NONTERM](gArg_1,gLoc_1); if gLoc_2 < gLoc_5 then break else gLoc_2 := gLoc_2-gLoc_5 fi od; gLoc_3 := iquo(gLoc_2,Count[B_NONTERM](gLoc_1),'gLoc_2'); doSubst([Draw2[A_NONTERM](gArg_1,gLoc_1,gLoc_3)], Draw[B_NONTERM](gLoc_1,gLoc_2)) end: makeSubst := proc(A,B,Phi) global Draw2; local flag1,flag2; flag1 := evalb(Phi=`Subst*`); flag2 := evalb(`gaia/sys`[2]='labelled'); Draw2[`gaia/sys`][A] := subs(A_NONTERM=A,LABELFLAG=flag2, COUNT=Count[`gaia/sys`], COUNT2=Count2[`gaia/sys`],DRAW=Draw[`gaia/sys`], DRAW2=Draw2[`gaia/sys`],op(`gaia/power/template`)); subs(SUBSFLAG=flag1,A_NONTERM=A,B_NONTERM=B,op(`gaia/Subst/template`)); end: `gaia/heuristic` := proc(sys,typ,A,B,n0) local i,sf,sba,sbo,n,st,b; # best way to generate Prod(A,B) between forwards, backwards, boustrophedon if typ='labelled' then b:=binomial(n,i) else b:=1 fi; sf:=0; sba:=0; sbo:=0; for n from n0 to 2*n0 while sf=0 do st := sum(b*'Count[sys,typ][A](i)*Count[sys,typ][B](n-i)',i=0..n); if st=0 then next fi; sf := sum(b*'(i+1)*Count[sys,typ][A](i)*Count[sys,typ][B](n-i)',i=0..n); sba := sum(b*'(i+1)*Count[sys,typ][B](i)*Count[sys,typ][A](n-i)',i=0..n); sbo := sum(b*'2*(i+1)*Count[sys,typ][A](i)*Count[sys,typ][B](n-i)',i=0..iquo(n+1,2)) +sum(b*'2*(i+1)*Count[sys,typ][B](i)*Count[sys,typ][A](n-i)',i=0..iquo(n,2)); od; userinfo(3,gaia,`forwards:`,evalf(sf/st),`backwards:`,evalf(sba/st), `boustrophedon:`,evalf(sbo/st)); if sf<=sba then if sf=` := proc(x) if type(x,integer=2 then assign(args[2],op(1,x)+1) fi; true elif type(x,integer<=identical(card)) then if nargs>=2 then assign(args[2],op(1,x)) fi; true else false fi end: `type/gaia/card=` := proc(x) if type(x,identical(card)=integer) then if nargs>=2 then assign(args[2],op(2,x)) fi; true else false fi end: `type/gaia/card<=` := proc(x) if type(x,identical(card)=2 then assign(args[2],op(2,x)-1) fi; true elif type(x,identical(card)<=anything) then if nargs>=2 then assign(args[2],op(2,x)) fi; true else false fi end: alias( count = gaia[count], count_CNF = `gaia/count_CNF`, makecount = `gaia/makecount`, mkproc = `gaia/mkproc`, mkproc0 = `gaia/mkproc0`, mkproc2 = `gaia/mkproc2` ): gaia[usegfun]:=false: count := proc(sys,typ,y,n) if not type(y,name) then ERROR(`invalid non terminal`,y) fi; if not type(n,nonnegint) then ERROR(`invalid size`,n) fi; compile(sys,typ); Count[sys,typ][y](n) end: # algorithme Calcul des Coefficients page 74 count_CNF:=proc (Sys,typ, val) global Count; local A,eq,sys,hpart,res; res := NULL; sys := gaia[standardform](Sys,typ); userinfo(1,gaia,`Generating counting procedures ...`); for A in indets(sys,epsilon) do Count[Sys,typ][A]:=makecount(Epsilon, val) od; for eq in sys do A := op(1,eq); Count[Sys,typ][A] := subs(Count=Count[Sys,typ],Count2=Count2[Sys,typ], makecount(op(2,eq), val)) od; if gaia['usegfun']=true and gfun_installed()=true then hpart := algsys_to_eqs(hol_to_algsys(holonomic_part(sys))); userinfo(2,gaia,`Holonomic part is:`,hpart); for eq in hpart do A := op(indets(eq) minus {'z'}); Count[Sys,typ][A]:=`makecount/holonomic`(eq) od fi; res end: # count gfun_installed := proc() global invborel; local l,res; option remember; res:=traperror(with(gfun)); if res=lasterror then # let's try to install it ourselves res:=traperror(with(share)); if res=lasterror then lprint(`Cannot find the share library, so cannot use gfun`); RETURN(false) else res:=traperror(readshare(gfun,calculus)); if res=lasterror then lprint(`The gfun package is not properly installed`); RETURN(false) fi; fi; fi; if type(Heaviside,procedure) then # Release 2 invborel:=exptoord fi; l := [algfuntodiffeq,diffeqtorec,rectoproc,invborel]; if has(map(type,l,procedure),false) then lprint(`The gfun version is not the good one`); RETURN(false) fi; true end: makecount := proc(rhs,val) local Phi; if type(rhs,nonterminal) then `makecount/name`(rhs) elif type(rhs,epsilon) then `makecount/epsilon`() elif type(rhs,atom) then `makecount/atom`() elif type(rhs,function) then Phi := op(0,rhs); if Phi=Union then `makecount/Union`([op(rhs)], val) elif type(Phi,product) then `makecount/Product`[`gaia/sys`[2]](op(1,rhs),op(2,rhs), val) elif type(Phi,subst) then `makecount/Subst`(op(rhs),Phi, val) elif type(Phi,delta) then `makecount/Delta`(op(rhs),op(Phi)) elif Phi=Int then `makecount/Int`(op(rhs)) elif Phi=Theta then `makecount/Theta`(op(rhs)) else ERROR(`unknown constructor`,Phi) fi else ERROR(`bad production`,rhs) fi; end: `makecount/epsilon` := proc() proc(gArg_1) option remember; if gArg_1=0 then 1 else 0 fi end end: `makecount/atom` := proc() proc(gArg_1) if gArg_1=1 then 1 else 0 fi end end: `makecount/name`:=proc(A) subs(A_NONTERM=A, proc(gArg_1) option remember; Count[A_NONTERM](gArg_1) end) end: `makecount/Int`:=proc(A) subs(A_NONTERM=A, proc(gArg_1) option remember; if gArg_1=0 then 0 else Count[A_NONTERM](gArg_1)/gArg_1 fi end) end: `makecount/Theta`:=proc(A) subs(A_NONTERM=A, proc(gArg_1) option remember; gArg_1*Count[A_NONTERM](gArg_1) end) end: `makecount/Delta` := proc(A,typ) local uofk,start,k; # warning: the 'gLoc_1' in this procedure is the same gLoc_1 that is local # in the unnamed procedure being 'subs' into, hence they must have # the same name - this ought to be changed eventually start:=1; if typ=Set then uofk:=1 elif typ=Cycle then uofk:='phi(gLoc_1)' elif type(typ,Set(integer)) then uofk:=1; start:=op(typ) elif type(typ,Cycle(integer)) then uofk:='phi(gLoc_1)'; start:=op(typ) elif type(typ,integer) or type(typ, InCycle(integer)) then if type(typ,integer) then k:= typ else k := op(typ) fi; RETURN(subs(A_NONTERM=A,INT_K=k, proc(gArg_1) if gArg_1 mod INT_K=0 then Count[A_NONTERM](gArg_1/INT_K) else 0 fi end)) else ERROR(`unknown Delta type`,typ) fi; subs(UOFK=uofk,A_NONTERM=A,START=start, proc(gArg_1) local gLoc_1,gLoc_2; option remember; gLoc_2 := 0; for gLoc_1 from START to gArg_1 do if gArg_1 mod gLoc_1 = 0 then gLoc_2:=gLoc_2+UOFK*Count[A_NONTERM](gArg_1/gLoc_1) fi od; gLoc_2 end) end: `count/Product/template/labelled` := proc(gArg_1) local gVar_1,gVar_2,gVar_3; option remember; gVar_2 := 0; gVar_3 := START; # gVar_3 = binomial(gArg_1,gVar_1) for gVar_1 from VALB to gArg_1-VALC do gVar_2 := gVar_2 + gVar_3*Count[B_NONTERM](gVar_1)*Count[C_NONTERM](gArg_1-gVar_1); gVar_3 := gVar_3 * (gArg_1-gVar_1) / (gVar_1+1) od; gVar_2 end: `count/Product/template/unlabelled` := proc(gArg_1) local gVar_1,gVar_2; option remember; gVar_2 := 0; for gVar_1 from VALB to gArg_1-VALC do gVar_2 := gVar_2 + Count[B_NONTERM](gVar_1)*Count[C_NONTERM](gArg_1-gVar_1); od; gVar_2 end: # the point of the k=gArg_1 substitution is that the default # function body is written in terms of k, but we must use # gArg_1 to avoid possible conflicts with non-terminal names # the assumption is that none of the other arguments are in terms of k mkproc0 := proc(k,default) local body; body := subs(k=gArg_1, default); subs(BODY=body,proc(gArg_1) option remember; BODY end); end: mkproc := proc(k0,v0,k,default) local body; body := subs(k=gArg_1, default); subs(BODY=body,INT_K0=k0,INT_V0=v0, proc(gArg_1) option remember; if gArg_1=INT_K0 then INT_V0 else BODY fi end); end: mkproc2 := proc(k0,v0,k1,v1,k,default) local body; body := subs(k=gArg_1, default); subs(BODY=body,INT_K0=k0,INT_V0=v0,INT_K1=k1,INT_V1=v1, proc(gArg_1) option remember; if gArg_1=INT_K0 then INT_V0 elif gArg_1=INT_K1 then INT_V1 else BODY fi end); end: `makecount/Union` := proc(l,val) local n0,def,n,v0,v1,x,eps,r; eps:=indets(l,epsilon); if has(l,Atom) then n0:=1 elif eps<>{} then n0:=0 else n0:=-1 fi; if n0>=0 then #subs(Atom=,l); r:=select(proc(x,val) evalb(val[x]=0) end,l,val); v0:=subs(seq(Count[x](0)=1,x=eps),convert([seq(Count[x](0),x=r)],`+`)); fi; if n0>=1 then r:=subs(seq(x=NULL,x=eps),Atom=,l); v1:=subs(Count[Atom](1)=1,convert([seq(Count[x](1),x=r)],`+`)); fi; def:=convert([seq(Count[x](n), x=eval(subs(seq(x=NULL,x=eps),Atom=,l)))],`+`); if n0=-1 then mkproc0(n,def) elif n0=0 then mkproc(0,v0,n,def) else mkproc2(0,v0,1,v1,n,def) fi end: `count/power/template` := proc(gArg_1,gArg_2) # labelled : nb of objects of Set(A,card=gArg_2) of size gArg_1 # unlabelled : nb of objects of Sequence(A,card=gArg_2) of size gArg_1 local gLoc_1,gLoc_2,gLoc_3,gLoc_4; option remember; if gArg_2=1 then COUNT[A_NONTERM](gArg_1) else gLoc_1:=0; if LABELFLAG then gLoc_3 := gLoc_2-1 else gLoc_3 := gArg_1-1 fi; # necessarily val(A)>=1 # Set(A,card=gArg_2) = A_min Set(A,card=gArg_2-1) for gLoc_2 to gArg_1 - gArg_2 + 1 do gLoc_4 := binomial(gArg_1-1,eval(gLoc_3)); gLoc_1 := gLoc_1 + gLoc_4 * COUNT[A_NONTERM](gLoc_2) * COUNT2[A_NONTERM](gArg_1-gLoc_2,gArg_2-1) od; gLoc_1 fi end: `count/Subst/template` := proc(gArg_1) local gLoc_1,gLoc_2,gLoc_3; option remember; gLoc_2:=0; if SUBSFLAG then gLoc_3:= gArg_1-1 else gLoc_3 := gArg_1 fi; if VALA>1 then gLoc_3 := min(gLoc_3, iquo(gArg_1, VALA)) fi; for gLoc_1 to gLoc_3 do # B of size gLoc_1 [ A^gLoc_1 of size gArg_1 ] gLoc_2 := gLoc_2 + Count[B_NONTERM](gLoc_1) * Count2[A_NONTERM](gArg_1,gLoc_1) od; gLoc_2 end: `makecount/Subst` := proc(A,B,Phi, val) # B where Z is replaced by A global Count2; local flag1,flag2; flag1 := evalb(Phi=`Subst*`); flag2 := evalb(`gaia/sys`[2]=labelled); Count2[`gaia/sys`][A] := subs(A_NONTERM=A,LABELFLAG=flag2, COUNT=Count[`gaia/sys`], COUNT2=Count2[`gaia/sys`], op(`count/power/template`)); subs(SUBSFLAG=flag1,A_NONTERM=A,B_NONTERM=B,VALA=val[A], op(`count/Subst/template`)) end: `makecount/Product`[labelled] := proc(B,C,val) local n,p; if type(B,atom) then if type(C,atom) then mkproc(2,2,n,0) else mkproc(0,0,n,n*Count[C](n-1)) fi elif type(C,atom) then mkproc(0,0,n,n*Count[B](n-1)) else p := op(`count/Product/template/labelled`); subs(B_NONTERM=B,C_NONTERM=C,VALB=val[B],VALC=val[C], START=binomial(op(1,op(p)),val[B]),op(p)) fi; end: `makecount/Product`[unlabelled] := proc(B,C,val) local n,p; if type(B,atom) then if type(C,atom) then mkproc(2,1,n,0) else mkproc(0,0,n,Count[C](n-1)) fi elif type(C,atom) then mkproc(0,0,n,Count[B](n-1)) else p := op(`count/Product/template/unlabelled`); subs(B_NONTERM=B,C_NONTERM=C,VALB=val[B],VALC=val[C],op(p)) fi; end: alias( MultiConstructors = `gaia/MultiConstructors`, N_Phi = `gaia/N_Phi`, properties = `gaia/properties` ): PutsInCNF := proc(expr) global `gaia/def`; local r; # returns a non terminal name for expr if type(expr,epsilon) then if nargs>=2 then `type/epsilon`(args[2]):=true; `gaia/def`(args[2]):=expr; RETURN(args[2]) else RETURN(expr) fi elif expr=Atom then if nargs=2 then `type/atom`(args[2]):=true; RETURN(args[2]) else RETURN(expr) fi elif type(expr,name) then # RETURN(expr) fixed PZ 150394 if nargs=2 and `gaia/user_spec`(args[2]) then r:=eval(expr) else RETURN(expr) fi elif type(expr,function) then if type(op(0,expr),std_constr) then r:=map(procname,expr) else # substitution (this syntax made illegal by EM, April 94) ERROR(`invalid construction`, expr); fi else ERROR(`invalid construction`,expr) fi; NameOf(r,args[2..nargs]) end: NameOf := proc(expr) # returns the name for expr in CNF, or a new name local n,r; n := op(4,op(procname))[expr]; if type(n, string) then # if there is a user-given name as well as a system name for # the same expr, keep the user name as well # if nargs>=2 and member(args[2], indets(`gaia/sys`[1])) # and not(evalb(n=args[2])) then if nargs>=2 and `gaia/user_spec`(args[2]) and not(evalb(n=args[2])) then procname(n) := args[2]; # redirect new name to old one else n; fi elif nargs>=2 then procname(expr) := args[2] else r:=newname(); userinfo(2,gaia,`Introducing new type`,r,`for`,expr); procname(expr) := r fi; end: `gaia/nb` := -1: newname := proc() global `gaia/nb`; local s,r; `gaia/nb`:=`gaia/nb`+1; r := {unames()}: s := T.`gaia/nb`; while s<>eval(s) or member(s,r) do `gaia/nb`:=`gaia/nb`+1; r := {unames()}: s := T.`gaia/nb` od; s end: Greene[Set,unlabelled] := proc(B,A) local j,init,k,i,l; if nargs=2 then B=Union(Epsilon,Int(Prod[Set](Delta[Set](Theta(A)),B))) else if type(args[3],`gaia/card>=`('k')) then if k=0 then RETURN(procname(B,A)) fi; l := array(0..k); l[0]:=newname(); init := procname(l[0],A); # sets with no restriction : B0=Set(A) for i to k do if i=k then l[i] := B else l[i] :=newname() fi; init := init,eval(l[i])=Int(Union(seq(Prod[Set](Delta[j](Theta(A)), l[i-j]),j=1..i-1),Prod[Set](Delta[Set(i)](Theta(A)),l[0]))) od elif type(args[3],`gaia/card=`(k)) then if k=0 then RETURN(B=Epsilon) elif k=1 then RETURN(B=Prod[Set](A,Epsilon)) fi; l := array(0..k); l[0]:=Epsilon; l[1]:=Prod[Set](A,Epsilon); init := NULL; for i from 2 to k do if i=k then l[i]:=B else l[i]:=newname() fi; init := init,l[i]=Int(Union(seq(Prod[Set](Delta[j]( Theta(A)),l[i-j]),j=1..i))) od elif type(args[3],`gaia/card<=`(k)) then if k=0 then RETURN(B=Epsilon) elif k=1 then RETURN(B=Union(Epsilon,A)) fi; l := array(0..k); l[0]:=Epsilon; l[1]:=Union(Epsilon,A); init := NULL; for i to k do if i=k then l[i]:=B else l[i]:=newname() fi; init := init,l[i]=Union(Int(Union(seq(Prod[Set]( Delta[j](Theta(A)),l[i-j]),j=1..i))),Epsilon) od else ERROR(`bad restriction on cardinality`) fi; init fi; end: Greene[Set,labelled] := proc(A,B) # A = Union(1,Int(Prod[Set](Theta(B),A))) local j,p,x,init,A0,toto; if nargs=2 then A=Union(Epsilon,Int(Prod[Set](Theta(B),A))) else if type(args[3],`gaia/card>=`(j)) then A0:=newname(); init := procname(A0,B); userinfo(2,gaia,`Introducing new type`,A0,'verifying',init); p := Int(Prod[Set](Theta(B),x)); elif type(args[3],`gaia/card=`(j)) then A0:=Epsilon; init := NULL; p := Int(Prod[Set](Theta(B),x)); elif type(args[3],`gaia/card<=`(j)) then A0:=Epsilon; init := NULL; p := Union(Epsilon,Int(Prod[Set](Theta(B),x))); else ERROR(`bad restriction on cardinality`) fi; if not type(j,integer) then ERROR(`should not happen`) fi; toto:=x; to j do toto:=subs(x=p,toto) od; A=subs(x=A0,toto),init fi; end: Greene[Sequence,unlabelled]:=proc(A,B) # Sequence(B) ==> Union(1,Prod(B,A)) local j,p,x,init,A0,toto; if nargs=2 then NameOf(Sequence(B)):=A; A=Union(Epsilon,Prod[Seq](B,A)) else if type(args[3],`gaia/card>=`('j')) then A0 := NameOf(Sequence(B)); if not type(A0,name) then A0:=newname() fi; init := procname(A0,B); p := Prod[Seq](B,x) elif type(args[3],`gaia/card=`(j)) then A0:=Epsilon; init := NULL; p := Prod[Seq](B,x); elif type(args[3],`gaia/card<=`(j)) then A0:=Epsilon; init := NULL; p := Union(Epsilon,Prod[Seq](B,x)); else ERROR(`bad restriction on cardinality`) fi; if not type(j,integer) then ERROR(`should not happen`) fi; toto:=x; to j do toto:=subs(x=p,toto) od; A=subs(x=A0,toto),init fi; end: Greene[Sequence,labelled]:=op(Greene)[Sequence,unlabelled]: Greene[Cycle,unlabelled] := proc(B,A) local k, j, r; if nargs=2 then r:=Delta[Cycle](Prod[Cyc](Theta(A),Sequence(A))) elif type(args[3],`gaia/card>=`(k)) then r:=Union(seq(Delta[InCycle(j)](Prod[Cyc](Theta(A), Sequence(A,card>=ceil(k/j-1))))$phi(j),j=1..k-1), Delta[Cycle(k)](Prod[Cyc](Theta(A),Sequence(A)))) elif type(args[3],`gaia/card=`(k)) then # ack! small kludge to fix printing bug if k=1 then r:=Prod[Cyc](Theta(A),Epsilon) else r:=Union(seq(Delta[InCycle(j)](Prod[Cyc](Theta(A), A$(k/j-1)))$phi(j),j=numtheory[divisors](k))) fi; elif type(args[3],`gaia/card<=`(k)) then r:=Union(seq(Delta[InCycle(j)](Prod[Cyc](Theta(A), Sequence(A,card<=iquo(k,j)-1)))$phi(j),j=1..k)) else ERROR(`invalid cardinality restriction`,args[3]) fi; B=Int(r) end: Greene[Cycle,labelled] := proc(A,B) # Cycle(B) ==> OProduct(B,Sequence(B)) local j,r; if nargs=2 then r:=Sequence(B) elif type(args[3],`gaia/card>=`(j)) then r:=Sequence(B,card>=j-1) elif type(args[3],`gaia/card=`(j)) then r:=Sequence(B,card=j-1) elif type(args[3],`gaia/card<=`(j)) then r:=Sequence(B,card<=j-1) else ERROR(`invalid cardinality restriction`,args[3]) fi; A=Int(Prod[Cyc](Theta(B),eval(r))) end: readlib(forget): # write a real function for this, that is capable of returning false `gaia/user_spec` := proc(A) local thisismyuser_specfunction; false end: GreeneTransform := proc(Sys,typ) local sys,x,rhs,t,n,wanted_types,constr,r; sys := Sys; forget(NameOf); for x in sys do Nameof(op(2,x)):=op(1,x); `gaia/user_spec`(op(1,x)):=true od; for constr in [Set,Cycle,Sequence] do wanted_types:=specfunc(anything,constr); r:=select(hastype,sys,wanted_types); while r<>{} do x := op(1,r); rhs := op(2,x); if type(rhs,wanted_types) then sys := subs(x=Greene[op(0,rhs),typ](op(1,x),op(rhs)),sys) else t:=op(1,indets(rhs,wanted_types)); n := NameOf(t); if n='NameOf(t)' then n:=newname() fi; sys := subs(x=subs(t=n,x),sys) union {Greene[op(0,t),typ](n,op(t))} fi; r:=select(hastype,sys,wanted_types) od; od; eval(subs(Union=proc(x) if nargs=1 then x else 'Union(args)' fi end, Set(1)=Set,Cycle(1)=Cycle,Delta[1]=,Delta[InCycle(1)]=, Theta=proc(x) if type(x,atom) then x else 'Theta(x)' fi end,sys)) end: `gaia/CNF/Prod` := proc() local p; if type(procname,indexed) then if op(1,procname)='Cyc' then p:=Prod['InCyc'] else p:=procname fi; else p:=Prod['Prd'] fi; if nargs=0 then Epsilon elif nargs=1 then args elif nargs=2 then 'procname'(args) else 'procname'(args[1],p(args[2..nargs])) fi end: gaia[standardform] := proc(Sys,typ) global `gaia/nb`, Prod; # puts in Chomsky Normal Form local x,sys,l,y,z,r; option remember; if type(Sys,set) then sys:=Sys else ERROR(`not a valid specification`,Sys) fi; if has(sys,Z) then sys:=sys union {Z=Atom} fi; `gaia/nb` := -1; sys := GreeneTransform(sys,typ); # transform n-ary into binary products Prod := op(`gaia/CNF/Prod`); sys := eval(sys); Prod := 'Prod'; # then define aliases such that each production is X -> Phi(Y) forget(NameOf); l:=NULL; for x in sys do r:=PutsInCNF(op(2,x),op(1,x)); if type(r,name) and r<>op(1,x) then l:=l,op(1,x)=r fi; od; sys:=op(4,op(NameOf)); if sys=NULL then # ERROR(`invalid specification`,Sys) fi; sys:={{{}}} fi; r:=eval(subs(l,op(op(sys)))); sys:={seq(op(2,x)=op(1,x),x=r)}; r := op(4,eval(`type/atom`)); if r<>NULL then sys := sys union {seq(op(x)=Atom,x=[indices(r)])} fi; r:=(op(4,eval(`type/epsilon`))); if r<>NULL then sys:=sys union {seq(op(x)=`gaia/def`(op(x)),x=[indices(r)])} fi; if sys={} then ERROR(`invalid specification`,Sys) fi; # look for x=Theta(y), y=Int(z) for y in map( proc(y,s) if hastype(s,name=Theta(identical(y))) and hastype(s,identical(y)=Int(name)) then y fi end, indets(sys,name),sys) do x := NameOf(Theta(y)); # ack! don't remove user-given non-terminals if not `gaia/user_spec`(x) and not `gaia/user_spec`(y) then z := op(op(2,op(1,select(type,sys,identical(y)=Int(name))))); sys := subs(z=x,sys minus {x=Theta(y),y=Int(z)}); if has(sys,y) then sys:=sys union {y=Int(x)} fi; fi; od; sys end: # Sys is the grammar in CNF, zprods is the boolean table of user-given # non-terminals along with whether or not they have valuation 0 # for internal use : argument is already in CNF valuation_CNF:=proc (Sys,zprods) # algorithm A page 28 of [Zimmermann91] local eq,x,v,changed,Phi,y,sys,T,k,r,inds,m; option remember; sys := Sys; # initialization v:=table([seq(x=0,x=indets(sys,epsilon)), seq(x=1,x=indets(sys,atom))]); userinfo(2,gaia,`Initializing valuations ...`); for eq in sys do if type(op(2,eq),epsilon) then v[op(1,eq)]:=0; sys:=sys minus {eq} elif type(op(2,eq),atom) then v[op(1,eq)]:=1; sys := sys minus {eq} elif zprods[op(1,eq)]=true then v[op(1,eq)]:=0 # we already know this one else v[op(1,eq)]:=infinity # non terminals fi od; changed:=true; while changed=true do changed:=false; for eq in sys do T := op(1,eq); x := op(2,eq); if type(x,name) then r:=v[x] elif type(x,function) then Phi := op(0,x); if Phi=Epsilon then r:=0 elif Phi=Union then r:=min(seq(v[y],y=x)) elif type(Phi,product) then r:=convert([seq(v[y],y=x)],`+`) elif Phi=Subst then r:=v[op(2,x)]*v[op(1,x)] elif Phi=`Subst*` then if v[op(1,x)]=1 then r:=v[op(2,x)]+1 else r:=v[op(2,x)]*v[op(1,x)] fi elif member(Phi,{Theta,Int}) then r:=v[op(x)] elif type(Phi,delta) then r:=v[op(x)] else ERROR(`unknown constructor`,Phi) fi; else ERROR(`bad production`) fi; x:=r; inds:=indets(x,v[name]); if inds<>{} then ERROR(`undefined non terminals:`,seq(op(x),x=inds)) fi; m:=readlib(evalr)(min(x,v[T])); if m<>v[T] then v[T]:=m; userinfo(3,gaia,`new valuation for type`,T,'is',v[T]); changed:=true fi; od; od; # fill in remember table for calls without the "zprods" parameter procname(Sys) := op(v); op(v) end: # valuation # according to definition 6 page 26 properties := table([ # f_Phi, 2a or 2b, [(P|Q)*] Sequence = [<0>,`2a`,[P]], Cycle = [,`2a`,[P]], Set = [<0>,`2a`,[P]], # multiset Int = [,`2b`,[Q,P]], Theta = [,`2b`,[Q,P]] ]): `gaia/MultiConstructors` := map(op,[indices(properties)]): N_Phi:=table(map(proc(x) x=nops(op(3,properties[x]))-1 end,MultiConstructors)): CoeffOneDelta := proc(Phi) if member(Phi,{Delta[Set],Delta[Cycle]}) then 1 elif type(Phi,Delta[integer]) or type(Phi,Delta[InCycle(integer)]) then if op(Phi)=1 then 1 else 0 fi elif type(Phi,{Delta[Set(integer)],Delta[Cycle(integer)]}) then if op(op(Phi))<=1 then 1 else 0 fi else ERROR(`unknown Delta constructor`,Phi) fi end: property := proc(Phi,n) # P or Q local l; l := op(3,properties[Phi]); op(min(n+1,nops(l)),l); end: has_zero_values := proc(sys) local zval,rhs,nonterm,eq,oldval,changed; # initialization for eq in sys do if type(op(2,eq),epsilon) then zval[op(1,eq)] := true else zval[op(1,eq)] := false; fi; od; changed := true; while changed=true do changed := false; for eq in sys do nonterm := op(1,eq); rhs := op(2,eq); oldval := zval[nonterm]; zval[nonterm] := has_zero(rhs,zval); if zval[nonterm] <> oldval then changed := true fi; od; od; eval(zval); end: # subfunction - given rhs of a production, determine if it can produce # elements of size 0 # we assume this function goes near the beginning, so it does lots of # error checking has_zero := proc(rhs,zval) local r,Phi,k; if type(rhs,epsilon) then r := true; elif type(rhs, atom) then r:= false; elif type(rhs,name) then if evalb('zval[rhs]'=zval[rhs]) then ERROR(`undefined non-terminal`, rhs) fi; r := zval[rhs] elif type(rhs, function) then Phi := op(0, rhs); if Phi=Union then r := convert(map(procname, [op(rhs)], zval), `or`) ; elif type(Phi,product) then r := convert(map(procname, [op(rhs)], zval), `and`) ; elif type(Phi, usecard) then if nops(rhs) > 2 then ERROR(`too many arguments`, rhs) fi; # Cycles cannot be of length 0 if Phi=Cycle then if nops(rhs) = 2 and ( type(op(2,rhs), `gaia/card<=`('k')) or type(op(2,rhs), `gaia/card=`('k')) ) and k=0 then ERROR(`there is no cycle of length 0`, rhs); fi; fi; r := procname(op(1,rhs),zval); if r=true then # the inside produces elements of size 0 # we reject all such cases - some because they derive # infinitely many objects,some because the code to handle # such bizarre cases isn't written yet if nops(rhs) = 1 or type(op(2,rhs), `gaia/card>=`) then ERROR(rhs, `derives infinitely many objects of size 0`); elif type(op(2,rhs), `gaia/card=`('k')) or type(op(2,rhs), `gaia/card<=`('k')) then ERROR(op(1,rhs), `produces an element of size 0 inside`, Phi, `which is not allowed`); else ERROR(`bad cardinality restriction`, rhs); fi; else # the inside production cannot produce elements of size 0 if nops(rhs) = 1 or type(op(2,rhs), `gaia/card<=`) then if Phi = Cycle then r := false; #there are no cycles of size 0 by convention else r := true; fi; elif type(op(2,rhs), `gaia/card>=`('k')) or type(op(2,rhs), `gaia/card=`('k')) then if Phi=Cycle then r := false else r:=evalb(k=0); fi; else ERROR(`bad cardinality restriction`, rhs); fi; fi; elif type(Phi, subst) then if nops(rhs) <> 2 then ERROR(`Subst takes two arguments`, rhs); fi; if procname(op(1,rhs),zval) = true then ERROR(`cannot substitute objects of size 0`, rhs); else r := procname(op(2,rhs),zval); fi; elif type(Phi, delta) then if nops(rhs) <> 1 then ERROR(`too many arguments`, rhs); fi; r := procname(op(rhs), zval); elif Phi=Theta or Phi=Int then if nops(rhs) <> 1 then ERROR(`too many arguments`, rhs); fi; r := procname(op(rhs), zval); if r = true then ERROR(Phi, `must not have arguments that produce elements of size 0`, rhs); fi; else ERROR(`unknown constructor`, Phi); fi; else ERROR(`invalid production`, rhs); fi; r; end: `gaia/is_well_founded` := proc(Sys,typ,zprods) # algorithm B page 33 # Sys is the original grammar # returns the valuation table if well-founded, false or InfiniteValuation # otherwise local t,Y,Phi,N,n,g,l,val,alpha,sys, eq,inds; sys := gaia[standardform](Sys,typ); userinfo(1,gaia,`Chomsky Normal Form=`,sys); # step 1 userinfo(2,gaia,`Computing valuations ...`); val:=valuation_CNF(sys,zprods); userinfo(2,gaia,`Looking for an infinite valuation ...`); if has(op(val),infinity) then inds:=indets(op(op(val)),name=identical(infinity)); RETURN(InfiniteValuation({seq(op(1,t),t=inds)})) fi; # step 2 userinfo(2,gaia,`Looking for a production like Sequence(Epsilon) ...`); for t in indets(Sys,function) do Phi := op(0,t); Y := op(t); if member(Phi,MultiConstructors) then if op(2,properties[Phi])=`2a` and val[Y]=0 then userinfo(2,gaia,t,` derives an infinity of objects of size 0`); RETURN(false) fi fi od; N:=0; for eq in sys do t := op(2,eq); if type(t,function) then Phi := op(0,t); Y := op(t); if member(Phi,MultiConstructors) then N := max(N,N_Phi[Phi]) elif member(Phi,{Union}) or type(Phi,product) then elif type(Phi,delta) then N := max(N,1) elif type(Phi,subst) then N := max(N,val[op(1,t)]) else ERROR(`unknown constructor`,Phi) fi; fi od; userinfo(2,gaia,`The dependency graph is invariant up from size`,N); alpha := array(0..N); if N>0 then userinfo(2,gaia,`Looking for loops for size <`,N) fi; # alpha[k] is a valid order between non terminals for size k # step 3 n:=0; while n=`,N, `and`,l,`derives an object with such size`); RETURN(false) else ERROR(`should not happen 2`) fi; fi; alpha[N]:=l; userinfo(3,gaia,`order of types for size`,N,'is',l); eval(val); end: derive := proc(U,n,sys) # does U derive an object of size n ? local val; val := valuation_CNF(sys); if n= n ? local val; val := valuation_CNF(sys); if n<=val[U] then true else ERROR(`not yet implemented`) fi; end: dependency_graph := proc(sys,n) # sys is in CNF local g,eq,Phi,rhs,U,V,W,val; val := valuation_CNF(sys); # needs valuations g := NULL; for eq in sys do rhs := op(2,eq); if type(rhs,function) then U := op(1,eq); Phi := op(0,rhs); if Phi=Union then g := g,seq(V0 then g := g,op(rhs){} do ll := ll minus {op(i,l)} od; if i=j then l := [op(l),U] # insert last else l := subsop(i=op([U,op(i,l)]),l) fi; # insert before l[i] od; l end: # all ProdPrd are part of a larger Prod, so we don't # actually want the word Prod to appear `gaia/print/ProdPrd` := proc(a,b) args; end: `gaia/print/ProdSet` := proc(A,b) local a; if type(A,list) then a:=op(A) else a:=A fi; if b=Epsilon then `gaia/print/Set`(a) elif type(b,specfunc(anything,Set)) then `gaia/print/Set`(a,op(b)) else `gaia/print/Set`(a,b) fi end: `gaia/print/ProdSeq`:=subs(Set=Sequence, `gaia/print/Set`=`gaia/print/Sequence`,op(`gaia/print/ProdSet`)): # what reaches here ought to be the outside level of a cycle, hence we want # the Cycle in front `gaia/print/ProdCyc` := proc(a,b) if type(b,epsilon) then `gaia/print/Cycle`(a) elif type(b,specfunc(anything,Sequence)) then `gaia/print/Cycle`(a,op(b)) else `gaia/print/Cycle`(args); fi end: # these elements are inside a cycle, hence should not get the word "Cycle" # in front `gaia/print/ProdInCyc` := proc(a,b) if type(b,epsilon) then a; elif type(b,specfunc(anything,Sequence)) then a,op(b) else args fi end: `gaia/print/DeltaSet` := proc(s,k) if k=1 then s else [s $ k] fi end: `gaia/print/DeltaCyc` := proc(c,k) if k=1 then c else `gaia/print/Cycle`('op(c)' $ k) fi end: `gaia/print/Prod` := proc() Prod(args) end: `gaia/print/Set` := proc() Set(args) end: `gaia/print/Sequence` := proc() Sequence(args) end: `gaia/print/Cycle` := proc() #local i; # {seq([args[i..nargs],args[1..i-1]],i=1..nargs)}; Cycle(op(op(1,"))) Cycle(args) end: `gaia/print/Epsilon` := () -> Epsilon: `gaia/print/constr`['Prd']:=Prod: `gaia/print/constr`['Set']:='Set': `gaia/print/constr`['Seq']:='Sequence': `gaia/print/constr`['Cyc']:='Cycle': `help/text/gaia` := TEXT(``, `HELP FOR: The Gaia package`, ``, `CALLING SEQUENCE:`, ` (args)`, ` gaia[](args)`, ``, `SYNOPSIS:`, `- The gaia package contains functions that count and generate random`, ` labelled or unlabelled combinatorial objects.`, ``, `- The functions available are:`, ``, ` compile count draw standardform`, ``, `- For help with a particular function do: ?gaia[], where`, ` is taken from the above list.`, ``, `- To see how to define a combinatorial class, see ?gaia[specification].`, ``, `- To draw a random object from a given specification, see ?gaia[draw]`, ``, `- Information about the computations that are being done can be obtained by`, ` setting infolevel['gaia'] to anything between 1 and 5.`, ` `, `SEE ALSO: with, gaia[specification], gaia[variables]`): `help/gaia/text/variables`:=TEXT(``, `HELP FOR: gaia variables`, ``, `CALLING SEQUENCE:`, ` Count[specif,typ][A](n) Draw[specif,typ][A](n,k)`, ``, `SYNOPSIS:`, `- These "global" variables are created by the gaia package. The variable`, ` Count contains counting procedures and Draw contains drawing procedures.`, ` If s is a specification, t its labelling type, A a non terminal of s and`, ` n a non-negative integer, Count[s,t][A](n) gives the number of objects of`, ` size n in A; if k is an integer in the range 0..Count[s,t][A](n)-1, then`, ` Draw[s,t][A](n,k) draws a random object of size n in A.`, ``, `SEE ALSO: gaia[count], gaia[draw], gaia[specification]` ): `help/gaia/text/compile` := TEXT(``, `FUNCTION: gaia[compile] - compile a combinatorial specification`, ``, `CALLING SEQUENCE:`, ` compile(spec,typ)`, ``, `PARAMETERS:`, ` spec - a combinatorial specification`, ` typ - its labelling type ('labelled' or 'unlabelled')`, ``, `SYNOPSIS:`, `- For each non terminal T of the specification, the function produces a`, ` counting procedure and a drawing procedure, that are assigned respec-`, ` tively to Count[spec,typ][T] and Draw[spec,typ][T].`, ``, `EXAMPLES:`, `> with(gaia):`, `> s := { A = Prod(Z, Set(A)) }, labelled:`, `> compile(s);`, `> seq(Count[s][A](n),n=0..10);`, ``, ` 0, 1, 2, 9, 64, 625, 7776, 117649, 2097152, 43046721, 1000000000`, ``, `> draw(s,A,4);`, ``, ` Prod(Z[4], Set(Prod(Z[1], Set(Prod(Z[2], Epsilon))), Prod(Z[3], Epsilon)))`, ``, `SEE ALSO: gaia[draw]` ): `help/gaia/text/draw` := TEXT(``, `FUNCTION: gaia[draw] - draws a random combinatorial object`, ``, `CALLING SEQUENCE:`, ` draw(spec, typ, A, n)`, ``, `PARAMETERS:`, ` spec - a combinatorial specification`, ` typ - its labelling type`, ` A - a non terminal of spec`, ` n - integer specifying the size of the object to be drawn`, ``, `SYNOPSIS:`, `- This function outputs a random object of size n in the class A defined by`, ` the specification 'spec', with uniform distribution among all objects of`, ` same size.`, ``, `EXAMPLES:`, `> with(gaia):`, `> bin := {B=Union(Z,Prod(B,B))}:`, `> draw(bin,labelled,B,7);`, ``, ` Prod(Prod(Prod(Prod(Z[4], Z[1]), Z[2]), Prod(Z[7], Prod(Z[6], Z[5]))), Z[3])`, ``, `> draw(bin,unlabelled,B,7);`, ``, ` Prod(Prod(Z, Prod(Prod(Z, Prod(Z, Z)), Prod(Z, Z))), Z)`, ``, `SEE ALSO: gaia[specification]` ): `help/gaia/text/specification` := TEXT(``, `HELP FOR: specification of a combinatorial class`, ``, `SYNOPSIS:`, `- A combinatorial class either contains only one object, or is built from`, ` simpler classes with "constructors". The elementary classes are Epsilon,`, ` which denotes an object of size zero, or Z, which denotes an object of`, ` size one. The available constructors are:`, ``, ` Atom object of size 1 (Z is a predefined atom)`, ` Union(A, B, ...) disjoint union of the classes A, B, ...`, ` Prod(A, B, ...) partitional product of the classes A, B, ...`, ` Set(A) all sets with repetitions whose elements are in A`, ` Sequence(A) all sequences of elements of A`, ` Cycle(A) all directed cycles of elements of A`, ` Subst(A,B) or B(A) B-objects whose atoms are replaced by A-objects`, ``, `- For the constructors Set, Sequence and Cycle, it is possible to add some`, ` restrictions on the cardinality: for example, Set(A,card>=1) means all`, ` non empty sets whose elements are in A, Sequence(A,card<=3) means all`, ` sequences of at most three elements of A, and Cycle(A,card=5) means all`, ` directed cycles of five elements from A.`, ``, `- A specification is a set of productions of the form A = , where A is`, ` the name of the class being defined, and is an expression involving`, ` elementary classes, constructors and other classes. For example, below are`, ` the specifications of some well-known combinatorial classes:`, ``, ` Labelling type = 'labelled':`, ``, ` {A = Prod(Z,Set(A))} non plane trees`, ` {B = Union(Z,Prod(B,B))} plane binary trees`, ` {C = Prod(Z,Sequence(C))} plane general trees`, ` {D = Set(Cycle(Z))} permutations`, ` {E = Set(Cycle(A)), A=Prod(Z,Set(A))} functional graphs`, ` {F = Set(Set(Z,card>=1))} set partitions`, ` {G = Union(Z,Prod(Z,Set(G,card=3)))} non plane ternary trees`, ` {H = Union(Z,Set(H,card>=2))} hierarchies`, ` {L = Set(Set(Set(Z,card>=1),card>=1))} 3-balanced hierarchies`, ` {M = Sequence(Set(Z,card>=1))} surjections`, ``, ` Labelling type = 'unlabelled':`, ``, ` {A = Set(Sequence(Z,card>=1))} integer partitions`, ` {B = Sequence(Union(Z,Z))} binary sequences`, ` {C = Cycle(Set(Z,card>=1))} necklaces`, ` {D = Prod(Z,Set(D))} rooted unlabelled trees`, ` {E = Set(Cycle(D)), D=Prod(Z,Set(D))} random mappings patterns`, ` {F = Union(Z,Set(F,card=2))} nonplane binary trees`, ` {G = Union(Z,Set(G,card=3))} nonplane ternary trees`, ` {H = Union(Z,Set(H,card>=2))} unlabelled hierarchies`, `` ): `help/gaia/text/count`:=TEXT(``, `FUNCTION: gaia[count] - counts the number of objects of a given size`, ``, `CALLING SEQUENCE:`, ` count(spec, typ, A, n)`, ``, `PARAMETERS:`, ` spec - a combinatorial specification`, ` typ - its labelling type`, ` A - a non terminal of spec`, ` n - integer specifying the size of the object to be drawn`, ``, `SYNOPSIS:`, `- This function outputs the number of objects of size n derived by the non`, ` terminal A in the specification spec.`, ``, `EXAMPLES:`, `> with(gaia):`, `> # number of rooted nonplane unlabelled trees of size 100 is:`, `> count({D=Prod(Z,Set(D))},unlabelled,D,100);`, ``, ` 51384328351659326880337136395054298255277970`, ``, `SEE ALSO: gaia[draw], gaia[specification]` ): `help/gaia/text/standardform`:=TEXT(``, `FUNCTION: gaia[standardform] - gives the standard form of a specification`, ``, `CALLING SEQUENCE:`, ` standardform(spec,typ)`, ``, `PARAMETERS:`, ` spec - a combinatorial specification`, ` typ - its labelling type ('labelled' or 'unlabelled')`, ``, `SYNOPSIS:`, `- This function transforms spec into a standard specification, that is a`, ` specification with only the following constructors: Atom, Union, Prod,`, ` Theta, Int and Delta. Theta is the pointing constructor, Int the formal`, ` inverse of Theta, and Delta is the generalized diagonal.`, ``, `EXAMPLES:`, `> with(gaia):`, `> standardform({A=Sequence(B)},labelled);`, ``, ` {A = Union(Epsilon, T0), T0 = Prod[Seq](B, A)}`, ``, `> standardform({A=Set(B)},labelled);`, ``, ` {T1 = Theta(B), T3 = Int(T2), A = Union(Epsilon, T3), T2 = Prod[Set](T1, A)}`, ``, `> standardform({A=Cycle(B,card>=2)},unlabelled);`, ``, ` {T11 = Union(T8, T10), T6 = Prod[Seq](B, T4), T8 = Prod[Cyc](T7, T6),`, ``, ` T7 = Theta(B), T4 = Union(Epsilon, T6), T9 = Prod[Cyc](T7, T4),`, ``, ` A = Int(T11), T10 = Delta[Cycle(2)](T9)}`, ``, `SEE ALSO: gaia[specification]` ): ############################################################################### holonomic_part := proc(Sys) local l,x,changed,eq,rhs,sys,Phi,hol,r,rnep; sys:=Sys; l := {seq(op(1,x),x=sys)}; do changed:=false; for eq in sys do rhs := op(2,eq); hol := true; if type(rhs,{epsilon,atom}) then elif type(rhs, nonterminal) then elif type(rhs,function) then Phi:=op(0,rhs); if Phi=Union or type(Phi,product) then r:={op(rhs)} minus l; rnep:=select(proc(x) not type(x,epsilon) end,r); if rnep <> {} then hol:=false fi; else hol:=false fi else ERROR(`bad production`,rhs) fi; if hol=false then l := l minus {op(1,eq)}; sys := sys minus {eq}; changed := true; fi od; if not changed then break fi; od; userinfo(2,gaia,`holonomic part of`,Sys,'is',sys); sys end: hol_to_algsys := proc(sys) local res,eq,rhs,x,r; res:=NULL; for eq in sys do rhs := op(2,eq); if type(rhs, nonterminal) then r := rhs elif type(rhs,epsilon) then r:=1 elif type(rhs,atom) then r:='z' elif type(rhs,specfunc(anything,Union)) then r:=convert(rhs,`+`) elif type(rhs,function) and type(op(0,rhs),product) then r:=convert(rhs,`*`) else ERROR(`non holonomic production`,rhs) fi; res:=res,op(1,eq)=r od; res:={res}; subs(seq(x=1,x=indets(res,epsilon)),res) end: algsys_to_eqs := proc(sys) local s,inds,r; s := map(proc(x) if type(x,`=`) then op(1,x)-op(2,x) else x fi end,sys); inds := indets(s,name) minus {z}; solve(s,inds); _EnvFloats:=false; # for version V.2 r:=frontend(`solve/gensys`,[s,inds,{},100],[{`+`,`*`,set},inds]); map(hol_simpl,r); end: hol_simpl := proc(eq) local rhs,r,T,s; if not has(eq,RootOf) then op(1,eq)-op(2,eq) else T := op(1,eq); rhs := simplify(op(2,eq),RootOf); r:=indets(rhs,RootOf); if nops(r)<>1 then ERROR(`too many RootOfs`,r) fi; r:=op(r); s:=subs(r=_Z,rhs-T); s:=resultant(s,op(r),_Z); primpart(s,T); fi end: # these three functions should be replaced by those of gfun # once the fix will be distributed `gaia/minindex` := proc(rec,u,n) local s; s:=select(has,indets(rec,u(anything)),n); if s={} then -infinity else min(op(map(op,s)))-n fi; end: `gaia/maxindex` := proc(rec,u,n) local s; s:=select(has,indets(rec,u(anything)),n); if s={} then infinity else max(op(map(op,s)))-n fi; end: `gaia/difforder` := proc(rec,u,n) `gaia/maxindex`(rec,u,n)-`gaia/minindex`(rec,u,n) end: `makecount/holonomic` := proc(eq) local y,a,n,deq,rec,recorder,b,p,i,inds,r; inds:=indets(eq,name) minus {'z'}; if nops(inds)<>1 then ERROR(`invalid arguments`,eq) else y:=op(inds) fi; if not type(y,{procedure,table}) and eval(subs(INDS=y,'assigned(INDS)')) then ERROR(`assigned variable`,y) fi; if degree(eq,y)=1 then p:=solve(eq,y) fi; if degree(eq,y)>=2 or not type(p,polynom) then userinfo(3,gaia,`algebraic degree of`,y,'is',degree(eq,y)); userinfo(4,gaia,`algebraic equation of`,y,'is',eq); deq:=algfuntodiffeq(numer(eq),y('z')); userinfo(3,gaia,`differential order of`,y,'is', `gfun/dorder`(deq,y('z'))); userinfo(4,gaia,`differential equation of`,y,'is',deq); rec:=diffeqtorec(deq,y('z'),a(n)); if hastype(rec,y(integer)) then rec := subs(seq(y(i)=Count[`gaia/sys`][y](i), i=map(op,indets(rec,y(integer)))),rec) fi; if `gaia/sys`[2]='labelled' then rec := subs(b=a,invborel(rec,a(n),b)) fi; recorder := `gaia/difforder`(rec,a,n); userinfo(3,gaia,`recurrence order for`,y,'is',recorder); userinfo(4,gaia,`recurrence for`,y,'is',rec); if not type(rec,set) then rec:={rec} fi; r:={seq(a(i)=Count[`gaia/sys`][y](i),i=0..recorder-1)}; subs(r,rectoproc(rec union r,a(n))) else # degree=1 and p is a polynom subs(POLY=p,proc(gArg_1) option remember; coeff(POLY,z,gArg_1) end) fi end: gaia[version]:=1.2: save `gaia.m`;quit