print(`This is AppsWZmulti, one of the two Maple packages`): print(`accompanying the article `): print(`Some Enumerative and Probabilstic Applications `): print(` of Wilf-Zeilberger Theory `): print( ` by Moa Apagodu and Doron Zeilberger.`): lprint(``): print(`The most current version is available on WWW at:`): print(` http://www.math.rutgers.edu/~zeilberg .`): print(`Please report all bugs to: zeilberg at math dot rutgers dot edu .`): lprint(``): print(`For general help, and a list of the available functions,`): print(` type "ezra();". For specific help type "ezra(procedure_name)" `): lprint(``): ###########Just MAZ############# ezraMAZ:=proc() if nargs=0 then print(`MultiAlmkvistZeilberger`): print(`A Maple package for finding recurrences for multi-integrals of `): print(` Hypergeometric/Hyperexponential Integrands.`): print(`For help with a specific procedure, type "ezra(procedure_name);"`): print(`Contains procedures: MAZ `): elif nargs=1 and args[1]=MAZ then print(`MAZ(POL,H,rat,x,n,N,pars): Given a polynomial, POL, an hyper-exponential function, H, and`): print(`a rational function, rat, of the continuous variables in the list x, and a discrete variable`): print(`n, and the shift-operator in n, N, and a set of auxiliary parameters, par`): print(`ouputs a recurrence operator, let's call it ope(n,n), and a multi-certificate, let's call it R`): print(`such that the function F(n;x):=POL*H*rat^n satisfies`): print(` ope(n,N)F=sum(D_{x[i]} (R[i]H*rat^n),i=1..nops(R)`): print(`In particular, try:`): print(`MAZ(1,((1-1/x)*(1-y))^b*((1-1/x/y)*(1-1/y))^c/x/y,(1-x)*(1-x*y),[x,y],a,A,{b,c});`): print(`MAZ(1,exp(-x^2/2-y^2/2),(x-y)^2,[x,y],n,N,{});`): else print(`There is no help for`, args): fi: end: ez:=proc():print(`MAZ1(POL,H,rat,x,n,N,L,pars), MAZ(POL,H,rat,x,n,N,pars)`): print(`FindDeg(POL,H,rat,x,n,L)`): print(`CheckMAZ1(POL,H,rat,x,n,N,L,pars)`): end: _EnvAllSolutions:=true: ####From MultiZeilberger #IV1(d,n): all the vectors of non-negative integres of length d whose #sum is n IV1:=proc(d,n) local gu,i,gu1,i1: if d=0 then if n=0 then RETURN({[]}): else RETURN({}): fi: fi: gu:={}: for i from 0 to n do gu1:=IV1(d-1,n-i): gu:=gu union {seq([op(gu1[i1]),i],i1=1..nops(gu1))}: od: gu: end: yafe:=proc(ope,N) local i: add(factor(coeff(ope,N,i))*N^i,i=ldegree(ope,N)..degree(ope,N)):end: #IV(d,n): all the vectors of non-negative integres of length d whose #sum is <=n IV:=proc(d,n) local i: {seq(op(IV1(d,i)),i=0..n)}: end: #GenPol(kList,a,deg): The generic polynomial in kList of #degree deg, using the indexed variable a, followed by the set #of coeffs. GenPol:=proc(kList,a,deg) local gu,i,i1: gu:=IV(nops(kList),deg): convert([seq(a[i]*convert([seq(kList[i1]^gu[i][i1],i1=1..nops(kList))],`*`), i=1..nops(gu))],`+`),{seq(a[i],i=1..nops(gu))}: end: #Extract1(POL,kList): extracts all the coeffs. of a POL in kList Extract1:=proc(POL,kList) local k1,kList1,i: k1:=kList[nops(kList)]: kList1:=[op(1..nops(kList)-1,kList)]: if nops(kList)=1 then RETURN({seq(coeff(POL,k1,i),i=0..degree(POL,k1))}): else RETURN({seq( op(Extract1(coeff(POL,k1,i),kList1)), i=0..degree(POL,k1))}): fi: end: ###########End from MultiZeilberger FindDeg:=proc(POL,H,rat,x,xSet,n,L) local s,t,h,e,i,Hbar,q,r,gu: s:=numer(rat): t:=denom(rat): h:=convert([seq(e[i]*subs(n=n+i,POL)*s^i*t^(L-i),i=0..L)],`+`): Hbar:=H*s^n/t^(n+L): gu:=normal(simplify(diff(Hbar,x)/Hbar)): q:=numer(gu): r:=denom(gu): max(degree(h,xSet)-degree( diff(r,x)+q,xSet)+degree(r,xSet) ,degree(h,xSet)-degree(r,xSet)+1+degree(r,xSet)): end: #MAZ1(POL,H,rat,x,n,N,L,pars): Given a polynomial POL of the discrete variable n #and the continuous variables in x #and a pure-hyperexponential function H of the variables in x, and a rational function #Rat of x, and a shift operator N, and a non-neg. integer L #returns FAIL or, if there is one, an operator ope(N,n) of order L and #a list of rational functions R, such that F:=POL*H*rat^n satisfies #ope(N,n)F=sum(diff(F*R[i],x[i]),i=1..nops(R)) . #For example, try MAZ1(1,exp(-x),x,[x],n,N,1,{}) MAZ1:=proc(POL,H,rat,x,n,N,L,pars) local deg,deg1,m,gu,i,i1: deg:=[]: for i from 1 to nops(x) do deg:=[op(deg),FindDeg(POL,H,rat,x[i],{seq(x[i1],i1=1..nops(x))},n,L)]: od: m:=min(op(deg)): for i from 0 to m do deg1:=[seq(deg[i1]-m+i,i1=1..nops(deg))]: gu:=MAZ1deg(POL,H,rat,x,n,N,L,pars,deg1): if gu<>FAIL then RETURN(gu): fi: od: FAIL: end: MAZ:=proc(POL,H,rat,x,n,N,pars) local gu,L: for L from 0 do gu:=MAZ1(POL,H,rat,x,n,N,L,pars): if gu<>FAIL then RETURN(gu): fi: od: end: CheckMAZ:=proc(POL,H,rat,x,n,N,pars) local F,gu,luL,luR,R,ope,i: F:=H*rat^n: gu:=MAZ(POL,H,rat,x,n,N,pars): if gu=FAIL then RETURN(FAIL): fi: ope:=gu[1]: R:=gu[2]: luL:=normal(add(subs(n=n+i,POL)*coeff(ope,N,i)*normal(simplify(subs(n=n+i,F)/F)),i=0..degree(ope,N))): luR:=normal(add( normal(diff(R[i]*F,x[i])/F ), i=1..nops(R))): evalb(normal(luL/luR)=1): end: #MAZ1deg(POL,H,rat,x,n,N,L,pars,deg): Given a polynomial POL of the discrete variable n #and the continuous variables in x #and a pure-hyperexponential function H of the variables in x, and a rational function #Rat of x, and a shift operator N, and a non-neg. integer L #returns FAIL or, if there is one, an operator ope(N,n) of order L and #a list of rational functions R, such that F:=POL*H*rat^n satisfies #ope(N,n)F=sum(diff(F*R[i],x[i]),i=1..nops(R)) . #For example, try MAZ1(1,exp(-x),x,[x],n,N,1,{}) MAZ1deg:=proc(POL,H,rat,x,n,N,L,pars,deg) local s,t,h,e,i,Hbar,gu,X,a,var,q,r,eq,ope,var1,ku,i1,eqN,var1N,opeN,bilti,meka: s:=numer(rat): t:=denom(rat): ope:=add(e[i]*N^i,i=0..L): h:=convert([seq(e[i]*subs(n=n+i,POL)*s^i*t^(L-i),i=0..L)],`+`): Hbar:=H*s^n/t^(n+L): gu:=[]: for i from 1 to nops(x) do gu:=[op(gu),normal(simplify(diff(Hbar,x[i])/Hbar))]: od: q:=[]: r:=[]: for i from 1 to nops(gu) do q:=[op(q),numer(gu[i])]: r:=[op(r),denom(gu[i])]: od: var:={}: X:=[]: for i from 1 to nops(x) do ku:=GenPol(x,a[i],deg[i]): X:=[op(X),ku[1]]: var:=var union ku[2]: od: var:=var union {seq(e[i],i=0..L)}: gu:=h: for i from 1 to nops(x) do gu:=expand(normal(gu-(diff(r[i],x[i])+q[i])*X[i]/r[i]-r[i]*diff(X[i]/r[i],x[i]))): od: eq:=Extract1(numer(gu),x): eqN:=subs(n=9/17,eq): eqN:=subs({seq(pars[i]=1/(i^2+3),i=1..nops(pars))},eqN): var1N:=solve(eqN,var): opeN:=subs(var1N,ope): if opeN=0 then RETURN(FAIL): fi: var1:=solve(eq,var): ope:=subs(var1,ope): X:=subs(var1,X): if ope=0 then RETURN(FAIL): fi: bilti:={}: for i from 1 to nops(var1) do if op(1,op(i,var1))=op(2,op(i,var1)) then bilti:=bilti union {op(1,op(i,var1))}: fi: od: for i from 1 to nops(bilti) do if coeff(ope,bilti[i],1)<>0 then ope:=coeff(ope,bilti[i],1): X:=[seq(coeff(X[i1],bilti[i],1),i1=1..nops(X))]: meka:=coeff(ope,N,degree(ope,N)): ope:=expand(normal(ope/meka)): ope:=yafe(ope,N): X:=[seq(normal(X[i1]/meka/t^L),i1=1..nops(X))]: RETURN(ope,X): fi: od: end: ##########End MAZ############ ####Begin Asymptotics with(SolveTools): Yafe:=proc(ope,N) local i,ope1,coe1,L: if ope=0 then RETURN(1,0): fi: ope1:=expand(ope): L:=degree(ope1,N): coe1:=coeff(ope1,N,L): ope1:=normal(ope1/coe1): ope1:=normal(ope1): ope1:= convert( [seq(factor(coeff(ope1,N,i))*N^i,i=ldegree(ope1,N)..degree(ope1,N))],`+`): factor(coe1),ope1: end: #Aope1(ope,n,N): the asymptotics of the difference operator with poly. coeffs. ope(n,N) Aope1:=proc(ope,n,N) local gu: gu:=expand(numer(normal(ope))): gu:=coeff(gu,n,degree(gu,n)): gu:=expand(gu/coeff(gu,N,degree(gu,N))): factor(gu): end: C3:=proc(a,b,c) option remember: if [a,b,c]=[0,0,0] then RETURN(1): fi: if a<0 or b<0 or c<0 then RETURN(0): fi: C3(a-1,b,c)+C3(a,b-1,c)+C3(a,b,c-1)-4*C3(a-1,b-1,c-1): end: #Aluf1(pit): is it positive dominant? Aluf1:=proc(pit) local aluf,shi,i,makom: shi:=evalf(abs(pit[1])): aluf:={evalf(evalc(pit[1]))}: makom:={1}: for i from 2 to nops(pit) do if evalf(abs(pit[i]))>evalf(shi) then shi:=abs(evalf(evalc(pit[i]))): aluf:={pit[i]}: makom:={i}: elif evalf(abs(pit[i]))=shi then aluf:=aluf union {evalf(evalc(pit[i]))}: makom:=makom union {i}: fi: od: aluf,shi,makom: end: #OneStepAS1(ope1,n,N,alpha,f,S1): Given the partial asymptotic #expansion of solutions of ope1(n,N) with exponent alpha #extends it by one term OneStepAS1:=proc(ope1,n,N,alpha,f,S1) local x,f1,L,F,A,mu,i,A1: f1:=subs(n=1/x,f): L:=degree(f1,x): F:=f1+A*x^(L+S1): mu:=add( subs(n=1/x,coeff(ope1,N,i))*(1+i*x)^alpha*subs(x=x/(1+i*x),F) ,i=ldegree(ope1,N)..degree(ope1,N)): mu:=normal(mu): mu:=taylor(mu,x=0,L+11): mu:=simplify(mu): for i from L+1 to L+9 while coeff(mu,x,i)=0 do od: if i=L+10 then RETURN(FAIL): fi: mu:=coeff(mu,x,i): A1:=subs(Linear({mu},{A}),A): if A1=0 then RETURN(FAIL): fi: subs({x=1/n,A=A1},F): end: #OneStepA(ope1,n,N,alpha,f): Given the partial asymptotic #expansion of solutions of ope1(n,N) with exponent alpha #extends it by one term OneStepA:=proc(ope1,n,N,alpha,f) local S1: for S1 from 1 while OneStepAS1(ope1,n,N,alpha,f,S1)=FAIL do od: OneStepAS1(ope1,n,N,alpha,f,S1): end: #Asy(ope,n,N,K): the asymptotic expansion of solutions #to ope(n,N)f(n)=0,where ope(n,N) is a recurrence operator #up to the K's term Asy:=proc(ope,n,N,K) local gu,pit,lu,makom,alpha,mu,ope1,ku,i,f,x,ka: gu:=Aope1(ope,n,N): if degree(gu,N)=0 then print(`Not of exponential growth, case not implemented`): RETURN(FAIL): fi: if not type(expand(normal(ope/gu)),`+`) then pit:=[solve(gu,N)]: makom:=Aluf1(pit)[3]: RETURN(pit[makom[1]]^n*expand((nops(makom)-1)!*binomial(nops(makom)-1+n,nops(makom)-1))): fi: pit:=[solve(gu,N)]: makom:=Aluf1(pit)[3]: if nops(makom)<>1 then lu:=evalf(evalc(pit[makom[1]])): if coeff(lu,I,1)>(1/10)^(Digits-2) then print(`Not only that but Dominant roots are complex`): RETURN(FAIL): fi: fi: lu:=evalf(evalc(pit[makom[1]])): if coeff(lu,I,1)>(1/10)^(Digits-2) then print(`Dominant root is complex`): RETURN(FAIL): fi: lu:=pit[makom[1]]: ope1:=Yafe(subs(N=lu*N,ope),N)[2]: mu:=add( subs(n=1/x,coeff(ope1,N,i))*(1+i*x)^alpha,i=ldegree(ope1,N)..degree(ope1,N)): mu:=normal(mu): ku:=factor(coeff(taylor(mu,x=0,nops(makom)+1),x,nops(makom))): ka:=simplify([solve(ku,alpha)]): alpha:=max(op(ka)): if normal(subs(N=1,ope1))=0 then RETURN(lu^n*n^alpha): fi: f:=1: for i from 1 to K do f:=OneStepA(ope1,n,N,alpha,f): od: lu^n*n^alpha*f: end: #SeqFromRec(ope,n,N,Ini,K): Given the first L-1 #terms of the sequence Ini=[f(1), ..., f(L-1)] #satisfied by the recurrence ope(n,N)f(n)=0 #extends it to the first K values SeqFromRec:=proc(ope,n,N,Ini,K) local ope1,gu,L,n1,j1: ope1:=Yafe(ope,N)[2]: L:=degree(ope1,N): if nops(Ini)<>L then ERROR(`Ini should be of length`, L): fi: ope1:=expand(subs(n=n-L,ope1)/N^L): gu:=Ini: for n1 from nops(Ini)+1 to K do gu:=[op(gu), -add(gu[nops(gu)+1-j1]*subs(n=n1,coeff(ope1,N,-j1)), j1=1..L)]: od: gu: end: #AsyC(ope,n,N,K,Ini,L): the asymptotic expansion #complete with the estimated constant, of solutions #to ope(n,N)f(n)=0, with initial conditions #Ini, where ope(n,N) is a recurrence operator #up to the K's term. The constant is estimated using up to L terms #For example, try #AsyC(N^2-N-1,n,N,1,[1,1],100); AsyC:=proc(ope,n,N,K,Ini,L) local lu,gu,C,dig,i: Digits:=100: if nops(Ini)1 then ERROR(`Bad input`): fi: if {seq(type(L[i][2],numeric),i=1..nops(L))}={true} and add(L[i][2],i=1..nops(L))<>1 then ERROR(`Bad input`): fi: n:=gu[1]: add(L[i][2]*mul(x[j]^L[i][1][j],j=1..n),i=1..nops(L)): end: #RandomWalkRecurrence(L,n,N): given a list L of the form [x,p] where #x is a step and p is the probability that it will happen, #finds a linear recurrence equation with constant coefficients #(with N being the shift in n) for the probab. of #coming back to the initial point after n steps. #For example, try: #RandomWalkRecurrence([[1,0],1/4],[[-1,0],1/4],[[0,1],1/4],[[0,-1],1/4] ],n,N); RandomWalkRecurrence:=proc(L,n,N) local gu,x,i,d: gu:={seq(nops(L[i][1]),i=1..nops(L))}: if nops(gu)<>1 then ERROR(`Bad input`): fi: if {seq(type(L[i][2],numeric),i=1..nops(L))}={true} and add(L[i][2],i=1..nops(L))<>1 then ERROR(`Bad input`): fi: d:=gu[1]: MAZ(1,1/mul(x[i],i=1..d), Mpgf(L,x) ,[seq(x[i],i=1..d)],n,N,{})[1]: end: #RandomWalkRecurrenceE(L,n,N): given a list L of the form [x,p] where #x is a step and p is the probability that it will happen, #finds a linear recurrence equation with constant coefficients #(with N being the shift in n) for the probab. of #coming back to the initial point after 2n steps. #For example, try: #RandomWalkRecurrenceE #([[1,0],1/4],[[-1,0],1/4],[[0,1],1/4],[[0,-1],1/4] ],n,N); RandomWalkRecurrenceE:=proc(L,n,N) local gu,x,i,d: gu:={seq(nops(L[i][1]),i=1..nops(L))}: if nops(gu)<>1 then ERROR(`Bad input`): fi: if {seq(type(L[i][2],numeric),i=1..nops(L))}={true} and add(L[i][2],i=1..nops(L))<>1 then ERROR(`Bad input`): fi: d:=gu[1]: MAZ(1,1/mul(x[i],i=1..d), Mpgf(L,x)^2 ,[seq(x[i],i=1..d)],n,N,{})[1]: end: #LatticePaths(Steps,m,M): Given a set of steps (all in the same #dimensional space of course) and letters m, and M, #outputs ordinary #recurrences (in each direction) #for the function F(m):=number of ways of walking #from the origin to the point m =(m[1],m[2], ...) #using the steps in Steps #For example, try: LatticePaths({[0,1],[1,0],[1,1]},m,M); LatticePaths:=proc(Steps,m,M) local gu,i,d,x,j,i1: gu:={seq(nops(Steps[i]),i=1..nops(Steps))}: if nops(gu)<>1 then ERROR(`Bad input`): fi: d:=gu[1]: gu:=1/(1-add(mul(x[j]^Steps[i][j],j=1..d),i=1..nops(Steps))): [seq( MAZ(1,gu/mul(x[i1],i1=1..d)/ mul(x[i1]^m[i1],i1=1..i-1)/mul(x[i1]^m[i1],i1=i+1..d) , 1/x[i] ,[seq(x[i1],i1=1..d)],m[i],M[i],{})[1],i=1..d)]: end: #LatticePathsDiagonal(Steps,n,N): Given a set of steps (all in the same #dimensional space of course) and letters n, and N, #outputs an ordinary #recurrenc #for the function F(n):=number of ways of walking #from the origin to the point (n,n,..., n) #using the steps in Steps #For example, try: LatticePathsDiagonal({[0,1],[1,0],[1,1]},n,N); LatticePathsDiagonal:=proc(Steps,n,N) local gu,i,d,x,j,i1,ope,L: gu:={seq(nops(Steps[i]),i=1..nops(Steps))}: if nops(gu)<>1 then ERROR(`Bad input`): fi: d:=gu[1]: gu:=1/(1-add(mul(x[j]^Steps[i][j],j=1..d),i=1..nops(Steps))): ope:=MAZ(1,gu/mul(x[i1],i1=1..d), 1/mul(x[i1],i1=1..d), [seq(x[i1],i1=1..d)],n,N,{})[1]: L:=degree(ope,N): #HERE ope,[seq(Mekadem(gu,x,[i1$d]),i1=0..L)]: end: #LatticePathsDiagonalStory(Steps,n): #Like LatticePathsDiagonal(Steps,n,N) but the output is #in narrative form #For example, try: LatticePathsDiagonalStory({[0,1],[1,0],[1,1]},n); LatticePathsDiagonalStory:=proc(Steps,n) local N,ope,d,gu,i,ope1: gu:={seq(nops(Steps[i]),i=1..nops(Steps))}: if nops(gu)<>1 then ERROR(`Bad input`): fi: d:=gu[1]: ope1:=LatticePathsDiagonal(Steps,n,N): ope:=ope1[1]: print(`Let F(n) be the number of ways of walking from `, [0$d], `to the `): print(`point `, [n$d], `in the `, d, `-dimensional cubic lattice `): print(`using the following allowed steps:`): print(Steps): print(): print(`F(n) satisfies the following linear recurrence equation`): print(`with polynomial coefficients`): print(add(coeff(ope,N,i)*F(n+i),i=0..degree(ope,N))=0): print(``): print(`subject to the initial conditions`): print(seq(F(i1-1)=ope1[2][i1],i1=1..nops(ope1[2]))): print(`This implies, thanks to Birkhoff-Trijinski, that F(n) is`): print(`asymptotically a constant times `): gu:=AsyC(ope,n,N,3,[op(2..nops(ope1[2]),ope1[2])],100): print(gu): print(`which is roughly equal to`): print(evalf(gu)): print(): print(`For the record, the first 31 terms of the sequence are`): print([1,op(SeqFromRec(ope,n,N,[op(2..nops(ope1[2]),ope1[2])],30))]): end: #LatticePathsStory(Steps,m): #Like LatticePaths(Steps,m,M) but the output is #in narrative form #For example, try: LatticePathsStory({[0,1],[1,0],[1,1]},m); LatticePathsStory:=proc(Steps,m) local M, d,gu,i,i1,j,Opes,mu,ope: gu:={seq(nops(Steps[i]),i=1..nops(Steps))}: if nops(gu)<>1 then ERROR(`Bad input`): fi: d:=gu[1]: Opes:=LatticePaths(Steps,m,M): print(`Let`, F([seq(m[i],i=1..d)]), ` be the number of ways of `): print(` of walking from `): print( [0$d] ): print(`to the point `, m=[seq(m[i],i=1..d)]): print(`in the `, d, `-dimensional cubic lattice `): print(`using the following allowed steps:`): print(Steps): print(): print(F(seq(m[i],i=1..d) ) , `satisfies the following linear recurrences`): print(` equation`): print(`with polynomial coefficients`): for i from 1 to d do print(`in the`, m[i], `direction, it satisfies `): ope:=Opes[i]: mu:=add( coeff(ope,M[i],i1)*F(seq(m[j],j=1..i-1),m[i]+i1,seq(m[j],j=i+1..d)), i1=0..degree(ope,M[i])): print(mu=0): od: end: #Mekadem(F,x,a): the coeff. of x^a in the Maclaurin expansion #of F, where F depends on x[1], ..., x[n] (n=nops(a)): #For example, try: Mekadem(1/(1-x[1]-x[2]),x,[1,1]); Mekadem:=proc(F,x,a) local n,F1,a1: n:=nops(a): if n=0 then RETURN(FAIL): fi: F1:=taylor(F,x[n]=0,a[n]+1): F1:= coeff(F1,x[n],a[n]): if n=1 then RETURN(F1): fi: a1:=[op(1..n-1,a)]: Mekadem(F1,x,a1): end: #LatticePathsDiagonalStory1(Steps,n): #Like LatticePathsDiagonalStory(Steps,n,N) but with #lprint LatticePathsDiagonalStory1:=proc(Steps,n) local N,ope,d,gu,i,ope1: gu:={seq(nops(Steps[i]),i=1..nops(Steps))}: if nops(gu)<>1 then ERROR(`Bad input`): fi: d:=gu[1]: ope1:=LatticePathsDiagonal(Steps,n,N): ope:=ope1[1]: lprint(`Let F(n) be the number of ways of walking from `, [0$d], `to the `): lprint(`point `, [n$d], `in the `, d, `-dimensional cubic lattice `): lprint(`using the following allowed steps:`): lprint(Steps): print(): lprint(`F(n) satisfies the following linear recurrence equation`): lprint(`with polynomial coefficients`): lprint(add(coeff(ope,N,i)*F(n+i),i=0..degree(ope,N))=0): print(``): lprint(`subject to the initial conditions`): lprint(seq(F(i1)=ope1[2][i1],i1=1..nops(ope1[2]))): print(`This implies, thanks to Birkhoff-Trijinski, that F(n) is`): print(`asymptotically a constant times `): lprint(`The first 30 terms of the sequence are`): lprint([1,op(SeqFromRec(ope,n,N,[op(2..nops(ope1[2]),ope1[2])],30))]): gu:=AsyC(ope,n,N,3,[op(2..nops(ope1[2]),ope1[2])],100): lprint(gu): print(`which is roughly equal to`): lprint(evalf(gu)): end: