#------------------------------------------------------------------------------ | | | The CS package | | | | Version 1.4 | | | | For MuPAD 1.3, 1.4 and 1.4.1 | | | | | | Available from: | | http://dept-info.labri.u-bordeaux.fr/~dutour/CS | | | | Copyright: Sylvie Corteel - LRI | | Alain Denise - LRI | | Isabelle Dutour - LaBRI | | Paul Zimmermann - INRIA | | | | Bugs Report: cs@labri.u-bordeaux.fr | | | ------------------------------------------------------------------------------- Modifications: 1.4 (Mar 99) Now, the C programs respect the ANSI norm for the declaration of functions. Add the option target=MODULE at cs::compile for generating automatically C dynamics modules. First running version for library=GMP for generating automatically C code which uses GMP library instead of DPE. 1.3 (Jul 98) Fixed a bug in the standardForm procedure (same name for different productions). 1.2 (May 98) Add error message when some non-terminals are undefined. Add tests on arguments for speciftoalgeq and speciftorec. All C compilation options are now in capital letters, and the C programs are now compatible also with AIX and LINUX. Implantation of the lazy multiplication algorithm designed by J. Van der Hoeven. Added a global variable cs::LAZY which says if this algorithm is used or not in compile for the Prod constructors. 1.1 (Apr 98) Now treat special case C=Z in cs::makecount2[Prod], both for MuPAD/C target languages (counting cost O(n) instead of O(n^2)). For C programs, the compilation option -Dwords permits to print just the sequence of letters and not the operator Prod. In progress: integration of Joris' algorithm. 1.0 (Mar 98) Fixed a conflict name in speciftoalgeq and algeqtodiffeq (when A was used for a non-terminal). Fixed a bug due to the initialization of the globale variable cs::type ==> standardForm1 takes now the type as argument and initializes cs::type. The linear recurrences are not computed any more for the non-terminals T=Prod(A,B) if A or B is Atom or Epsilon. Several other small modifications before some bigger ones !... 0.92 (Mar 98) The printing post-treatment for Set (without card conditions) is done, in MuPAD and in C. Fixed a bug for the optimizations on counting procedures done in the previous version. Renamed cs::version in cs::VERSION. 0.91 (Feb 98) Sequence(B,card=k) is now supported. The standardForm procedure have been changed in oder to permit the printing post-treatment for Prod and Sequence in MuPAD and in C. For counting procedures, optimizations have been down for not computing several times the same procedure. 0.9 (Jan 98) The call cs::compile(NIL) deletes the whole remember table. Added remember option for all count functions for T=Prod(A,B). The linear recurrences are only computed for the non-terminals T=Prod(A,B) where the both A and B are not Atom or Epsilon. The C programs are compatible also with HP. For the moment, the C programs generated with HOLONOMIC=TRUE produce right calculations only if the recurrences do not implicate negative coefficients. 0.8 (Jan 98) The C programs are compatible with Solaris, SunOs, DecAlpha and IRIX. Fixed a bug in C count functions (Delta). Added remember option for the standard form and added the variable cs::version. 0.7 (Jan 98) For each non-terminal T=Prod(A,B), the count function is changed to use the linear recurrence, also for target=C. 0.6 (Jan 98) Functions produced by rectoproc are now "derecursived" in order to not go out the limit of the stack (MAXDEPTH). 0.5 (Jan 98) Added a global variable cs::HOLONOMIC which says if speciftorec and rectoproc are used or not in compile. Modification of the use of the remember table for compile. Fixed the bug for the number of initial values useful for the recurrences. 0.4 (Dec 97) Added speciftorec and rectoproc (for MuPAD). compile call speciftorec and rectoproc if the specification is context-free. For each non-terminal T=Prod(A,B), the count function is changed to use the linear recurrence. There are yet some questions about the number of initial values useful for the recurrence... 0.3 (Dec 97) First running version for target=C. Modification of the interface for compile, count and draw. The count and draw procedures in MuPAD are stored in two tables indexed by the specifications. 0.2 (Dec 97) Added algeqtodiffeq. Fixed a bug in diffeqtorec. Compatibility with MuPAD 1.3 and 1.4. 0.1 (Dec 97) Added speciftoalgeq and diffeqtorec. 0.0 (Oct 97) First version of compile, count and draw generating procedures in MuPAD and in C (for unlabelled structures). ------------------------------------------------------------------------------- Binary trees: >> sys:={B=Union(Z,Prod(B,B))}: >> cs::compile(sys): >> cs::count([B,sys],size=i)$i=0..10; 0, 1, 1, 2, 5, 14, 42, 132, 429, 1430, 4862 Integer partitions: >> sys:={P=Set(N),N=Sequence(Z,card>=1)}: >> cs::compile(sys): >> cs::count([P,sys],size=i)$i=0..10; 1, 1, 2, 3, 5, 7, 11, 15, 22, 30, 42 Idem with Predefined: >> entier:=fun(1):entier(0):=0: >> sys:={P=Set(N),N=Predefined(entier)}: >> cs::compile(sys): >> cs::draw([P,sys],size=10); Set(N(1), N(1), N(1), N(1), N(3), N(1), N(2)) Motzkin trees: >> sys:={M=Union(Z,Prod(Z,M),Prod(Z,M,M))}: >> cs::compile(sys): >> cs::count([M,sys],size=i)$i=0..10; 0, 1, 1, 2, 4, 9, 21, 51, 127, 323, 835 Schroeder trees: >> sys:={S=Union(Z,Prod(S,U)), U=Union(Prod(Z,S),Prod(S,U))}: >> cs::compile(sys): >> cs::count([S,sys],size=i)$i=0..10; 0, 1, 0, 1, 1, 3, 6, 15, 36, 91, 232 Integer partitions with odd summands: >> isodd:=fun(args(1) mod 2): >> sys:={P=Set(N),N=Predefined(isodd)}: >> cs::compile(sys): >> cs::count([P,sys],size=i)$i=0..10; 1, 1, 1, 2, 2, 3, 4, 5, 6, 8, 10 Binary trees (in C): >> cs::compile({B=Union(Z,Prod(B,B))}, target=C, file="B.c"): >> system("cc -DSUN4 B.c -lm; a.out 10"): Prod(Z,Prod(Prod(Z,Z),Prod(Z,Prod(Prod(Prod(Z,Prod(Z,Z)),Z),Prod(Z,Z))))) Red-blue trees: >> sys:={R=Union(Z,Prod(Z,B,B)),B=Union(Z,Prod(Z,R,R,R))}: >> cs::speciftoalgeq(sys,R); 3 3 3 6 3 z - R + z + 2 R z + R z >> cs::speciftoalgeq(sys,B); 4 2 4 4 4 6 4 z - B + z + 3 B z + 3 B z + B z Motzkin trees: >> sys:={M=Union(Z,Prod(Z,M),Prod(Z,M,M))}: >> eq:=cs::speciftoalgeq(sys,M); 2 z - M + M z + M z >> deq:=cs::algeqtodiffeq(eq,M(z)); 2 3 2 z - M(z) + z M(z) - z D(M)(z) + 2 z D(M)(z) + 3 z D(M)(z) >> r:=cs::diffeqtorec(deq, M(z), v(n)); v(n - 1) - v(n) - n v(n) + 2 (n - 1) v(n - 1) + 3 (n - 2) v(n - 2) or directly: >> r:=cs::speciftorec(sys,M,v(n)); v(n - 1) - v(n) - n v(n) + 2 (n - 1) v(n - 1) + 3 (n - 2) v(n - 2) >> f:=cs::rectoproc(r,v(n)); proc(n) name f; option remember; begin (n*(-1) + (-1))^(-1)*(level(procname)(n + (-1))*(n*2 + (-1)) + level(pro\ cname)(n + (-2))*(n*3 + (-6)))*(-1) end_proc >> cs::LIM; 2 >> f(0):=0: f(1):=1: f(2):=1: f(i)$i=0..10; 0, 1, 1, 2, 4, 9, 21, 51, 127, 323, 835 C code generation: ------------------ WITH DPE LIBRARY ---------------- Generates 5 binary trees of size 10: >> cs::compile({B=Union(Z,Prod(B,B))}, target=C, file="B.c", main=B): >> system("cc -DSUN4 B.c -lm; a.out 10 5"): Prod(Z,Prod(Prod(Z,Z),Prod(Z,Prod(Prod(Z,Z),Prod(Z,Prod(Prod(Z,Z),Z)))))) Prod(Prod(Z,Prod(Prod(Prod(Z,Prod(Z,Z)),Prod(Z,Z)),Prod(Z,Prod(Z,Z)))),Z) Prod(Prod(Z,Prod(Prod(Prod(Z,Prod(Prod(Z,Z),Z)),Prod(Prod(Z,Z),Z)),Z)),Z) Prod(Prod(Prod(Z,Prod(Prod(Prod(Z,Prod(Prod(Z,Z),Z)),Prod(Z,Z)),Z)),Z),Z) Prod(Prod(Prod(Prod(Prod(Z,Prod(Prod(Z,Z),Prod(Z,Z))),Prod(Z,Z)),Z),Z),Z) Compilation on various systems: Solaris: cc -DSUN4 B.c -lm Sunos: cc -DSUNOS B.c -lm Alpha: cc -fprm d -DALPHA B.c -lm IRIX: cc -DIRIX64 B.c -lm HP: cc -DHP700 B.c -lm AIX: cc -DAIX B.c -lm LINUX: cc -DLINUX B.c -lm A regular grammar: (-DWORDS option) >> spec:={S=Union(Epsilon,Prod(a,S,b),Prod(c,S)),a=Atom,b=Atom,c=Atom}: >> cs::compile(spec, target=C, file="S.c", main=S): >> system("cc -DSUN4 S.c -lm; a.out 9"): Prod(c,Prod(a,Prod(a,Prod(a,Prod(c,Prod(c,Epsilon)),b),b),b)) >> system("cc -DSUN4 -DWORDS S.c -lm; a.out 9"): acccccabb Motzkin trees: >> spec:={M=Union(Z,Prod(Z,M),Prod(Z,M,M))}: >> cs::compile(spec, target=C, file="M.c", main=M): Red-blue trees: >> sys:={R=Union(Z,Prod(Z,B,B)),B=Union(Z,Prod(Z,R,R,R))}: >> cs::compile(sys, target=C, file="RB.c", main=R): Integer partitions: >> spec:={P=Set(Sequence(Z,card>=1))}: >> cs::compile(spec, target=C, file="P.c", main=P): >> system("cc -DSUN4 P.c -lm; a.out 4"): Set(Sequence(Z),Sequence(Z,Z,Z)) Binary sequences: >> spec:={S=Sequence(Union(A,B)),A=Atom,B=Atom}: >> cs::compile(spec, target=C, file="S.c", main=S): >> system("cc -DSUN4 S.c -lm; a.out 9"): Sequence(A,A,A,B,B,A,A,A,B) Rooted unlabelled trees: >> cs::compile({A=Prod(Z,Set(A))}, target=C, file="D.c", main=A): >> system("cc -DSUN4 D.c -lm; a.out 3"): Prod(Z,Set(Prod(Z,Epsilon),Prod(Z,Epsilon))) And not yet implemented: Nonplane binary trees: >> spec:={F=Union(Z,Set(F,card=2))}: >> cs::compile(spec, target=C, file="F.c", main=F): Nonplane ternary trees: >> spec:={G=Union(Z,Set(G,card=3))}: >> cs::compile(spec, target=C, file="G.c", main=G): Hierarchies: >> spec:={H=Union(Z,Set(H,card>=2))}: >> cs::compile(spec, target=C, file="H.c", main=H): Necklaces: >> spec:={C=Cycle(Sequence(Z,card>=1))}: >> cs::compile(spec, target=C, file="C.c", main=C): WITH GMP LIBRARY ---------------- >> spec:={T=Prod(Z,Sequence(T))}: >> cs::compile(spec, target=C, library=GMP, file="PTrgmp.c", main=T): >> system("gcc PTrgmp.c -lgmp; a.out 9"): C dynamic modules generation: ----------------------------- WITH DPE LIBRARY ---------------- >> spec:={T=Prod(Z,Sequence(T))}: >> cs::compile(spec, target=MODULE, file="PTrmoddpe.c"): >> system("mmg -DSUN4 PTrmoddpe.c"): >> module("PTrmoddpe"); >> PTrmod::count(T,10); 4.8619999999999957*10^3..4.8620000000000055*10^3 [4862.0, 4862.0] >> PTrmod::draw(T,10); WITH GMP LIBRARY ---------------- not yet implemented # case version() of [1,4,1] do cs:=domain("cs"); DOM_FUN:="fun"; break of [1,4,0] do cs:=domain("cs"); DOM_FUN:="fun"; break of [1,3,0] do cs:=domain("cs"); DOM_FUN:=DOM_EXEC; break otherwise cs:=domain(); DOM_FUN:=DOM_EXEC; end_case: # constants and types # cs::VERSION:="1.4": cs::HOLONOMIC:=FALSE: # using or not linear recurrences # cs::LAZY:=FALSE: # using or not lazy product algorithm # cs::Q:=Dom::Fraction(Dom::DistributedPolynomial([z],Dom::Rational)): cs::MQ:= Dom::Matrix(): # s is a set of productions of the form A = Epsilon A = Atom A = Union(B,C,...) A = Prod(B,C,...) A = Set(B) A = Sequence(B) A = Predefined(function) compile(s [, type=labelled] [, target=C] [, library=GMP] [, file="B.c"] [, main=NT]) - type is either unlabelled (default) or labelled - target is either MuPAD (default) or C or MODULE - library is either DPE (default) or GMP - file (string) is the name of the output file (in C or MODULE mode) - main is a non-terminal (in C mode) # cs::compile := proc(s) local type_option,target_option,library_option,file_option,nt_option,arguments; begin if testargs() then if args(0)=0 then error("invalid number of arguments"); end_if; # delete the remember table of compile1 # if args(1)=NIL then cs::compile1:=subsop(cs::compile1,5=NIL); return(NIL); end_if; if not testtype(s,DOM_SET) then error("invalid grammar specification"); end_if; end_if; # get or set all the options # type_option:=subs(type, args(2..args(0)), type=unlabelled); if type_option=labelled then error("Not yet implemented: ".expr2text(type_option)) end_if; target_option:=subs(target, args(2..args(0)), target=MuPAD); library_option:=subs(library, args(2..args(0)), library=DPE); file_option:=subs(file, args(2..args(0))); nt_option:=subs(main, args(2..args(0)), main=op(op(s,1),1)); arguments:=s,type_option,target_option,library_option,file_option,nt_option; if cs::compile1(arguments) <> [cs::HOLONOMIC,cs::LAZY] then userinfo(1,"new compile for these arguments"); cs::compile1(arguments):=NIL; cs::compile1(arguments); end_if; TRUE: end_proc: cs::compile1 := proc(s,type_option,target_option,library_option,file_option,nt_option) option remember; local sys,p; begin userinfo(1,"compile1 called with ",args(2..args(0))); # index for the tables tCount and tDraw # cs::user_spec_NF:=s, type_option; cs::tCount[cs::user_spec_NF]:=table(); cs::tDraw[cs::user_spec_NF]:=table(); sys:=cs::standardForm1(s,type_option); userinfo(1,"standard form is: ",sys); cs::valuations(sys); # detect non-terminals that will have the same count function # cs::equalcount(sys); if target_option<>MuPAD then target_option:=target_option,library_option; end_if; cs::init[target_option](s,sys,file_option,nt_option); if cs::HOLONOMIC and cs::testcontextfreespecif(sys) then userinfo(1,"using linear recurrences..."); cs::makecount[MuPAD](sys); cs::makecountprod[target_option](sys,HOLONOMIC); elif cs::LAZY then userinfo(1,"using lazy products..."); cs::makecount[MuPAD](sys); cs::makecountprod[target_option](sys,LAZY); else userinfo(1,"using naive products..."); cs::makecount[target_option](sys); end_if; cs::makedraw[target_option](sys); cs::reset[target_option](); return([cs::HOLONOMIC,cs::LAZY]) end_proc: cs::count:=proc() local y, sys, typ, s; begin ([y,sys,typ,s]):=cs::parse_args(args()); cs::compile(sys,type=typ); cs::tCount[sys,typ][y](s) end_proc: cs::draw:=proc() local y, sys, typ, s; begin ([y,sys,typ,s]):=cs::parse_args(args()); cs::compile(sys,type=typ); if cs::tCount[sys,typ][y](s)=0 then error("there is no structure of this size") end_if; cs::tDraw[sys,typ][y](s) end_proc: cs::parse_args:=proc(Spec,Size) local s, spec; begin if args(0)<>2 then error("invalid number of arguments"); else if (testtype(Size,"_equal") and testtype(op(Size,2),DOM_INT) and op(Size,1)=hold(size)) then s:=op(Size,2) else error("invalid size specification ".expr2text(Size)) end_if end_if; spec := Spec; if testtype(spec,DOM_LIST) then if nops(spec)=2 then spec := [op(spec), hold(unlabelled)] end_if; if nops(spec)<>3 then error("wrong number of arguments in grammar specification") end_if; # check that the main non-terminal is really in the grammar (in standard form) # if not has(indets(cs::standardForm1(spec[2],spec[3])),spec[1]) then error("non-terminal ".expr2text(spec[1])." does not appear in grammar specification (in standard form)") end_if; else error("invalid grammar specification ".expr2text(spec)) end_if; [spec[1],spec[2],spec[3],s] end_proc: # p=Epsilon(A) or Atom(A) or Union(B,C,...) or Prod(B,C) or ProdPrd(B,C) or ProdSeq(B,C) or ProdPrdSeq(B,C) or ProdSet(B,C) or ProdPrdSet(B,C) or Theta(B) or Int(B) or Delta(B,u) # cs::typeofoper:=proc(p) begin if contains(cs::prodlist,op(p,0)) then return(Prod) else return(op(p,0)) end_if end_proc: cs::prodlist:={Prod,ProdPrd,ProdSeq,ProdPrdSeq,ProdSet,ProdPrdSet}: cs::next:=table(Prod=ProdPrd, ProdPrd=ProdPrd, ProdSeq=ProdPrdSeq, ProdPrdSeq=ProdPrdSeq): cs::contextfreelist:={Atom,Epsilon,Union,Prod}: # sys is supposed in standard form return FALSE if sys contains other operators than cs::contextfreelist # cs::testcontextfreespecif:=proc(sys) local s; begin s:=map(sys,fun(cs::typeofoper(op(args(1),2)))); if (s union cs::contextfreelist)=cs::contextfreelist then TRUE else FALSE end_if end_proc: # sys is NOT supposed in standard form return FALSE if sys contains other operators than cs::contextfreelist # cs::testalgebraicspecif:=proc(sys) local s; begin s:=map(sys,op,2); while map(s,type)<>{DOM_IDENT} do soper:=map(s,cs::typeofoper) minus {FAIL}; if (soper union cs::contextfreelist)<>cs::contextfreelist then return(FALSE) end_if; s:=map(s,op); end_while; return(TRUE) end_proc: # useful functions for getting the list of indices or values of a table t # cs::listofindices:=proc(t) begin map([op(t)], op, 1) end_proc: cs::listofvalues:=proc(t) begin map([op(t)], op, 2) end_proc: # s is in standard form detect non-terminals (Prod productions) that will have the same count function and then, if cs::lequalcount[A]=B, countA <-- countB # cs::equalcount:=proc(s) local lprod,p,i; begin lprod:=[op(select(s,func((if cs::typeofoper(op(p,2))=Prod then TRUE else FALSE end_if),p)))]; cs::lequalcount:=table(); while nops(lprod)>1 do p:=lprod[1]; i:=2; while i<=nops(lprod) do if {op(op(p,2))}={op(op(lprod[i],2))} then cs::lequalcount[op(lprod[i],1)]:= op(p,1); lprod[i]:=NIL; i:=i-1; end_if; i:=i+1; end_while; lprod[1]:=NIL; end_while; end_proc: # eq is an expression C=f(A,B) of a specification s in standard form returns TRUE if type of f is Prod and if A or B are not Atom or Epsilon returns FALSE if THE BOTH A and B are Atom or Epsilon # cs::isProd_1:=proc(eq,s) local tmp1,tmp2; begin if cs::typeofoper(op(eq,2))<>Prod then return(FALSE) end_if; tmp1:=select(s,has,op(op(eq,2),1)); tmp2:=select(s,has,op(op(eq,2),2)); if (has(tmp1,Atom) or has(tmp1,Epsilon)) and (has(tmp2,Atom) or has(tmp2,Epsilon)) then FALSE else TRUE end_if end_proc: # eq is an expression C=f(A,B) of a specification s in standard form returns TRUE if type of f is Prod and if THE BOTH A and B are not Atom or Epsilon returns FALSE if A or B are Atom or Epsilon # cs::isProd_2:=proc(eq,s) local tmp1,tmp2; begin if cs::typeofoper(op(eq,2))<>Prod then return(FALSE) end_if; tmp1:=select(s,has,op(op(eq,2),1)); tmp2:=select(s,has,op(op(eq,2),2)); if (has(tmp1,Atom) or has(tmp1,Epsilon)) or (has(tmp2,Atom) or has(tmp2,Epsilon)) then FALSE else TRUE end_if end_proc: cs::isProd[HOLONOMIC]:=cs::isProd_1: cs::isProd[LAZY]:=cs::isProd_2: #----------------------------------------------------- Standard Form and Valuations -----------------------------------------------------# # types allowed for nonterminals # cs::types:={DOM_IDENT,DOM_FUNC_ENV}: # puts a specification s (set of productions) into standard form, i.e. only with productions of the form A=

where

=Epsilon(A) or Atom(A) or Union(B,C,...) or Prod(B,C) or Theta(B) or Int(B) or Delta(B,u) The second argument designes the type: labelled or unlabelled. Initialize the global variable cs::type Example: >> cs::standardForm1({M=Union(Z,Prod(Z,M),Prod(Z,M,M))},unlabelled); # cs::standardForm1 := proc(s,typ) option remember; local t,p; begin if testargs() then if args(0)=1 then error("wrong number of arguments"); end_if end_if; cs::type:=typ; cs::aliases:=table(); t:={}; for p in s do if type(p)<>"_equal" or not contains(cs::types,type(op(p,1))) then error("invalid production: ".expr2text(p)) end_if; t:=t union cs::standardForm2(op(p)) end_for; # Z is a predefined Atom # if has(t,Z) then t:=t union {Z=Atom} end_if; # transform A = Atom into A = Atom(A) and A = Epsilon into A = Epsilon(A) # map(t, func((if op(p,2)=Atom then op(p,1)=Atom(op(p,1)) elif op(p,2)=Epsilon then op(p,1)=Epsilon(op(p,1)) else p end_if),p)); end_proc: # returns a set of productions for the standard form of {A=p} # cs::standardForm2 := proc(A,p) local q,i; begin if contains(cs::types,type(p)) then {A=p} elif type(p)="function" then case cs::typeofoper(p) of Union do return(cs::standardForm3([op(p)],{A=p})) of Prod do return(( case nops(p) of 0 do {A=Epsilon}; break of 1 do cs::standardForm2(A,op(p)); break of 2 do cs::standardForm3([op(p)],{A=p},A,p); break otherwise # first look for existing aliases # # not used any more because of the decomposition with ProdPrd, ... to be changed... for i from 2 to nops(p) do if type((q:=cs::aliases[op(p,0)(op(p,i-1..i))]))<>"_index" then return(cs::standardForm2(A,subsop(p,i=q,i-1=null()))) end_if end_for; # cs::standardForm2(A,op(p,0)(op(p,1),cs::next[op(p,0)](op(p,2..nops(p))))) end_case)) of Theta do of Int do if nops(p)<>1 then error("invalid right-hand side: ".expr2text(p)) end_if; return(cs::standardForm3([op(p)],{A=p})) of Delta do if nops(p)<>2 then error("invalid right-hand side: ".expr2text(p)) end_if; return(cs::standardForm3([op(p,1)],{A=p})) otherwise q:=cs::standardForm[op(p,0)]; if contains({DOM_PROC,DOM_FUN},type(q)) then q:=q(A,op(p)); return((if q=p then {A=p} else cs::standardForm2(A,q) end_if)) end_if; error("invalid right-hand side of production: ".expr2text(p)) end_case; else error("invalid right-hand side of production: ".expr2text(p)) end_if; end_proc: # for each p in l, if it is not an identifier, creates a new nonterminal and substitutes p by it in s # cs::standardForm3 := proc(l,s) local p,A,t,m,oldl; begin # first puts everything in l in the form A or f(A,B,...) where A,B,... are identifiers # l:={op(l)}; t:=null(); repeat oldl:=l; l:=map(l,proc(x) begin if type(x)=DOM_IDENT then x else x,op(select([op(x)],testtype,"function")) end_if end_proc); until l=oldl end_repeat; # then sort by increasing size # l:=sort([op(l)],proc(x,y) begin length(x)[] do p:=l[1]; l[1]:=NIL; if contains({Atom,Epsilon},p) or not contains (cs::types,type(p)) then if type(cs::aliases[p])="_index" then A:=genident("T"); cs::aliases[p]:=A else A:=cs::aliases[p] # already an alias # end_if; t:=t,p=A; s:=subs(s,p=A) union cs::standardForm2(A,p); l:=subs(l,p=A); end_if end_while; if args(0)>=4 then cs::aliases[subs(args(4),t)]:=args(3) end_if; # remove trivial expressions # select(s,not bool); end_proc: # standard form for A=Sequence(B,c) # cs::standardForm[Sequence]:=proc(A,B,c) local k; begin case args(0) of 0 do of 1 do break of 2 do return(Union(Epsilon,ProdSeq(B,Sequence1(B)))) otherwise # c=card=k or c=card>=k or card<=k # if type(c)="_equal" then # card=k # k:=op(c,2); if testtype(k,Type::NonNegInt) then c:=Epsilon; while k>1 do c:=ProdPrdSeq(B,c); k:=k-1 end_while; if k>0 then c:=ProdSeq(B,c) end_if; return(c) end_if elif type(c)="_leequal" then if op(c,2)=hold(card) then # k<=card # k:=op(c,1); if testtype(k,Type::NonNegInt) then return(ProdSeq(B$k,Sequence1(B))) end_if elif op(c,1)=hold(card) then # card<=k # k:=op(c,2); if testtype(k,Type::NonNegInt) then c:=Epsilon; while k>1 do c:=Union(Epsilon,ProdPrdSeq(B,c)); k:=k-1 end_while; if k>0 then c:=Union(Epsilon,ProdSeq(B,c)) end_if; return(c) end_if end_if end_if end_case; error("invalid production: ".expr2text(Sequence(B,c))) end_proc: cs::standardForm[Sequence1]:=proc(A,B) begin return(Union(Epsilon,ProdPrdSeq(B,A))) end_proc: # standard form for A=Set(B,c) # cs::standardForm[Set]:=proc(A,B,c) local k,C,i; begin C:=Theta(B); if cs::type=unlabelled then C:=Delta(C,fun(1)) end_if; case args(0) of 0 do of 1 do break of 2 do return(Union(Epsilon,Int(ProdSet(C,Set1(B))))) otherwise # c=card>=k or card<=k # if type(c)="_leequal" then if op(c,2)=hold(card) then # k<=card # k:=op(c,1); if testtype(k,Type::NonNegInt) then error("not yet implemented") end_if elif op(c,1)=hold(card) then # card<=k # k:=op(c,2); if testtype(k,Type::NonNegInt) then error("not yet implemented") end_if end_if end_if end_case; error("invalid production: ".expr2text(Set(B,c))) end_proc: cs::standardForm[Set1]:=proc(A,B) local C; begin C:=Theta(B); if cs::type=unlabelled then C:=Delta(C,fun(1)) end_if; return(Union(Epsilon,Int(ProdPrdSet(C,A)))) end_proc: # A = Predefined(B) # cs::standardForm[Predefined] := fun(Predefined(args(2..args(0)))): # s is in standard form put values in cs::val and cs::max # cs::valuations := proc(s) local v,p,changed,r,A,i,w,t; begin # v[A] <= size(element of A) <= w[A] # for p in s do if op(p,[2,0])=Epsilon then s:=s minus {p}; v[op(p,1)]:=0; w[op(p,1)]:=0 elif op(p,[2,0])=Atom then s:=s minus {p}; v[op(p,1)]:=1; w[op(p,1)]:=1 else v[op(p,1)]:=infinity; w[op(p,1)]:=infinity end_if end_for; p:=indets(s) minus {op(cs::listofindices(v))}; if p<>{} then error("undefined non-terminals: ".expr2text(p)) end_if; repeat changed:=FALSE; for p in s do A:=op(p,1); r:=op(p,2); t:=cs::typeofoper(r); if t=Union then min(v[op(r,i)]$i=1..nops(r)),max(w[op(r,i)]$i=1..nops(r)) elif t=Prod then v[op(r,1)]+v[op(r,2)],w[op(r,1)]+w[op(r,2)] elif contains({Theta,Int},t) then max(1,v[op(r,1)]),w[op(r,1)] elif t=Delta then v[op(r,1)],infinity else if type((p:=cs::valuation[op(r,0)]))=DOM_PROC then p(op(r)),infinity else error("unknown constructor: ".expr2text(r)) end_if end_if; p:=%; if p[1] don't print Epsilon # Sequence(args(1)) else Sequence(args(1..args(0))) end_if end_proc: cs::print[ProdPrdSeq]:=proc() begin if args(args(0))=Epsilon then # only 2 args : A,Epsilon => don't print Epsilon # args(1) else args(1..args(0)) end_if end_proc: cs::print[ProdSet]:=proc() begin if args(args(0))=Epsilon then # don't print Epsilon # Set(args(1..args(0)-1)) else Set(args(1..args(0))) end_if end_proc: cs::print[ProdPrdSet]:=proc() begin if args(args(0))=Epsilon then # don't print Epsilon # args(1..args(0)-1) else args(1..args(0)) end_if end_proc: #----------------------------------------------------- Products => linear recurrences -----------------------------------------------------# cs::makecountprod[MuPAD]:=proc(sys,algo) local p; begin for p in select(sys,cs::isProd[algo],sys) do # test if count function for op(p,1) will not be deduced from another one # if type(cs::lequalcount[op(p,1)])="_index" then cs::makecountprod1[algo](op(p),sys) end_if end_for; # count functions deduced from other ones # for p in cs::listofindices(cs::lequalcount) do cs::tCount[cs::user_spec_NF][p]:=subsop(proc(n) begin B end_proc, 4=hold(cs::tCount)[cs::user_spec_NF][cs::lequalcount[p]](hold(n)), 6=hold(cs::tCount)[cs::user_spec_NF][p]) end_for; end_proc: cs::makecountprod1[HOLONOMIC] := proc(A,rhs,sys) local r,u,n,f,i; begin r:=cs::speciftorec(sys,A,u(n)); userinfo(1,"recurrence for ",A," is: ",expand(r)); f:=cs::rectoproc(r,u(n)); for i from 0 to cs::LIM do f(i):=cs::tCount[cs::user_spec_NF][A](i) end_for; userinfo(1,"initial values: ",op(op(f),5)); (cs::tCount)[cs::user_spec_NF][A]:=subsop(f,6=hold(cs::tCount)[cs::user_spec_NF][A]); end_proc: # s is a set of productions of the form A = Epsilon A = Atom A = Union(B,C,...) A = Prod(B,C,...) NT is a non ternimal # cs::speciftoalgeq := proc(s,NT) local g,vars,matA,V,nb_lines,pNT,pol,linear_depend,ens,i,T,tc; begin if testargs() then if args(0)<2 then error("wrong number of arguments") end_if; if not cs::testalgebraicspecif(s) then error("specification in wrong form") end_if; end_if; # computation of groebner basis # ([g,vars]) := cs::groebner_calculus(s); # successive normal forms up to linear dependencies # matA := [[1]]; V := [1]; # involved terms # nb_lines := 1; pNT := poly(NT,vars,cs::Q); pol := poly(1,vars,cs::Q); linear_depend := FALSE; while not linear_depend do nb_lines := nb_lines + 1; pol := groebner::normalf(pol*pNT,g,DegInvLexOrder); userinfo(2,"polynomial is: ",expr(pol)); ens := {expr(nthterm(pol,i)) $ hold(i)=1..nterms(pol)}; if nops(ens union {op(V)}) < nb_lines then linear_depend := TRUE; end_if; ens := ens minus {op(V)}; for T in ens do V := append(V, T); matA := map (matA, append, 0); end_for; userinfo(2,"involved terms are: ",V); for T in V do tc[T] := 0; end_for; for i from 1 to nterms(pol) do tc[expr(nthterm(pol,i))] := nthcoeff(pol,i); end_for; matA := append(matA, [ seq(tc[T], T=V) ]); end_while; # matrix building and gauss elimination # for i from 1 to nops(matA) do matA[i] := append(matA[i], NT^(i-1)); end_for; matA:=cs::MQ(matA); matA:=linalg::gaussElim(matA); # take the algebraic equation (numer of the right coefficient) # i := 1; while not has(matA[i,i], NT) do i := i+1; end_while; numer(matA[i,i]); end_proc: cs::groebner_calculus := proc(s) local hasA,hasE,vars,g; option remember; begin # get variables # if has(s,Z) then s:=s union {Z=Atom} end_if; hasA := select(s,has,Atom); hasE := select(s,has,Epsilon); s := [op(s minus hasA minus hasE)]; vars := map(s, op, 1); userinfo(1,"vars are: ",vars); # transform into polynomials # # the following is better than subs(hasA,Atom=z) since hasA may contain both A=Atom and A=Atom(A), idem for Epsilon # s := subs(s, map(hasA, fun(op(args(1),1)=z))); s := subs(s, map(hasE, fun(op(args(1),1)=1))); s := subs(s, Union=_plus); s := subs(s, op(map(cs::prodlist,fun(args(1)=_mult)))); s := map(s, fun(op(args(1),2)-op(args(1),1))); userinfo(1,"polynomials are: ",s); s := map(s, poly, vars, cs::Q); # groebner basis # g := groebner::gbasis(s, DegInvLexOrder); userinfo(1,"groebner basis is: ",map(g,expr)); [g,vars]; end_proc: cs::algeqtodiffeq := proc(algeq, yofx) local y,x,pol,g,u,v,deg,deq,Y,d,matA,V,T,i; begin if testargs() then if args(0)<2 then error("wrong number of arguments") end_if; end_if; y := op(yofx,0); x := op(yofx,1); pol := poly(algeq, [y]); ([g,u,v]) := [gcdex(diff(pol,y),pol)]; userinfo(2,"g is: ",expr(g)); if has(expr(g),y) then return (cs::algeqtodiffeq(expr(normal(pol/g)),y(x),inits)); end_if; deg := degree(pol); if deg <= 1 then deq := subs(expr(pol),y=y(x)); else Y[1] := map(divide(-u * g * diff(pol,x), pol, Rem), normal); for d from 2 to deg-1 do Y[d] := map(divide(diff(Y[d-1],x)+diff(Y[d-1],y)*Y[1], pol, Rem), normal); end_for; matA := [[1, 0 $ hold(i)=1..deg-1, 1], [0,1, 0 $ hold(i)=1..deg-2, y(x)]]; V := [y^i $ hold(i)=0..deg-1]; # involved terms # for d from 1 to deg-1 do for T in V do tc[T] := 0; end_for; for i from 1 to nterms(pol) do tc[expr(nthterm(Y[d],i))] := nthcoeff(Y[d],i); end_for; matA := append(matA, [ seq(tc[T], T=V) ]); end_for; # matrix building and gauss elimination # for i from 3 to nops(matA) do matA[i] := append(matA[i], (D@@(i-2))(y)(x)); end_for; matA:=cs::MQ(matA); matA:=linalg::gaussElim(matA); # take the differential equation (numer of the right coefficient) # i := 1; while not has(matA[i,i], y) do i := i+1; end_while; deq := numer(matA[i,i]); end_if; deq; end_proc: cs::diffeqtorec := proc(deq, yofx, uofn) begin if testargs() then if args(0)<3 then error("wrong number of arguments") end_if; end_if; deq:=expand(deq); # minimal value of n for which the recurrence is valide # cs::LIM:=degree(select(deq,not has,op(yofx,0)),op(yofx,1))+1; case type(deq) of DOM_INT do of DOM_IDENT do of "function" do of "_mult" do break of "_plus" do deq:=op(deq); break otherwise error("invalid argument") end_case; _plus(op(map([deq],cs::diffeqtorec2,op(yofx,0..1),op(uofn,0..1)))) end_proc: # t is of the form c*x^j*(D@@k)(y)(x), returns the coefficient of x^n in the g.f. of t in terms of u(n)=[x^n] y(x) # cs::diffeqtorec2 := proc(t, y, x, u, n) local c,k,j,i; begin if type(t)="_mult" then ([t,c,k]):=split(t,has,{x,y}) else c:=1 end_if; j:=degree(t,[x]); t:=t/x^j; # now t = 1 or y(x) or D(y)(x) or (D@@k)(y)(x) # if not has(t,y) then 0 else # t = y(x) or D(y)(x) or D(D(y))(x) or ... # if type(t)<>"function" or [op(t)]<>[x] then error("unexpected term: ".expr2text(t)) end_if; t:=op(t,0); # now y or D(y) or D(D(y)) or ... # k:=0; while type(t)="D" do k:=k+1; t:=op(t) end_while; if t<>y then error("unexpected term: ".expr2text(t)) end_if; c*_mult(n+i-j$hold(i)=1..k)*u(n+k-j) end_if end_proc: # speciftorec(s,N,u(n)) computes the linear recurrence for the coefficients u(n) of the generating function associated to the nonterminal N in the (context-free) specification s. Example: >> cs::speciftorec({B = Union(Z, Prod(B, B))}, B, u(n)); u(n)*(-2) + n*u(n)*4 + (n + 1)*u(n + 1)*(-1) # cs::speciftorec := proc(s,N,uofn) local algeq,deq; begin if testargs() then if args(0)<3 then error("wrong number of arguments") end_if; if not cs::testalgebraicspecif(s) then error("specification in wrong form") end_if; end_if; algeq:=cs::speciftoalgeq(s,N); userinfo(1,"algeq is: ",algeq); deq:=cs::algeqtodiffeq(algeq,N(z)); userinfo(1,"diffeq is: ",deq); cs::diffeqtorec(deq,N(z),uofn) end_proc: # returns the smallest i such that u(n+i) appears in a recurrence # cs::minindex:= proc(r,u,n) begin min(op(map(indets(r,PolyExpr) minus {n},fun(op(args(1))-n)))); end_proc: # returns the largest i such that u(n+i) appears in a recurrence # cs::maxindex := proc(r,u,n) begin max(op(map(indets(r,PolyExpr) minus {n},fun(op(args(1))-n)))); end_proc: # r is a linear recurrence with polynomial coefficients and its unknown u, n names to be assigned the unknown function and variable Ouput: a list of polynomials in n: [b(n),p_0(n),...,p_d(n)] meaning r=p_0(n)u(n)+...+p_d(n)u(n+d)+b(n) # cs::formatrec:=proc(r,u,n) local mi,ma,cr,i; begin mi:=cs::minindex(r,u,n); ma:=cs::maxindex(r,u,n); cr:=collect(r,[u(n+i)$i=mi..ma],normal); if mi<>0 then cr:=subs(cr,n=n-mi); ma:=ma-mi; cs::LIM:=cs::LIM+mi; end_if; # I don't know if this test is really useful... # if has(map({op(r)},denom),n) then cr:=collect(numer(normal(cr)),[u(n+i)$i=0..ma],normal) end_if; [subs(cr,[u(n+i)=0$i=0..ma]), coeff(cr,u(n+i),1)$i=0..ma]; end_proc: cs::rectoproc:=proc(r,uofn) local u,n,R,ORDER,lroots,i,FF,Rec; begin u:=op(uofn,0); n:=op(uofn,1); R:=cs::formatrec(r,u,n); ORDER:=nops(R)-2; R:=subs(R,n=n-ORDER); cs::LIM:=cs::LIM+ORDER; userinfo(1,"normalized coefficients are: ",R); # list of integer solutions of the equation coeff(u(n))=0 # lroots:=op(sharelib::iroots(R[nops(R)],n)); # we will have to compute initial values up to this limit # userinfo(1,"LIM,ORDER,ROOTS: ",cs::LIM, ORDER, lroots); cs::LIM:=max(cs::LIM, ORDER, lroots); R:=- _plus(R[1], R[i+2]*u(n-ORDER+i)$i=0..ORDER-1) / R[nops(R)]; subs(proc(n) local nmax,j,g,ff; begin nmax := max(op(cs::listofindices(op(FF,5)))); for j from nmax+1 to n do ff:=FF; g:=procname(j); evalassign(g,Rec,1); end_for; end_proc, FF=hold(level(procname)), Rec=subs(R,u=hold(ff),n=hold(j)) ); end_proc: #----------------------------------------------------- Products => lazy algorithm of Joris VdH -----------------------------------------------------# # n is an integer return p such that p=2^k and 2^(k-1) < n <= 2^k # cs::poweroftwohi:=proc(n) local p; begin p:=1; while n>p do p:=p*2 end_while; p end_proc: cs::makecountprod1[LAZY] := proc(A,rhs,sys) local f,i,p; begin # initialisation of the lazy-incremental table # # cs::tProdlazy[cs::user_spec_NF][A]:=table(i=0$i=0..2000);# p:=cs::poweroftwohi(cs::val[A]); cs::tProdlazy[cs::user_spec_NF][A]:=table(i=0$i=0..2*p); f:=cs::makecountprod2[LAZY](op(rhs),A); for i from 0 to cs::val[A]-1 do f(i):=0 end_for; (cs::tCount)[cs::user_spec_NF][A]:=subsop(f,6=hold(cs::tCount)[cs::user_spec_NF][A]) end_proc: cs::makecountprod2[LAZY] := proc(B,C) # A = Prod(B,C) # begin subs(( proc(n) local nmax,t,d,i; option remember; begin # userinfo(1,"A Prod: ",AA,n); userinfo(1,"B: ",BB,n-1); cB(n-1); userinfo(1,"C: ",CC,n-1); cC(n-1); userinfo(1,"rec A: ",AA,n-1); cA(n-1); # userinfo(1,"A Prod: ",AA,n); nmax := max(op(cs::listofindices(op(cA,5)))); for i from nmax+1 to n-1 do userinfo(1,"rec A: ",AA,i); cA(i) end_for; userinfo(1,"B: ",BB,n-1); cB(n-1); userinfo(1,"C: ",CC,n-1); cC(n-1); if n=0 then T0[0]:=cB(0)*cC(0); return(T0[0]) end_if; if cB(0)<>0 then T0[n]:=T0[n] + cB(0)*cC(n); end_if; if cC(0)<>0 then T0[n]:=T0[n] + cB(n)*cC(0); end_if; if n>1 then # n=(k+1)*2^p, t=2^p and d=2^(p+1) ; Begin with p=0 # t:=1: d:=2: repeat userinfo(2,"Prod t,d: ",t,d); # P(2^p,k*2^p) # cs::karatsuba(0,cB,t,cC,n-t,t); for i from 0 to d-2 do T0[n+i]:=T0[n+i]+cs::dest[i]; end_for; # k=1 and t=n-t=2^p # if t=n-t then break end_if; # P(k*2^p,2^p) # cs::karatsuba(0,cB,n-t,cC,t,t); for i from 0 to d-2 do T0[n+i]:=T0[n+i]+cs::dest[i]; end_for; t:=t*2; d:=d*2; until not cs::iseven(n/(t/2)) end_repeat; end_if; T0[n]; end_proc ), hold(AA)=args(3),hold(BB)=B,hold(CC)=C, hold(cB)=hold(cs::tCount)[cs::user_spec_NF][B], hold(cC)=hold(cs::tCount)[cs::user_spec_NF][C], hold(T0)=hold(cs::tProdlazy)[cs::user_spec_NF][args(3)], hold(cA)=hold(level(procname)) ) end_proc: cs::makecountprod2[LAZY] := proc(B,C) # A = Prod(B,C) # begin subs(( proc(n) local nmax,nmin,l,i,f,FF; begin f:=proc(m) local t,d,i; begin cB(m-1); cC(m-1); if m=0 then T0[0]:=cB(0)*cC(0); return(T0[0]) end_if; if cB(0)<>0 then T0[m]:=T0[m] + cB(0)*cC(m); end_if; if cC(0)<>0 then T0[m]:=T0[m] + cB(m)*cC(0); end_if; if m>1 then # m=(k+1)*2^p, t=2^p and d=2^(p+1) ; Begin with p=0 # t:=1: d:=2: repeat userinfo(2,"Prod t,d: ",t,d); # P(2^p,k*2^p) # cs::karatsuba(0,cB,t,cC,m-t,t); for i from 0 to d-2 do T0[m+i]:=T0[m+i]+cs::dest[i]; end_for; # k=1 and t=m-t=2^p # if t=m-t then break end_if; # P(k*2^p,2^p) # cs::karatsuba(0,cB,m-t,cC,t,t); for i from 0 to d-2 do T0[m+i]:=T0[m+i]+cs::dest[i]; end_for; t:=t*2; d:=d*2; until not cs::iseven(m/(t/2)) end_repeat; end_if; T0[m]; end_proc: # expand the table T0 if nessecary # l:=cs::listofindices(T0); nmax:=max(op(l)); (T0[i]:=0) $ hold(i)=nmax+1..2*cs::poweroftwohi(n); # principal loop # nmax := max(op(cs::listofindices(op(cA,5)))); for i from nmax+1 to n do FF:=procname(i); evalassign(FF,f(i),1); end_for; # garbage collecting of the table T0 # nmin:=min(op(l)); (T0[i]:=NIL) $ hold(i)=nmin..n-1; T0[n]; end_proc ), hold(AA)=args(3),hold(BB)=B,hold(CC)=C, hold(cB)=hold(cs::tCount)[cs::user_spec_NF][B], hold(cC)=hold(cs::tCount)[cs::user_spec_NF][C], hold(T0)=hold(cs::tProdlazy)[cs::user_spec_NF][args(3)], hold(cA)=hold(level(procname)) ) end_proc: cs::iseven:=proc(n) begin if type(n)=DOM_INT and (n mod 2)=0 then TRUE else FALSE end_if end_proc: cs::dest:=table(): # Karatsuba algorithm: n always a power of 2 cs::dest will contain the result of the product destzero: 1st indice of the result in cs::dest t1 and t2: functions to be multiplied t1zero and t2zero: indices from which the multiplication starts n: number of coefficients taken into account example: f1:=proc() begin end_proc: f1:=subsop(f1,5=table(0=1,1=2,2=3,3=4,4=1,5=2,6=3,7=4)): f2:=proc() begin end_proc: f2:=subsop(f2,5=table(0=5,1=6,2=7,3=8,4=9,5=10,6=11,7=12)): cs::karatsuba(0,f1,0,f2,0,8): cs::dest[i]$hold(i)=0..14; expand((1+2*z+3*z^2+4*z^3+1*z^4+2*z^5+3*z^6+4*z^7)*(5+6*z+7*z^2+8*z^3+9*z^4+10*z^5+11*z^6+12*z^7)); # cs::karatsuba:=proc(destzero,t1,t1zero,t2,t2zero,n) local i,h,somme1,somme2,magic; begin userinfo(2,"Begin Karatsuba: ",destzero,t1,t1zero,t2,t2zero,n); if n=1 then cs::dest[destzero]:=t1(t1zero)*t2(t2zero) else h:=n/2; cs::karatsuba(destzero,t1,t1zero,t2,t2zero,h); cs::karatsuba(destzero+n,t1,t1zero+h,t2,t2zero+h,h); cs::dest[destzero+n-1]:=0; for i from 0 to h-1 do somme1(i):=t1(t1zero+i)+t1(t1zero+i+h); somme2(i):=t2(t2zero+i)+t2(t2zero+i+h); end_for; # indice sufficiently large to not disturb the rest of the calculus # magic:=destzero+2*n; cs::karatsuba(magic,somme1,0,somme2,0,h); for i from 0 to n-2 do cs::dest[magic+i]:=cs::dest[magic+i] -cs::dest[destzero+i] -cs::dest[destzero+i+n]; end_for; for i from 0 to n-2 do cs::dest[destzero+i+h]:=cs::dest[destzero+i+h]+cs::dest[magic+i]; end_for; end_if: userinfo(2,"End Karatsuba: ",destzero,t1,t1zero,t2,t2zero,n); end_proc: #----------------------------------------------------- Count and Draw Functions for MuPAD target=MuPAD -----------------------------------------------------# #----------------------------------------------------- C code generation : initialisations -----------------------------------------------------# alias(Print=cs::fprint): cs::fp:=0: cs::fprint:=fun(fprint(Unquoted,cs::fp,args())): cs::filehead:=proc(s,sys,fname,lib) begin PRETTY_PRINT:=FALSE; if type(fname)=DOM_STRING then cs::fp:=fopen(Text,fname,Write); end_if; Print("/* This (".lib.") C source file was automatically generated by CS version ".cs::VERSION.",\n with cs::HOLONOMIC = ".expr2text(cs::HOLONOMIC)." and cs::LAZY = ".expr2text(cs::LAZY).",\n from the following specification:\n ".expr2text(s).",\n in the standard form:\n ".expr2text(sys).".\n You can get CS from http://dept-info.labri.u-bordeaux.fr/~dutour/CS.\n Please report any bug, improvement to the authors of CS:\n cs@labri.u-bordeaux.fr\n Sylvie Corteel - LRI\n Alain Denise - LRI\n Isabelle Dutour - LaBRI\n Paul Zimmermann - INRIA\n*/\n"); end_proc: cs::filetail:=proc() begin fclose(cs::fp); cs::fp:=0; PRETTY_PRINT:=TRUE end_proc: cs::error:=proc(target,lib,st) begin cs::reset[target,lib](); error(st) end_proc: cs::printwords:=proc() begin Print("#ifndef WORDS"); Print("#define PROD_OPEN printf(\"Prod(\")"); Print("#define PROD_SEP printf(\",\")"); Print("#define PROD_CLOSE printf(\")\")"); Print("#define EPSILON printf(\"Epsilon\")"); Print("#else"); Print("#define PROD_OPEN {}"); Print("#define PROD_SEP {}"); Print("#define PROD_CLOSE {}"); Print("#define EPSILON {}"); Print("#endif\n"); end_proc: cs::printdpe:=proc() begin Print("#include \n#include \n#include \n#ifdef SUN4\n#include \n#define TOINFP fpsetround(FP_RP)\n#define TOINFM fpsetround(FP_RM)\n#define Double double\n#elif IRIX64\n#include \n#define TOINFP swapRM(ROUND_TO_PLUS_INFINITY)\n#define TOINFM swapRM(ROUND_TO_MINUS_INFINITY)\n#define Double double\n#elif ALPHA\n#include \n#define TOINFP write_rnd(FP_RND_RP)\n#define TOINFM write_rnd(FP_RND_RM)\n#define Double volatile double\n#elif AIX\n#include \n#define TOINFP fp_swap_rnd(FP_RND_RP)\n#define TOINFM fp_swap_rnd(FP_RND_RM)\n#define Double double\n#elif SUNOS\n#include \nchar *out;\n#define TOINFP ieee_flags(\"set\",\"direction\",\"positive\",&out)\n#define TOINFM ieee_flags(\"set\",\"direction\",\"negative\",&out)\n#define Double double\n#elif HP700\n#define TOINFP fpsetround(FP_RP)\n#define TOINFM fpsetround(FP_RM)\n#define Double double\n#elif LINUX\n#include \n#define TOINFP __setfpucw(0x1b7f)\n#define TOINFM __setfpucw(0x177f)\n#define Double double\n#endif\n\ntypedef struct {\n Double lo,hi;\n int e;\n} DPE;\n\n#define MAXPREC 16 /* IEEE 754 mantissa length in digits */\n\nDouble eps;\n\nvoid DPEinit()\n{\n eps=1.0; while (eps/2.0 != 0.0) eps/=2.0;\n}\n\nvoid DPEprint(DPE *x)\n{\n printf(\"%1.16f*10^%d..%1.16f*10^%d\\n\",x->lo,x->e,x->hi,x->e);\n}\n\n/* x <- y */\nvoid DPEset(DPE *x, double y)\n{\n x->lo=x->hi=y; x->e=0;\n}\n\n/* x <- y */\nvoid DPEcopy(DPE *x, DPE *y)\n{\n x->lo=y->lo; x->hi=y->hi; x->e=y->e;\n}\n\nvoid DPEnormal(DPE *x)\n{\n if (fabs(x->hi)>=10.0) {\n while (fabs(x->hi)>=10.0) {\n TOINFM;\n x->lo /= 10;\n TOINFP;\n x->hi /= 10;\n x->e++; \n }\n } else\n while (x->lo!=0.0 && fabs(x->lo)<=1.0) {\n TOINFM;\n x->lo *= 10;\n TOINFP;\n x->hi *= 10;\n x->e--; \n }\n}\n\n/* x <- y*z, ONLY for y and z positive */\nvoid DPEmul(DPE *x, DPE *y, DPE *z)\n{\n TOINFM;\n x->lo = y->lo * z->lo;\n TOINFP;\n x->hi = y->hi * z->hi;\n x->e = y->e + z->e;\n DPEnormal(x);\n}\n\n/* x <- y*r */\nvoid DPEmulr(DPE *x, DPE *y, double r)\n{\n Double tmp;\n if (r<0) {\n TOINFM;\n tmp = x->lo;\n x->lo = r*x->hi;\n TOINFP;\n x->hi = r*tmp;\n }\n else {\n TOINFM;\n x->lo = r*y->lo; \n TOINFP;\n x->hi = r*y->hi; \n }\n x->e = y->e;\n DPEnormal(x);\n}\n\n/* x <- n*x */\nvoid DPEmuli(DPE *x, long n)\n{\n Double tmp;\n if (n<0) {\n TOINFM;\n tmp = x->lo;\n x->lo = n*x->hi;\n TOINFP;\n x->hi = n*tmp;\n }\n else {\n TOINFM;\n x->lo *= n;\n TOINFP;\n x->hi *= n;\n }\n DPEnormal(x);\n}\n\n/* x <- x/n */\nvoid DPEdivi(DPE *x, long n)\n{\n Double tmp;\n if (n<0) {\n TOINFM;\n tmp = x->lo;\n x->lo = x->hi/n;\n TOINFP;\n x->hi = tmp/n;\n }\n else {\n TOINFM;\n x->lo /= n;\n TOINFP;\n x->hi /= n;\n }\n DPEnormal(x);\n}\n\n/* x <- x+y, destructive on y */\nvoid DPEadd(DPE *x, DPE *y)\n{\n int d=x->e - y->e;\n if (d>0) {\n if (d>MAXPREC) { y->lo=0.0; y->hi=eps; } /* y too small */\n else while (d-->0) {\n TOINFM;\n y->lo /= 10.0;\n TOINFP;\n y->hi /= 10.0;\n }\n }\n else if (d<0) {\n if (d<-MAXPREC) /* x too small */\n { x->lo=y->lo; x->hi=y->hi; y->lo=0.0; y->hi=eps;} \n else while (d++<0) {\n TOINFM;\n x->lo /= 10.0; \n TOINFP;\n x->hi /= 10.0;\n }\n x->e = y->e; \n }\n else {} /* d=0 */\n TOINFM;\n x->lo += y->lo;\n TOINFP;\n x->hi += y->hi;\n DPEnormal(x);\n}\n\nint DPEcmp(DPE *x, DPE *y) /* supposes normalized and positive */\n{\n int d=x->e - y->e;\n if (d>0) return(1);\n else if (d<0) return(-1);\n else /* d=0 */\n if (x->lo>y->hi) return(1);\n else if (x->hilo) return(-1);\n else return(0);\n}\n\n"); Print("extern double drand48();\n"); Print("#define Nan (-1.0)"); Print("#define NOINIT(A) ((A)->lo==Nan)"); Print("#define INIT(A) ((A)->lo!=Nan)\n"); end_proc: cs::printgmp:=proc() begin Print("#include \n#include \n#include "); Print("#include \"gmp.h\"\n"); Print("extern double drand48();\n"); Print("int NOINIT(mpz_t p) { return(mpz_cmp_si(p,-1)); }"); Print("int INIT(mpz_t p) { return(mpz_cmp_si(p,0)); }\n"); end_proc: #----------------------------------------------------- Count and Draw Functions for C with DPE target=C, library=DPE -----------------------------------------------------# cs::init[C,DPE]:=proc(s,sys,fname,main) local A,eq; begin print(Unquoted,"CS warning: ".main." is taken as main non-terminal for C code."); cs::filehead(s,sys,fname,"DPE"); cs::printdpe(): cs::printwords(): Print("DPE s,tmp;"); for eq in sys do A:=expr2text(op(eq,1)); Print("DPE *".A.",*count".A."(int);"); Print("void draw".A."(int);"); end_for; if cs::HOLONOMIC and cs::testcontextfreespecif(sys) then for A in select(sys,cs::isProd[HOLONOMIC],sys) do if type(cs::lequalcount[op(A,1)])="_index" then Print("void initrec".op(A,1)."(int);"); end_if; end_for; end_if; Print("\n"); cs::makemain[C,DPE](sys,main); end_proc: cs::reset[C,DPE]:=cs::filetail: cs::makemain[C,DPE]:=proc(sys,main) local A; begin Print("main(int argc,char *argv[]) {"); Print(" int i,n,k;"); Print(" if (argc<=1) {"); Print(" printf(\"Usage: a.out n [k seed] to generate k objects of size n\\n\");"); Print(" exit(1);}"); Print(" n=atoi(argv[1]);"); Print(" k=(argc>=3) ? atoi(argv[2]) : 1;"); Print(" srand48((argc>=4) ? atoi(argv[3]) : getpid());"); Print(" DPEinit();"); for A in sys do A:=op(A,1); # test if count function for A will not be deduced from another one # if type(cs::lequalcount[A])="_index" then Print(" ".A."=(DPE*)malloc((n+1)*sizeof(DPE));"); Print(" for (i=0;i<".cs::val[A]." && i<=n;i++) DPEset(".A."+i,0.0);"); Print(" for (;i<=n;i++) ".A."[i].lo=Nan;"); end_if; end_for; # count functions deduced from other ones # for A in cs::listofindices(cs::lequalcount) do Print(" ".A."=".expr2text(cs::lequalcount[A]).";"); end_for; if cs::HOLONOMIC and cs::testcontextfreespecif(sys) then for A in select(sys,cs::isProd[HOLONOMIC],sys) do if type(cs::lequalcount[op(A,1)])="_index" then Print(" initrec".op(A,1)."(n);"); end_if; end_for; end_if; A:=main; Print(" count".A."(n);"); Print(" if (".A."[n].lo<=0.0) {printf(\"there is no such structure of this size\\n\"); exit(1);}"); Print(" if (k==0) DPEprint(".A."+n);"); Print(" for (i=1;i<=k;i++) { draw".A."(n); putchar('\\n'); }"); Print("}\n"); end_proc: cs::makecount[C,DPE] := proc(sys) local p; begin for p in sys do cs::makecount1[C,DPE](op(p)); end_for; end_proc: cs::makecount1[C,DPE] := proc(A,rhs) local f,i,NT; begin Print("/* ".expr2text(A=rhs)." */"); Print("DPE* count".A."(int n) { int k;"); # count function for A is the same as NT ? # NT:=cs::lequalcount[A]; if type(NT)<>"_index" then Print(" return(count".expr2text(NT)."(n));"); else # otherwise # Print(" if (NOINIT(".A."+n)) {"); f:=cs::makecount2[C,DPE,cs::typeofoper(rhs)]; if not contains({DOM_FUN,DOM_PROC},type(f)) then error("no C counting rule for constructor ".op(rhs,0)) end_if; f(op(rhs),A); Print(" }"); Print(" return(".A."+n);"); end_if; Print("}\n"); end_proc: cs::makedraw[C,DPE] := proc(sys) local p; begin for p in sys do cs::makedraw1[C,DPE](op(p)); end_for; end_proc: cs::makedraw1[C,DPE] := proc(A,rhs) local f,i; begin Print("/* ".expr2text(A=rhs)." */"); Print("void draw".A."(int n) { int k,c;"); f:=cs::makedraw2[C,DPE,op(rhs,0)]; if not contains({DOM_FUN,DOM_PROC},type(f)) then error("no C drawing rule for constructor ".op(rhs,0)) end_if; f(op(rhs),A); Print("}\n"); end_proc: # Epsilon # cs::makecount2[C,DPE,Epsilon] := proc() begin Print(" DPEset(".args(2)."+n,(n==0) ? 1.0 : 0.0);") end_proc: cs::makedraw2[C,DPE,Epsilon]:=proc() begin Print(" EPSILON;") end_proc: # Atom # cs::makecount2[C,DPE,Atom] := proc() begin Print(" DPEset(".args(2)."+n,(n==1) ? 1.0 : 0.0);") end_proc: cs::makedraw2[C,DPE,Atom]:=proc(A) begin Print(" printf(\"".A."\");") end_proc: # Union # cs::makecount2[C,DPE,Union] := proc() local i,k; begin k:=args(0); Print(" DPEcopy(".args(k)."+n,count".args(1)."(n));"); for i from 2 to k-1 do Print(" DPEcopy(&tmp,count".args(i)."(n));"); Print(" DPEadd(".args(k)."+n,&tmp);"); end_for end_proc: cs::makedraw2[C,DPE,Union] := proc() local i,k; begin k:=args(0); Print(" DPEmulr(&s,count".args(k)."(n),drand48());"); for i from 1 to k-1 do if i>1 then Print(" else {") end_if; Print(" DPEcopy(&tmp,count".args(i)."(n));"); Print(" DPEadd(&s,&tmp);"); Print(" if ((c=DPEcmp(&s,".args(k)."+n))==0) {fprintf(stderr,\"FAIL\\n\"); exit(1);}"); Print(" else if (c==1) { draw".args(i)."(n); return; }"); end_for; Print(" else {fprintf(stderr,\"FAIL\\n\"); exit(1);}"); Print(" "._concat("}"$(k-2))) end_proc: # B=Int(A) ==> B(n) = A(n)/n # cs::makecount2[C,DPE,Int] := proc(A,B) begin Print(" DPEcopy(".B."+n,count".A."(n));"); Print(" DPEdivi(".B."+n,n);"); end_proc: # B=Int(A) ==> B(n) = A(n)/n # cs::makedraw2[C,DPE,Int] := proc(A,B) begin Print(" draw".A."(n);"); end_proc: # B=Theta(A) ==> B(n) = n*A(n) # cs::makecount2[C,DPE,Theta] := proc(A,B) begin Print(" DPEcopy(".B."+n,count".A."(n));"); Print(" DPEmuli(".B."+n,n);"); end_proc: cs::makedraw2[C,DPE,Theta] := proc(A,B) begin Print(" draw".A."(n);"); end_proc: # A = Delta(B,u) # cs::makecount2[C,DPE,Delta]:=proc(B,u,A) begin if u=fun(1) then Print(" DPEcopy(".A."+n,count".B."(n));"); Print(" for (k=2;k<=n;k++)"); Print(" if (n%k==0) {"); Print(" DPEcopy(&tmp,count".B."(n/k));"); Print(" DPEadd(".A."+n,&tmp);"); Print(" }"); else error("not yet implemented for u(n)<>1") end_if end_proc: cs::makedraw2[C,DPE,Delta]:=proc(B,u,A) begin if u=fun(1) then Print(" DPEmulr(&s,count".A."(n),drand48());"); Print(" for (k=1;k<=n;k++)"); Print(" if (n%k==0) {"); Print(" DPEcopy(&tmp,count".B."(n/k));"); Print(" DPEadd(&s,&tmp);"); Print(" if ((c=DPEcmp(&s,".A."+n))==0) {fprintf(stderr,\"FAIL\\n\"); exit(1);}"); Print(" else if (c==1) { draw".B."(n/k); if (k>1) printf(\"$%u\",k); return; }"); Print(" }"); Print(" printf(\"Error: you should not get there\\n\"); exit(1);"); else error("not yet implemented for u(n)<>1") end_if end_proc: #----------------------------------------------------- Naive Products -----------------------------------------------------# cs::makecount2[C,DPE,Prod] := proc(B,C,A) local k0,s0; begin if cs::val[C]=cs::max[C] then # exchange B and C, as counting is symetrical # B; B:=C; C:=%2 end_if; k0:=cs::val[B]; s0:=(if k0=0 then "" else "-".k0 end_if); Print(" DPEmul(".A."+n,count".B."(".k0."),count".C."(n".s0."));"); if k0<>cs::max[B] then s0:=(if cs::val[C]=0 then "" else "-".cs::val[C] end_if); Print(" for (k=".(k0+1).";k<=n".s0.";k++) {"); Print(" DPEmul(&tmp,count".B."(k),count".C."(n-k));"); Print(" DPEadd(".A."+n,&tmp);"); Print(" }"); end_if end_proc: cs::makedraw2[C,DPE,Prod] := proc(B,C,A) local k0,k1,s0; begin k0:=cs::val[B]; k1:=cs::val[C]; if k0=cs::max[B] then s0:=(if k0=0 then "" else "-".k0 end_if); Print(" PROD_OPEN; draw".B."(".k0."); PROD_SEP; draw".C."(n".s0."); PROD_CLOSE;") elif k1=cs::max[C] then s0:=(if k1=0 then "" else "-".k1 end_if); Print(" PROD_OPEN; draw".B."(n".s0."); PROD_SEP; draw".C."(".k1."); PROD_CLOSE;") else Print(" DPEmulr(&s,count".A."(n),drand48());"); Print(" for (k=".cs::val[B].";;k++) {"); Print(" DPEmul(&tmp,count".B."(k),count".C."(n-k));"); Print(" DPEadd(&s,&tmp);"); Print(" if ((c=DPEcmp(&s,".A."+n))==0) {fprintf(stderr,\"FAIL\\n\"); exit(1);}"); Print(" else if (c==1) break;"); k0:=k0-cs::val[C]; s0:=(if k0=0 then "" elif k0>0 then "+".k0 else "-".(-k0) end_if); Print(" k=n".s0."-k;"); Print(" DPEmul(&tmp,count".B."(k),count".C."(n-k));"); Print(" DPEadd(&s,&tmp);"); Print(" if ((c=DPEcmp(&s,".A."+n))==0) {fprintf(stderr,\"FAIL\\n\"); exit(1);}"); Print(" else if (c==1) break;"); Print(" k=n".s0."-k;"); Print(" }"); Print(" PROD_OPEN; draw".B."(k); PROD_SEP; draw".C."(n-k); PROD_CLOSE;"); end_if end_proc: cs::makedraw2[C,DPE,ProdPrd] := proc(B,C,A) local k0,k1,s0; begin k0:=cs::val[B]; k1:=cs::val[C]; if k0=cs::max[B] then s0:=(if k0=0 then "" else "-".k0 end_if); Print(" draw".B."(".k0."); PROD_SEP; draw".C."(n".s0.");") elif k1=cs::max[C] then s0:=(if k1=0 then "" else "-".k1 end_if); Print(" draw".B."(n".s0."); PROD_SEP; draw".C."(".k1.");") else Print(" DPEmulr(&s,count".A."(n),drand48());"); Print(" for (k=".cs::val[B].";;k++) {"); Print(" DPEmul(&tmp,count".B."(k),count".C."(n-k));"); Print(" DPEadd(&s,&tmp);"); Print(" if ((c=DPEcmp(&s,".A."+n))==0) {fprintf(stderr,\"FAIL\\n\"); exit(1);}"); Print(" else if (c==1) break;"); k0:=k0-cs::val[C]; s0:=(if k0=0 then "" elif k0>0 then "+".k0 else "-".(-k0) end_if); Print(" k=n".s0."-k;"); Print(" DPEmul(&tmp,count".B."(k),count".C."(n-k));"); Print(" DPEadd(&s,&tmp);"); Print(" if ((c=DPEcmp(&s,".A."+n))==0) {fprintf(stderr,\"FAIL\\n\"); exit(1);}"); Print(" else if (c==1) break;"); Print(" k=n".s0."-k;"); Print(" }"); Print(" draw".B."(k); PROD_SEP; draw".C."(n-k);"); end_if end_proc: cs::makedraw2[C,DPE,ProdSeq] := proc(B,C,A) local k0,k1,s0; begin k0:=cs::val[B]; k1:=cs::val[C]; if k0=cs::max[B] then s0:=(if k0=0 then "" else "-".k0 end_if); Print(" printf(\"Sequence(\"); draw".B."(".k0.");"); Print(" if ((n".s0.")>0) {putchar(','); draw".C."(n".s0.");}"); Print(" putchar(')');") elif k1=cs::max[C] then s0:=(if k1=0 then "" else "-".k1 end_if); Print(" printf(\"Sequence(\"); draw".B."(n".s0.");"); Print(" if ((".k1.")>0) {putchar(','); draw".C."(".k1.");}"); Print(" putchar(')');") else Print(" DPEmulr(&s,count".A."(n),drand48());"); Print(" for (k=".cs::val[B].";;k++) {"); Print(" DPEmul(&tmp,count".B."(k),count".C."(n-k));"); Print(" DPEadd(&s,&tmp);"); Print(" if ((c=DPEcmp(&s,".A."+n))==0) {fprintf(stderr,\"FAIL\\n\"); exit(1);}"); Print(" else if (c==1) break;"); k0:=k0-cs::val[C]; s0:=(if k0=0 then "" elif k0>0 then "+".k0 else "-".(-k0) end_if); Print(" k=n".s0."-k;"); Print(" DPEmul(&tmp,count".B."(k),count".C."(n-k));"); Print(" DPEadd(&s,&tmp);"); Print(" if ((c=DPEcmp(&s,".A."+n))==0) {fprintf(stderr,\"FAIL\\n\"); exit(1);}"); Print(" else if (c==1) break;"); Print(" k=n".s0."-k;"); Print(" }"); Print(" printf(\"Sequence(\"); draw".B."(k);"); Print(" if ((n-k)>0) {putchar(','); draw".C."(n-k);}"); Print(" putchar(')');"); end_if end_proc: cs::makedraw2[C,DPE,ProdPrdSeq] := proc(B,C,A) local k0,k1,s0; begin k0:=cs::val[B]; k1:=cs::val[C]; if k0=cs::max[B] then s0:=(if k0=0 then "" else "-".k0 end_if); Print(" draw".B."(".k0.");"); Print(" if ((n".s0.")>0) {putchar(','); draw".C."(n".s0.");}"); elif k1=cs::max[C] then s0:=(if k1=0 then "" else "-".k1 end_if); Print(" draw".B."(n".s0.");"); Print(" if ((".k1.")>0) {putchar(','); draw".C."(".k1.");}"); else Print(" DPEmulr(&s,count".A."(n),drand48());"); Print(" for (k=".cs::val[B].";;k++) {"); Print(" DPEmul(&tmp,count".B."(k),count".C."(n-k));"); Print(" DPEadd(&s,&tmp);"); Print(" if ((c=DPEcmp(&s,".A."+n))==0) {fprintf(stderr,\"FAIL\\n\"); exit(1);}"); Print(" else if (c==1) break;"); k0:=k0-cs::val[C]; s0:=(if k0=0 then "" elif k0>0 then "+".k0 else "-".(-k0) end_if); Print(" k=n".s0."-k;"); Print(" DPEmul(&tmp,count".B."(k),count".C."(n-k));"); Print(" DPEadd(&s,&tmp);"); Print(" if ((c=DPEcmp(&s,".A."+n))==0) {fprintf(stderr,\"FAIL\\n\"); exit(1);}"); Print(" else if (c==1) break;"); Print(" k=n".s0."-k;"); Print(" }"); Print(" draw".B."(k);"); Print(" if ((n-k)>0) {putchar(','); draw".C."(n-k);}"); end_if end_proc: # idem cs::makedraw2[C,DPE,ProdSeq] with Sequence --> Set # cs::makedraw2[C,DPE,ProdSet] := proc(B,C,A) local k0,k1,s0; begin k0:=cs::val[B]; k1:=cs::val[C]; if k0=cs::max[B] then s0:=(if k0=0 then "" else "-".k0 end_if); Print(" printf(\"Set(\"); draw".B."(".k0.");"); Print(" if ((n".s0.")>0) {putchar(','); draw".C."(n".s0.");}"); Print(" putchar(')');") elif k1=cs::max[C] then s0:=(if k1=0 then "" else "-".k1 end_if); Print(" printf(\"Set(\"); draw".B."(n".s0.");"); Print(" if ((".k1.")>0) {putchar(','); draw".C."(".k1.");}"); Print(" putchar(')');") else Print(" DPEmulr(&s,count".A."(n),drand48());"); Print(" for (k=".cs::val[B].";;k++) {"); Print(" DPEmul(&tmp,count".B."(k),count".C."(n-k));"); Print(" DPEadd(&s,&tmp);"); Print(" if ((c=DPEcmp(&s,".A."+n))==0) {fprintf(stderr,\"FAIL\\n\"); exit(1);}"); Print(" else if (c==1) break;"); k0:=k0-cs::val[C]; s0:=(if k0=0 then "" elif k0>0 then "+".k0 else "-".(-k0) end_if); Print(" k=n".s0."-k;"); Print(" DPEmul(&tmp,count".B."(k),count".C."(n-k));"); Print(" DPEadd(&s,&tmp);"); Print(" if ((c=DPEcmp(&s,".A."+n))==0) {fprintf(stderr,\"FAIL\\n\"); exit(1);}"); Print(" else if (c==1) break;"); Print(" k=n".s0."-k;"); Print(" }"); Print(" printf(\"Set(\"); draw".B."(k);"); Print(" if ((n-k)>0) {putchar(','); draw".C."(n-k);}"); Print(" putchar(')');"); end_if end_proc: cs::makedraw2[C,DPE,ProdPrdSet] := subsop(cs::makedraw2[C,DPE,ProdPrdSeq], 6=hold(cs::makedraw2[C,DPE,ProdPrdSet])): #----------------------------------------------------- Products => linear recurrences -----------------------------------------------------# cs::makecountprod[C,DPE]:=proc(sys,algo) begin cs::makecountprod[C,DPE,algo](sys) end_proc: cs::makecountprod[C,DPE,HOLONOMIC]:=proc(sys) local s,pr,npr,p,r,u,n,NT; begin print(Unquoted,"BE CAREFUL !"); print(Unquoted,"Verify the recurrences (in the C file) don't have negative coefficients."); ([pr,npr,s]):=split(sys,cs::isProd[HOLONOMIC],sys); cs::makecount[C,DPE](npr); for p in pr do Print("/* ".expr2text(p)." */"); # count function for A is the same as NT ? # NT:=cs::lequalcount[op(p,1)]; if type(NT)<>"_index" then Print("DPE* count".expr2text(op(p,1))."(int n) {"); Print(" return(count".NT."(n));"); Print("}\n"); else # otherwise # r:=cs::speciftorec(sys,op(p,1),u(n)); userinfo(1,"recurrence for ",p," is: ",expand(r)); cs::makecountprod1[C,DPE,HOLONOMIC](op(p),r,u(n)); Print("}\n"); cs::makeinitrec[C,DPE](op(p,1)); end_if; end_for; end_proc: cs::makecountprod1[C,DPE,HOLONOMIC] := proc(A,rhs,r,uofn) local u,n,R,ORDER,lroots,i,k,den,v,NZ; begin u:=op(uofn,0); n:=op(uofn,1); R:=cs::formatrec(r,u,n); ORDER:=nops(R)-2; R:=subs(R,n=n-ORDER); cs::LIM:=cs::LIM+ORDER; userinfo(1,"normalized coefficients are: ",R); # list of integer solutions of the equation coeff(u(n))=0 # lroots:=op(sharelib::iroots(R[nops(R)],n)); # we will have to compute initial values up to this limit # userinfo(1,"LIM,ORDER,ROOTS: ",cs::LIM, ORDER, lroots); cs::LIM:=max(cs::LIM, ORDER, lroots); R:=subs(R,n=k); # Printing of the recurrence # den:=expr2text(-R[nops(R)]); v:="".A."(k) = "; NZ:=FALSE; if R[1]<>0 then NZ:=TRUE; v:=v.expr2text(R[1])."/(".den.")"; end_if; for i from 0 to ORDER-1 do if R[i+2]<>0 then if NZ then v:=v." + "; end_if; v:=v.A."(k".expr2text(-ORDER+i).")*(".expr2text(R[i+2]).")/(".den.")"; NZ:=TRUE; end_if; end_for; Print("/* ".v." */"); # Print("/* ".expr2text(R)." */");# Print("DPE* count".A."(int n) { int k,nmax;"); Print(" if (NOINIT(".A."+n)) {"); Print(" nmax=0;"); Print(" while (INIT(".A."+nmax)) nmax++;"); Print(" for (k=nmax;k<=n;k++) {"); Print(" DPEset(".A."+".k.",(double)(".expr2text(R[1])."));"); for i from 0 to ORDER-1 do if R[i+2]<>0 then Print(" DPEcopy(&tmp,count".A."(k+(".expr2text(-ORDER+i).")));"); Print(" DPEmuli(&tmp,".expr2text(R[i+2]).");"); Print(" DPEadd(".A."+".k.",&tmp);"); end_if; end_for; Print(" DPEdivi(".A."+".k.",".expr2text(-R[nops(R)]).");"); Print(" }"); Print(" }"); Print(" return(".A."+n);"); end_proc: cs::makeinitrec[C,DPE]:=proc(A) begin Print("void initrec".A."(int n) {"); for i from cs::val[A] to cs::LIM do Print(" if (".i."<=n)"); Print(" DPEset(".A."+".i.",".expr2text(cs::tCount[cs::user_spec_NF][A](i)).".0);"); end_for; Print("}\n"); end_proc: #----------------------------------------------------- Products => lazy algorithm of Joris VdH -----------------------------------------------------# cs::makecountprod[C,DPE,LAZY]:=proc(sys) local p; begin cs::error(C,DPE,"lazy products in C with DPE: not yet implemented. Put cs::LAZY to FALSE"); end_proc: #----------------------------------------------------- Count and Draw Functions for C with GMP target=C, library=GMP -----------------------------------------------------# cs::init[C,GMP]:=proc(s,sys,fname,main) local A,eq; begin print(Unquoted,"CS warning: ".main." is taken as main non-terminal for C code."); cs::filehead(s,sys,fname,"GMP"); cs::printgmp(); cs::printwords(): Print("mpz_t s,tmp;"); for eq in sys do A:=expr2text(op(eq,1)); Print("mpz_t *".A.";"); Print("void count".A."();"); Print("void draw".A."();"); end_for; if cs::HOLONOMIC and cs::testcontextfreespecif(sys) then for A in select(sys,cs::isProd[HOLONOMIC],sys) do if type(cs::lequalcount[op(A,1)])="_index" then Print("void initrec".op(A,1)."();"); end_if; end_for; end_if; Print("\n"); cs::makemain[C,GMP](sys,main); end_proc: cs::reset[C,GMP]:=cs::filetail: cs::makemain[C,GMP]:=proc(sys,main) begin Print("main(int argc,char *argv[]) {"); Print(" int i,n,k;"); Print(" if (argc<=1) {"); Print(" printf(\"Usage: a.out n [k seed] to generate k objects of size n\\n\");"); Print(" exit(1);}"); Print(" n=atoi(argv[1]);"); Print(" k=(argc>=3) ? atoi(argv[2]) : 1;"); Print(" srand48((argc>=4) ? atoi(argv[3]) : getpid());"); Print(" mpz_init(tmp);"); Print(" mpz_init(s);"); for A in sys do A:=op(A,1); # test if count function for A will not be deduced from another one # if type(cs::lequalcount[A])="_index" then Print(" ".A."=(mpz_t *)malloc((n+1)*sizeof(mpz_t));"); Print(" for (i=0;i<".cs::val[A]." && i<=n;i++) mpz_init_set_si(".A."[i],0);"); Print(" for (;i<=n;i++) mpz_init_set_si(".A."[i],-1);"); end_if; end_for; # count functions deduced from other ones # for A in cs::listofindices(cs::lequalcount) do Print(" ".A."=".expr2text(cs::lequalcount[A]).";"); end_for; if cs::HOLONOMIC and cs::testcontextfreespecif(sys) then for A in select(sys,cs::isProd[HOLONOMIC],sys) do if type(cs::lequalcount[op(A,1)])="_index" then Print(" initrec".op(A,1)."(n);"); end_if; end_for; end_if; A:=main; Print(" count".A."(n);"); Print(" if (mpz_cmp_ui(".A."[n],0)<=0) {printf(\"there is no such structure of this size\\n\"); exit(1);}"); Print(" if (k==0) {mpz_out_str (stdout, 10, ".A."[n]); puts(\"\");}"); Print(" for (i=1;i<=k;i++) { draw".A."(n); putchar('\\n'); }"); Print(" mpz_clear(tmp);"); Print(" mpz_clear(s);"); Print("}\n"); end_proc: cs::makecount[C,GMP] := proc(sys) local p; begin for p in sys do cs::makecount1[C,GMP](op(p)); end_for; end_proc: cs::makecount1[C,GMP] := proc(A,rhs) local f,i,NT; begin Print("/* ".expr2text(A=rhs)." */"); Print("void count".A."(int n) { int k;"); # count function for A is the same as NT ? # NT:=cs::lequalcount[A]; if type(NT)<>"_index" then Print(" count".expr2text(NT)."(n);"); else # otherwise # Print(" if (NOINIT(".A."[n])==0) {"); f:=cs::makecount2[C,GMP,cs::typeofoper(rhs)]; if not contains({DOM_FUN,DOM_PROC},type(f)) then error("no C counting rule for constructor ".op(rhs,0)) end_if; f(op(rhs),A); Print(" }"); end_if; Print("}\n"); end_proc: #..................................................................... DRAW DRAW DRAW DRAW DRAW DRAW DRAW DRAW DRAW DRAW DRAW DRAW DRAW DRAW 12/02/99 pour draw comparer a 0 et faire des soustractions .....................................................................# cs::makedraw[C,GMP] := proc(sys) local p; begin for p in sys do cs::makedraw1[C,GMP](op(p)); end_for; end_proc: cs::makedraw1[C,GMP] := proc(A,rhs) local f,i; begin Print("/* ".expr2text(A=rhs)." */"); Print("void draw".A."(int n) { int k,c;"); f:=cs::makedraw2[C,GMP,op(rhs,0)]; if not contains({DOM_FUN,DOM_PROC},type(f)) then error("no C drawing rule for constructor ".op(rhs,0)) end_if; f(op(rhs),A); Print("}\n"); end_proc: # Epsilon # cs::makecount2[C,GMP,Epsilon] := proc() begin Print(" if (n==0)"); Print(" mpz_set_si(".args(2)."[n],1);"); Print(" else"); Print(" mpz_set_si(".args(2)."[n],0);"); end_proc: cs::makedraw2[C,GMP,Epsilon]:=proc() begin Print(" EPSILON;") end_proc: # Atom # cs::makecount2[C,GMP,Atom] := proc() begin Print(" if (n==1)"); Print(" mpz_set_si(".args(2)."[n],1);"); Print(" else"); Print(" mpz_set_si(".args(2)."[n],0);"); end_proc: cs::makedraw2[C,GMP,Atom]:=proc(A) begin Print(" printf(\"".A."\");") end_proc: # Union # cs::makecount2[C,GMP,Union] := proc() local i,k; begin k:=args(0); Print(" count".args(1)."(n);"); Print(" mpz_set(".args(k)."[n],".args(1)."[n]);"); for i from 2 to k-1 do Print(" count".args(i)."(n);"); Print(" mpz_add(".args(k)."[n],".args(k)."[n],".args(i)."[n]);"); end_for end_proc: cs::makedraw2[C,GMP,Union] := proc() local i,k; begin k:=args(0); Print(" count".args(k)."(n);"); Print(" mpz_random(s,mpz_size(".args(k)."[n]));"); Print(" mpz_mod(s,s,".args(k)."[n]);"); for i from 1 to k-1 do if i>1 then Print(" else {") end_if; Print(" count".args(i)."(n);"); Print(" mpz_add(s,s,".args(i)."[n]);"); Print(" if (mpz_cmp(s,".args(k)."[n])>=0)"); Print(" { draw".args(i)."(n); return; }"); end_for; for i from 1 to k-2 do Print(" }"); end_for; end_proc: # B=Int(A) ==> B(n) = A(n)/n # cs::makecount2[C,GMP,Int] := proc(A,B) begin Print(" count".A."(n);"); Print(" mpz_fdiv_q_ui(".B."[n],".A."[n],n);"); end_proc: cs::makedraw2[C,GMP,Int] := proc(A,B) begin Print(" draw".A."(n);"); end_proc: # B=Theta(A) ==> B(n) = n*A(n) # cs::makecount2[C,GMP,Theta] := proc(A,B) begin Print(" count".A."(n);"); Print(" mpz_mul_ui(".B."[n],".A."[n],n);"); end_proc: cs::makedraw2[C,GMP,Theta] := proc(A,B) begin Print(" draw".A."(n);"); end_proc: # A = Delta(B,u) # cs::makecount2[C,GMP,Delta]:=proc(B,u,A) begin if u=fun(1) then Print(" count".B."(n);"); Print(" mpz_set(".A."[n],".B."[n]);"); Print(" for (k=2;k<=n;k++)"); Print(" if (n%k==0) {"); Print(" count".B."(n/k);"); Print(" mpz_add(".A."[n],".A."[n],".B."[n/k]);"); Print(" }"); else error("not yet implemented for u(n)<>1") end_if end_proc: cs::makedraw2[C,GMP,Delta]:=proc(B,u,A) begin if u=fun(1) then Print(" count".A."(n);"); Print(" mpz_random(s,mpz_size(".A."[n]));"); Print(" mpz_mod(s,s,".A."[n]);"); Print(" for (k=1;k<=n;k++)"); Print(" if (n%k==0) {"); Print(" count".B."(n/k);"); Print(" mpz_add(s,s,".B."[n/k]);"); Print(" if (mpz_cmp(s,".A."[n])>=0)"); Print(" { draw".B."(n/k); if (k>1) printf(\"$%u\",k); return; }"); Print(" }"); Print(" printf(\"Error: you should not get there\\n\"); exit(1);"); else error("not yet implemented for u(n)<>1") end_if end_proc: #----------------------------------------------------- Naive Products -----------------------------------------------------# cs::makecount2[C,GMP,Prod] := proc(B,C,A) local k0,s0; begin if cs::val[C]=cs::max[C] then # exchange B and C, as counting is symetrical # B; B:=C; C:=%2 end_if; k0:=cs::val[B]; s0:=(if k0=0 then "" else "-".k0 end_if); Print(" count".B."(".k0."); count".C."(n".s0.");"); Print(" mpz_mul(".A."[n],".B."[".k0."],".C."[n".s0."]);"); if k0<>cs::max[B] then s0:=(if cs::val[C]=0 then "" else "-".cs::val[C] end_if); Print(" for (k=".(k0+1).";k<=n".s0.";k++) {"); Print(" count".B."(".k."); count".C."(n-k);"); Print(" mpz_mul(tmp,".B."[k],".C."[n-k]);"); Print(" mpz_add(".A."[n],".A."[n],tmp);"); Print(" }"); end_if end_proc: cs::makedraw2[C,GMP,Prod] := proc(B,C,A) local k0,k1,s0; begin k0:=cs::val[B]; k1:=cs::val[C]; if k0=cs::max[B] then s0:=(if k0=0 then "" else "-".k0 end_if); Print(" PROD_OPEN; draw".B."(".k0."); PROD_SEP; draw".C."(n".s0."); PROD_CLOSE;") elif k1=cs::max[C] then s0:=(if k1=0 then "" else "-".k1 end_if); Print(" PROD_OPEN; draw".B."(n".s0."); PROD_SEP; draw".C."(".k1."); PROD_CLOSE;") else Print(" count".A."(n);"); Print(" mpz_random(s,mpz_size(".A."[n]));"); Print(" mpz_mod(s,s,".A."[n]);"); Print(" for (k=".cs::val[B].";;k++) {"); Print(" count".B."(".k."); count".C."(n-k);"); Print(" mpz_mul(tmp,".B."[".k."],".C."[n-k]);"); Print(" mpz_add(s,s,tmp);"); Print(" if (mpz_cmp(s,".A."[n])>=0) break;"); k0:=k0-cs::val[C]; #boustrophedon# s0:=(if k0=0 then "" elif k0>0 then "+".k0 else "-".(-k0) end_if); Print(" k=n".s0."-k;"); Print(" count".B."(".k."); count".C."(n-k);"); Print(" mpz_mul(tmp,".B."[".k."],".C."[n-k]);"); Print(" mpz_add(s,s,tmp);"); Print(" if (mpz_cmp(s,".A."[n])>=0) break;"); Print(" k=n".s0."-k;"); Print(" }"); Print(" PROD_OPEN; draw".B."(k); PROD_SEP; draw".C."(n-k); PROD_CLOSE;"); end_if end_proc: cs::makedraw2[C,GMP,ProdPrd] := proc(B,C,A) local k0,k1,s0; begin k0:=cs::val[B]; k1:=cs::val[C]; if k0=cs::max[B] then s0:=(if k0=0 then "" else "-".k0 end_if); Print(" PROD_OPEN; draw".B."(".k0."); PROD_SEP; draw".C."(n".s0."); PROD_CLOSE;") elif k1=cs::max[C] then s0:=(if k1=0 then "" else "-".k1 end_if); Print(" PROD_OPEN; draw".B."(n".s0."); PROD_SEP; draw".C."(".k1."); PROD_CLOSE;") else Print(" count".A."(n);"); Print(" mpz_random(s,mpz_size(".A."[n]));"); Print(" mpz_mod(s,s,".A."[n]);"); Print(" for (k=".cs::val[B].";;k++) {"); Print(" count".B."(".k."); count".C."(n-k);"); Print(" mpz_mul(tmp,".B."[".k."],".C."[n-k]);"); Print(" mpz_add(s,s,tmp);"); Print(" if (mpz_cmp(s,".A."[n])>=0) break;"); k0:=k0-cs::val[C]; #boustrophedon# s0:=(if k0=0 then "" elif k0>0 then "+".k0 else "-".(-k0) end_if); Print(" k=n".s0."-k;"); Print(" count".B."(".k."); count".C."(n-k);"); Print(" mpz_mul(tmp,".B."[".k."],".C."[n-k]);"); Print(" mpz_add(s,s,tmp);"); Print(" if (mpz_cmp(s,".A."[n])>=0) break;"); Print(" k=n".s0."-k;"); Print(" }"); Print(" PROD_OPEN; draw".B."(k); PROD_SEP; draw".C."(n-k); PROD_CLOSE;"); end_if end_proc: cs::makedraw2[C,GMP,ProdSeq] := proc(B,C,A) local k0,k1,s0; begin k0:=cs::val[B]; k1:=cs::val[C]; if k0=cs::max[B] then s0:=(if k0=0 then "" else "-".k0 end_if); Print(" printf(\"Sequence(\"); draw".B."(".k0.");"); Print(" if ((n".s0.")>0) {putchar(','); draw".C."(n".s0.");}"); Print(" putchar(')');") elif k1=cs::max[C] then s0:=(if k1=0 then "" else "-".k1 end_if); Print(" printf(\"Sequence(\"); draw".B."(n".s0.");"); Print(" if ((".k1.")>0) {putchar(','); draw".C."(".k1.");}"); Print(" putchar(')');") else Print(" count".A."(n);"); Print(" mpz_random(s,mpz_size(".A."[n]));"); Print(" mpz_mod(s,s,".A."[n]);"); Print(" count".A."(n);"); Print(" for (k=".cs::val[B].";;k++) {"); Print(" count".B."(".k."); count".C."(n-k);"); Print(" mpz_mul(tmp,".B."[".k."],".C."[n-k]);"); Print(" mpz_add(s,s,tmp);"); Print(" if ((c=mpz_cmp(s,".A."[n]))>=0) break;"); k0:=k0-cs::val[C]; s0:=(if k0=0 then "" elif k0>0 then "+".k0 else "-".(-k0) end_if); Print(" k=n".s0."-k;"); Print(" count".B."(".k."); count".C."(n-k);"); Print(" mpz_mul(tmp,".B."[".k."],".C."[n-k]);"); Print(" mpz_add(s,s,tmp);"); Print(" if ((c=mpz_cmp(s,".A."[n]))>=0) break;"); Print(" k=n".s0."-k;"); Print(" }"); Print(" printf(\"Sequence(\"); draw".B."(k);"); Print(" if ((n-k)>0) {putchar(','); draw".C."(n-k);}"); Print(" putchar(')');"); end_if end_proc: cs::makedraw2[C,GMP,ProdPrdSeq] := proc(B,C,A) local k0,k1,s0; begin k0:=cs::val[B]; k1:=cs::val[C]; if k0=cs::max[B] then s0:=(if k0=0 then "" else "-".k0 end_if); Print(" draw".B."(".k0.");"); Print(" if ((n".s0.")>0) {putchar(','); draw".C."(n".s0.");}") elif k1=cs::max[C] then s0:=(if k1=0 then "" else "-".k1 end_if); Print(" draw".B."(n".s0.");"); Print(" if ((".k1.")>0) {putchar(','); draw".C."(".k1.");}") else Print(" count".A."(n);"); Print(" mpz_random(s,mpz_size(".A."[n]));"); Print(" mpz_mod(s,s,".A."[n]);"); Print(" count".A."(n);"); Print(" for (k=".cs::val[B].";;k++) {"); Print(" count".B."(".k."); count".C."(n-k);"); Print(" mpz_mul(tmp,".B."[".k."],".C."[n-k]);"); Print(" mpz_add(s,s,tmp);"); Print(" if (mpz_cmp(s,".A."[n])>=0) break;"); k0:=k0-cs::val[C]; s0:=(if k0=0 then "" elif k0>0 then "+".k0 else "-".(-k0) end_if); Print(" k=n".s0."-k;"); Print(" count".B."(".k."); count".C."(n-k);"); Print(" mpz_mul(tmp,".B."[".k."],".C."[n-k]);"); Print(" mpz_add(s,s,tmp);"); Print(" if (mpz_cmp(s,".A."[n])>=0) break;"); Print(" k=n".s0."-k;"); Print(" }"); Print(" draw".B."(k);"); Print(" if ((n-k)>0) {putchar(','); draw".C."(n-k);}"); end_if end_proc: # idem cs::makedraw2[C,GMP,ProdSeq] with Sequence --> Set # cs::makedraw2[C,GMP,ProdSet] := proc(B,C,A) local k0,k1,s0; begin k0:=cs::val[B]; k1:=cs::val[C]; if k0=cs::max[B] then s0:=(if k0=0 then "" else "-".k0 end_if); Print(" printf(\"Set(\"); draw".B."(".k0.");"); Print(" if ((n".s0.")>0) {putchar(','); draw".C."(n".s0.");}"); Print(" putchar(')');") elif k1=cs::max[C] then s0:=(if k1=0 then "" else "-".k1 end_if); Print(" printf(\"Set(\"); draw".B."(n".s0.");"); Print(" if ((".k1.")>0) {putchar(','); draw".C."(".k1.");}"); Print(" putchar(')');") else Print(" count".A."(n);"); Print(" mpz_random(s,mpz_size(".A."[n]));"); Print(" mpz_mod(s,s,".A."[n]);"); Print(" count".A."(n);"); Print(" for (k=".cs::val[B].";;k++) {"); Print(" count".B."(".k."); count".C."(n-k);"); Print(" mpz_mul(tmp,".B."[".k."],".C."[n-k]);"); Print(" mpz_add(s,s,tmp);"); Print(" if (mpz_cmp(s,".A."[n])>=0) break;"); k0:=k0-cs::val[C]; s0:=(if k0=0 then "" elif k0>0 then "+".k0 else "-".(-k0) end_if); Print(" k=n".s0."-k;"); Print(" count".B."(".k."); count".C."(n-k);"); Print(" mpz_mul(tmp,".B."[".k."],".C."[n-k]);"); Print(" mpz_add(s,s,tmp);"); Print(" if (mpz_cmp(s,".A."[n])>=0) break;"); Print(" k=n".s0."-k;"); Print(" }"); Print(" printf(\"Set(\"); draw".B."(k);"); Print(" if ((n-k)>0) {putchar(','); draw".C."(n-k);}"); Print(" putchar(')');"); end_if end_proc: cs::makedraw2[C,GMP,ProdPrdSet]:= subsop(cs::makedraw2[C,GMP,ProdPrdSeq], 6=hold(cs::makedraw2[C,GMP,ProdPrdSet])): #----------------------------------------------------- Products => linear recurrences -----------------------------------------------------# cs::makecountprod[C,GMP]:=proc(sys,algo) begin cs::makecountprod[C,GMP,algo](sys) end_proc: cs::makecountprod[C,GMP,HOLONOMIC]:=proc(sys) begin cs::error(C,GMP,"recurrences in C with GMP: not yet implemented."); end_proc: #----------------------------------------------------- Products => lazy algorithm of Joris VdH -----------------------------------------------------# cs::makecountprod[C,GMP,LAZY]:=proc(sys) local p; begin cs::error(C,GMP,"lazy products in C with GMP: not yet implemented. Put cs::LAZY to FALSE"); end_proc: #----------------------------------------------------- Count and Draw Functions for C Modules with DPE target=MODULE, library=DPE -----------------------------------------------------# cs::init[MODULE,DPE]:=proc(s,sys,fname,main) local A,eq,l,i,ss; begin cs::filehead(s,sys,fname,"DPE"); Print("/* This is a C Dynamic Module for MuPAD */\n"); cs::printdpe(): Print("DPE s,tmp;"); for eq in sys do A:=expr2text(op(eq,1)); Print("DPE *".A.",*count".A."(int);"); Print("MTcell draw".A."(int);"); end_for; if cs::HOLONOMIC and cs::testcontextfreespecif(sys) then for A in select(sys,cs::isProd[HOLONOMIC],sys) do if type(cs::lequalcount[op(A,1)])="_index" then Print("void initrec".op(A,1)."();"); end_if; end_for; end_if; Print("\nint Nmax = -1;\n"); l:=[op(map(sys,op,1))]; for i from 0 to nops(l)-1 do Print("#define nt".expr2text(op(l,i+1))." ".i); end_for; ss:="\nDPE *(*count[])(int) = { "; map(l,fun(hold(count).expr2text(args(1)))); ss:=ss.expr2text(op(%))." };"; Print(ss); ss:="MTcell (*draw[])(int) = { "; map(l,fun(hold(draw).expr2text(args(1)))); ss:=ss.expr2text(op(%))." };"; Print(ss."\n"); cs::makemain[MODULE,DPE](sys); end_proc: cs::reset[MODULE,DPE]:=cs::filetail: cs::makemain[MODULE,DPE]:=proc(sys) begin cs::makeDPE2list[MODULE,DPE](); cs::makeinit[MODULE,DPE](sys); cs::makereinit[MODULE,DPE](sys); cs::makecountint[MODULE,DPE](sys); cs::makeMFcount[MODULE,DPE](sys); cs::makeMFdraw[MODULE,DPE](sys); end_proc: cs::makeDPE2list[MODULE,DPE]:=proc() begin Print("MTcell DPE2list(DPE *x) {"); Print("//MFeval(MFtext2expr(\"old_DIGITS:=DIGITS\"));"); Print("//MFeval(MFtext2expr(\"DIGITS:=16\"));"); Print(" MTcell e = MFlong( x->e );"); Print(" MTcell l = MFnewList(2);"); Print(" MTcell p = MFcall(\"_power\",2,MFlong(10),MFcopy(e));"); Print(" MFsetList(&l,0,MFcall(\"_mult\",2,MFdouble( x->lo ),MFcopy(p)));"); Print(" MFsetList(&l,1,MFcall(\"_mult\",2,MFdouble( x->hi ),MFcopy(p))); "); Print(" MFsig(l);"); Print("//MFeval(MFtext2expr(\"DIGITS:=old_DIGITS\"));"); Print("//MFeval(MFtext2expr(\"unassign(old_DIGITS)\"));"); Print(" return(l);"); Print("}\n"); end_proc: cs::makeinit[MODULE,DPE]:=proc(sys) local valmax,A; begin # cs::valuations has been called before (in cs::compile1) # valmax:=max(op(cs::listofvalues(cs::val))); Print("void init(int n) {"); Print(" int i;"); Print(" srand48(getpid());"); Print(" DPEinit();"); Print("//MFprintf(\"init n=%d\\n\",n);"); Print(" if (n<".valmax.") n=".valmax.";"); for A in sys do A:=op(A,1); # test if count function for A will not be deduced from another one # if type(cs::lequalcount[A])="_index" then Print(" ".A."=(DPE*)malloc((n+1)*sizeof(DPE));"); Print(" for (i=0;i<".cs::val[A].";i++) DPEset(".A."+i,0.0);"); Print(" for (;i<=n;i++) ".A."[i].lo=Nan;"); end_if; end_for; # count functions deduced from other ones # for A in cs::listofindices(cs::lequalcount) do Print(" ".A."=".expr2text(cs::lequalcount[A]).";"); end_for; Print(" Nmax=n;"); if cs::HOLONOMIC and cs::testcontextfreespecif(sys) then for A in select(sys,cs::isProd[HOLONOMIC],sys) do if type(cs::lequalcount[op(A,1)])="_index" then Print(" initrec".op(A,1)."();"); end_if; end_for; end_if; Print("}\n"); end_proc: cs::makereinit[MODULE,DPE]:=proc(sys) local A; begin Print("void reinit(int n) {"); Print(" int i;"); Print("//MFprintf(\"reinit n=%d\\n\",n);"); for A in sys do A:=op(A,1); # test if count function for A will not be deduced from another one # if type(cs::lequalcount[A])="_index" then Print(" ".A."=(DPE*)realloc(".A.", (n+1)*sizeof(DPE));"); Print(" for (i=Nmax+1;i<=n;i++) ".A."[i].lo=Nan;"); end_if; end_for; # count functions deduced from other ones # for A in cs::listofindices(cs::lequalcount) do Print(" ".A."=".expr2text(cs::lequalcount[A]).";"); end_for; Print(" Nmax=n;"); Print("}\n"); end_proc: cs::makecountint[MODULE,DPE]:=proc(sys) begin Print("void countint(int i, int n) {"); Print(" if (n<0) "); Print(" MFerror(\"Bad size\\n\");"); Print("//MFprintf(\"Nmax=%d n=%d\\n\",Nmax,n);"); Print(" if (Nmax==-1) init(n);"); Print(" if (Nmaxlo<=0.0) /* se bloque quand on fait B[n].lo !! */"); Print(" MFerror(\"There is no such structure of this size\\n\");"); Print(" MFreturn( MFeval( draw[i](n) ) ); "); Print("} MFEND \n"); end_proc: cs::makecount[MODULE,DPE] := proc(sys) local p; begin for p in sys do cs::makecount1[C,DPE](op(p)); end_for; end_proc: cs::makedraw[MODULE,DPE] := proc(sys) local p; begin for p in sys do cs::makedraw1[MODULE,DPE](op(p)); end_for; end_proc: cs::makedraw1[MODULE,DPE] := proc(A,rhs) local f,i; begin Print("/* ".expr2text(A=rhs)." */"); Print("MTcell draw".A."(int n) { int k,c;"); f:=cs::makedraw2[MODULE,DPE,op(rhs,0)]; if not contains({DOM_FUN,DOM_PROC},type(f)) then error("no C drawing rule for constructor ".op(rhs,0)) end_if; f(op(rhs),A); Print("} \n"); end_proc: # Epsilon # cs::makedraw2[MODULE,DPE,Epsilon]:=proc() begin Print(" return( MFident(\"Epsilon\") );") end_proc: # Atom # cs::makedraw2[MODULE,DPE,Atom]:=proc(A) begin Print(" return( MFident(\"".A."\") );") end_proc: # Union # cs::makedraw2[MODULE,DPE,Union] := proc() local i,k; begin k:=args(0); Print(" DPEmulr(&s,count".args(k)."(n),drand48());"); for i from 1 to k-1 do if i>1 then Print(" else {") end_if; Print(" DPEcopy(&tmp,count".args(i)."(n));"); Print(" DPEadd(&s,&tmp);"); Print(" if ((c=DPEcmp(&s,".args(k)."+n))==0) {MFerror(\"FAIL\\n\");}"); Print(" else if (c==1) { return( draw".args(i)."(n) ); }"); end_for; Print(" else {MFerror(\"FAIL\\n\");}"); Print(" "._concat("}"$(k-2))) end_proc: # B=Int(A) ==> B(n) = A(n)/n # cs::makedraw2[MODULE,DPE,Int] := proc(A,B) begin Print(" return( draw".A."(n) );"); end_proc: # B=Theta(A) ==> B(n) = n*A(n) # cs::makedraw2[MODULE,DPE,Theta] := proc(A,B) begin Print(" return( draw".A."(n) );"); end_proc: # A = Delta(B,u) # cs::makedraw2[MODULE,DPE,Delta]:=proc(B,u,A) begin Print(" MTcell res;"); if u=fun(1) then Print(" DPEmulr(&s,count".A."(n),drand48());"); Print(" for (k=1;k<=n;k++)"); Print(" if (n%k==0) {"); Print(" DPEcopy(&tmp,count".B."(n/k));"); Print(" DPEadd(&s,&tmp);"); Print(" if ((c=DPEcmp(&s,".A."+n))==0) {MFerror(\"FAIL\\n\");}"); Print(" else if (c==1) {"); Print(" res=MFnewExpr(3, MFident(\"_seqgen\"), draw".B."(n/k), MFint(k) );"); Print(" MFsig( res );"); Print(" return( res );"); Print(" }"); Print(" }"); Print(" MFerror(\"you should not get there\\n\");"); else MFerror("not yet implemented for u(n)<>1") end_if end_proc: #----------------------------------------------------- Naive Products -----------------------------------------------------# cs::makedraw2[MODULE,DPE,Prod] := proc(B,C,A) local k0,k1,s0; begin Print(" MTcell res;"); k0:=cs::val[B]; k1:=cs::val[C]; if k0=cs::max[B] then s0:=(if k0=0 then "" else "-".k0 end_if); Print(" res=MFnewExpr(3, MFident(\"Prod\"), draw".B."(".k0."), draw".C."(n".s0.") );"); elif k1=cs::max[C] then s0:=(if k1=0 then "" else "-".k1 end_if); Print(" res=MFnewExpr(3, MFident(\"Prod\"), draw".B."(n".s0."), draw".C."(".k1.") );"); else Print(" DPEmulr(&s,count".A."(n),drand48());"); Print(" for (k=".cs::val[B].";;k++) {"); Print(" DPEmul(&tmp,count".B."(k),count".C."(n-k));"); Print(" DPEadd(&s,&tmp);"); Print(" if ((c=DPEcmp(&s,".A."+n))==0) {MFerror(\"FAIL\\n\");}"); Print(" else if (c==1) break;"); k0:=k0-cs::val[C]; s0:=(if k0=0 then "" elif k0>0 then "+".k0 else "-".(-k0) end_if); Print(" k=n".s0."-k;"); Print(" DPEmul(&tmp,count".B."(k),count".C."(n-k));"); Print(" DPEadd(&s,&tmp);"); Print(" if ((c=DPEcmp(&s,".A."+n))==0) {MFerror(\"FAIL\\n\");}"); Print(" else if (c==1) break;"); Print(" k=n".s0."-k;"); Print(" }"); Print(" res=MFnewExpr(3, MFident(\"Prod\"), draw".B."(k), draw".C."(n-k) );"); end_if; Print(" MFsig( res );"); Print(" return( res );"); end_proc: cs::makedraw2[MODULE,DPE,ProdPrd] := proc(B,C,A) local k0,k1,s0; begin Print(" MTcell res;"); k0:=cs::val[B]; k1:=cs::val[C]; if k0=cs::max[B] then s0:=(if k0=0 then "" else "-".k0 end_if); Print(" res=MFnewExprSeq(2, draw".B."(".k0."), draw".C."(n".s0.") );"); elif k1=cs::max[C] then s0:=(if k1=0 then "" else "-".k1 end_if); Print(" res=MFnewExprSeq(2, draw".B."(n".s0."), draw".C."(".k1.") );"); else Print(" DPEmulr(&s,count".A."(n),drand48());"); Print(" for (k=".cs::val[B].";;k++) {"); Print(" DPEmul(&tmp,count".B."(k),count".C."(n-k));"); Print(" DPEadd(&s,&tmp);"); Print(" if ((c=DPEcmp(&s,".A."+n))==0) {MFerror(\"FAIL\\n\");}"); Print(" else if (c==1) break;"); k0:=k0-cs::val[C]; s0:=(if k0=0 then "" elif k0>0 then "+".k0 else "-".(-k0) end_if); Print(" k=n".s0."-k;"); Print(" DPEmul(&tmp,count".B."(k),count".C."(n-k));"); Print(" DPEadd(&s,&tmp);"); Print(" if ((c=DPEcmp(&s,".A."+n))==0) {MFerror(\"FAIL\\n\");}"); Print(" else if (c==1) break;"); Print(" k=n".s0."-k;"); Print(" }"); Print(" res=MFnewExprSeq(2, draw".B."(k), draw".C."(n-k) );"); end_if; Print(" MFsig( res );"); Print(" return( res );"); end_proc: cs::makedraw2[MODULE,DPE,ProdSeq] := proc(B,C,A) local k0,k1,s0; begin Print(" MTcell res;"); k0:=cs::val[B]; k1:=cs::val[C]; if k0=cs::max[B] then s0:=(if k0=0 then "" else "-".k0 end_if); Print(" if ((n".s0.")>0)"); Print(" res=MFnewExpr(3, MFident(\"Sequence\"), draw".B."(".k0."), draw".C."(n".s0.") );"); Print(" else"); Print(" res=MFnewExpr(2, MFident(\"Sequence\"), draw".B."(".k0.") );"); elif k1=cs::max[C] then s0:=(if k1=0 then "" else "-".k1 end_if); Print(" if ((".k1.")>0)"); Print(" res=MFnewExpr(3, MFident(\"Sequence\"), draw".B."(n".s0."), draw".C."(".k1.") );"); Print(" else"); Print(" res=MFnewExpr(2, MFident(\"Sequence\"), draw".B."(n".s0.") );"); else Print(" DPEmulr(&s,count".A."(n),drand48());"); Print(" for (k=".cs::val[B].";;k++) {"); Print(" DPEmul(&tmp,count".B."(k),count".C."(n-k));"); Print(" DPEadd(&s,&tmp);"); Print(" if ((c=DPEcmp(&s,".A."+n))==0) {MFerror(\"FAIL\\n\");}"); Print(" else if (c==1) break;"); k0:=k0-cs::val[C]; s0:=(if k0=0 then "" elif k0>0 then "+".k0 else "-".(-k0) end_if); Print(" k=n".s0."-k;"); Print(" DPEmul(&tmp,count".B."(k),count".C."(n-k));"); Print(" DPEadd(&s,&tmp);"); Print(" if ((c=DPEcmp(&s,".A."+n))==0) {MFerror(\"FAIL\\n\");}"); Print(" else if (c==1) break;"); Print(" k=n".s0."-k;"); Print(" }"); Print(" if ((n-k)>0)"); Print(" res=MFnewExpr(3, MFident(\"Sequence\"), draw".B."(k), draw".C."(n-k) );"); Print(" else"); Print(" res=MFnewExpr(2, MFident(\"Sequence\"), draw".B."(k) );"); end_if; Print(" MFsig( res );"); Print(" return( res );"); end_proc: cs::makedraw2[MODULE,DPE,ProdPrdSeq] := proc(B,C,A) local k0,k1,s0; begin Print(" MTcell res;"); k0:=cs::val[B]; k1:=cs::val[C]; if k0=cs::max[B] then s0:=(if k0=0 then "" else "-".k0 end_if); Print(" if ((n".s0.")>0)"); Print(" res=MFnewExprSeq(2, draw".B."(".k0."), draw".C."(n".s0.") );"); Print(" else"); Print(" res=draw".B."(".k0.");"); elif k1=cs::max[C] then s0:=(if k1=0 then "" else "-".k1 end_if); Print(" if ((".k1.")>0)"); Print(" res=MFnewExprSeq(2, draw".B."(n".s0."), draw".C."(".k1.") );"); Print(" else"); Print(" res=draw".B."(n".s0.");"); else Print(" DPEmulr(&s,count".A."(n),drand48());"); Print(" for (k=".cs::val[B].";;k++) {"); Print(" DPEmul(&tmp,count".B."(k),count".C."(n-k));"); Print(" DPEadd(&s,&tmp);"); Print(" if ((c=DPEcmp(&s,".A."+n))==0) {MFerror(\"FAIL\\n\");}"); Print(" else if (c==1) break;"); k0:=k0-cs::val[C]; s0:=(if k0=0 then "" elif k0>0 then "+".k0 else "-".(-k0) end_if); Print(" k=n".s0."-k;"); Print(" DPEmul(&tmp,count".B."(k),count".C."(n-k));"); Print(" DPEadd(&s,&tmp);"); Print(" if ((c=DPEcmp(&s,".A."+n))==0) {MFerror(\"FAIL\\n\");}"); Print(" else if (c==1) break;"); Print(" k=n".s0."-k;"); Print(" }"); Print(" if ((n-k)>0)"); Print(" res=MFnewExprSeq(2, draw".B."(k), draw".C."(n-k) );"); Print(" else"); Print(" res=draw".B."(k);"); end_if; Print(" MFsig( res );"); Print(" return( res );"); end_proc: # idem cs::makedraw2[MODULE,DPE,ProdSeq] with Sequence --> Set # cs::makedraw2[MODULE,DPE,ProdSet] := proc(B,C,A) local k0,k1,s0; begin Print(" MTcell res;"); k0:=cs::val[B]; k1:=cs::val[C]; if k0=cs::max[B] then s0:=(if k0=0 then "" else "-".k0 end_if); Print(" if ((n".s0.")>0)"); Print(" res=MFnewExpr(3, MFident(\"Set\"), draw".B."(".k0."), draw".C."(n".s0.") );"); Print(" else"); Print(" res=MFnewExpr(2, MFident(\"Set\"), draw".B."(".k0.") );"); elif k1=cs::max[C] then s0:=(if k1=0 then "" else "-".k1 end_if); Print(" if ((".k1.")>0)"); Print(" res=MFnewExpr(3, MFident(\"Set\"), draw".B."(n".s0."), draw".C."(".k1.") );"); Print(" else"); Print(" res=MFnewExpr(2, MFident(\"Set\"), draw".B."(n".s0.") );"); else Print(" DPEmulr(&s,count".A."(n),drand48());"); Print(" for (k=".cs::val[B].";;k++) {"); Print(" DPEmul(&tmp,count".B."(k),count".C."(n-k));"); Print(" DPEadd(&s,&tmp);"); Print(" if ((c=DPEcmp(&s,".A."+n))==0) {MFerror(\"FAIL\\n\");}"); Print(" else if (c==1) break;"); k0:=k0-cs::val[C]; s0:=(if k0=0 then "" elif k0>0 then "+".k0 else "-".(-k0) end_if); Print(" k=n".s0."-k;"); Print(" DPEmul(&tmp,count".B."(k),count".C."(n-k));"); Print(" DPEadd(&s,&tmp);"); Print(" if ((c=DPEcmp(&s,".A."+n))==0) {MFerror(\"FAIL\\n\");}"); Print(" else if (c==1) break;"); Print(" k=n".s0."-k;"); Print(" }"); Print(" if ((n-k)>0)"); Print(" res=MFnewExpr(3, MFident(\"Set\"), draw".B."(k), draw".C."(n-k) );"); Print(" else"); Print(" res=MFnewExpr(2, MFident(\"Set\"), draw".B."(k) );"); end_if; Print(" MFsig( res );"); Print(" return( res );"); end_proc: cs::makedraw2[MODULE,DPE,ProdPrdSet] := subsop(cs::makedraw2[MODULE,DPE,ProdPrdSeq], 6=hold(cs::makedraw2[MODULE,DPE,ProdPrdSet])): #----------------------------------------------------- Products => linear recurrences -----------------------------------------------------# cs::makecountprod[MODULE,DPE]:=proc(sys,algo) begin cs::makecountprod[MODULE,DPE,algo](sys) end_proc: cs::makecountprod[MODULE,DPE,HOLONOMIC]:=proc(sys) local s,pr,npr,p,r,u,n,NT; begin print(Unquoted,"BE CAREFUL !"); print(Unquoted,"Verify the recurrences (in the C file) don't have negative coefficients."); ([pr,npr,s]):=split(sys,cs::isProd[HOLONOMIC],sys); cs::makecount[MODULE,DPE](npr); for p in pr do Print("/* ".expr2text(p)." */"); # count function for A is the same as NT ? # NT:=cs::lequalcount[op(p,1)]; if type(NT)<>"_index" then Print("DPE* count".expr2text(op(p,1))."(int n) {"); Print(" return(count".NT."(n));"); Print("}\n"); else # otherwise # r:=cs::speciftorec(sys,op(p,1),u(n)); userinfo(1,"recurrence for ",p," is: ",expand(r)); cs::makecountprod1[C,DPE,HOLONOMIC](op(p),r,u(n)); Print("}\n"); cs::makeinitrec[MODULE,DPE](op(p,1)); end_if; end_for; end_proc: cs::makeinitrec[MODULE,DPE]:=proc(A) begin Print("void initrec".A."() {"); Print(" if (Nmax<".cs::LIM.") reinit(".cs::LIM.");"); for i from cs::val[A] to cs::LIM do Print(" DPEset(".A."+".i.",".expr2text(cs::tCount[cs::user_spec_NF][A](i)).".0);"); end_for; Print("}\n"); end_proc: #----------------------------------------------------- Products => lazy algorithm of Joris VdH -----------------------------------------------------# cs::makecountprod[MODULE,DPE,LAZY]:=proc(sys) local p; begin cs::error(MODULE,DPE,"lazy products in MODULE with DPE: not yet implemented. Put cs::LAZY to FALSE"); end_proc: #----------------------------------------------------- Count and Draw Functions for C Modules with GMP target=MODULE, library=GMP -----------------------------------------------------# cs::init[MODULE,GMP]:=proc(s,sys,fname,main) local A,eq; begin cs::filehead(s,sys,fname,"GMP"); Print("/* This is a C Dynamic Module for MuPAD */\n"); cs::printgmp(): end_proc: cs::reset[MODULE,GMP]:=cs::filetail: cs::makecount[MODULE,GMP] := proc(sys) local p; begin cs::error(MODULE,GMP,"count MODULE with GMP: not yet implemented."); end_proc: cs::makedraw[MODULE,GMP] := proc(sys) local p; begin cs::error(MODULE,GMP,"draw MODULE with GMP: not yet implemented."); end_proc: #----------------------------------------------------- Products => linear recurrences -----------------------------------------------------# cs::makecountprod[MODULE,GMP]:=proc(sys,algo) begin cs::makecountprod[MODULE,GMP,algo](sys) end_proc: cs::makecountprod[MODULE,GMP,HOLONOMIC]:=proc(sys) begin cs::error(MODULE,GMP,"recurrences MODULE with GMP: not yet implemented."); end_proc: #----------------------------------------------------- Products => lazy algorithm of Joris VdH -----------------------------------------------------# cs::makecountprod[MODULE,GMP,LAZY]:=proc(sys) local p; begin cs::error(MODULE,GMP,"lazy products MODULE with GMP: not yet implemented. Put cs::LAZY to FALSE"); end_proc: