#------------------------------------------------------------------------------
| |
| 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: