%% Created by Maple 9.00, Mac OS X
%% Source Worksheet: gf-sequences.latte.mla.mw - /Users/zaslav/Tarticles/SLS/Maple/
%% Generated: Wed Mar 10 22:24:44 EST 2010
\documentclass{article}
\usepackage{maple2e}
\def\emptyline{\vspace{12pt}}
\DefineParaStyle{maplegroup}
\DefineParaStyle{Author}
\DefineParaStyle{subtitle}
\DefineParaStyle{Bullet Item}
\DefineParaStyle{Dash Item}
\DefineParaStyle{Diagnostic}
\DefineParaStyle{Error}
\DefineParaStyle{Heading 1}
\DefineParaStyle{Heading 2}
\DefineParaStyle{Heading 3}
\DefineParaStyle{Heading 4}
\DefineParaStyle{Normal}
\DefineParaStyle{Text Output}
\DefineParaStyle{Title}
\DefineParaStyle{Warning}
\DefineCharStyle{2D Math}
\DefineCharStyle{Help Heading}
\DefineCharStyle{Maple Input}
\DefineCharStyle{FixedWidth}
\DefineCharStyle{endFixedWidth}
\begin{document}
\pagestyle{empty}
L, S:  magilatin, semimagic squares (all positive values).ml:  magilatin, except in g.f.'s.l, s:   normalized squares (symmetry types).R:   reduced squares (least element is 0).r:   reduced normalized squares (reduced symmetry types).n:  semimagic r.gf:  generating function in some form.gfsum:  generating function as a sum of simple terms.c:   Cubic (fixed strict upper bound; weak upper bound for reduced).a:   Affine (fixed magic sum).p:  Period of the quasipolynomial (known from geometry).  (Period of the truncated quasipolynomial, in the affine count.)
d:  Dimension of the geometry = degree of the quasipolynomials.
RtoLfactor:  the rational function that multiplies Rgf to Lgf and rgf to lgf.enddegree:  The number of terms desired in the sequences, from degree 1 (but initial zeros will be omitted). \begin{mapleinput}
\mapleinline{active}{1d}{enddegree:=100;}{}
\end{mapleinput}This is for : set up main constants.\begin{mapleinput}
\mapleinline{active}{1d}{d:=4; p:=120;
RtoLfactor:=x\symbol{94}3/(1-x\symbol{94}3);}{}
\end{mapleinput}We start by recomputing R_s=rsgf from the semimagic count.  From the Latte results we get the closed Ehrhart g.f. of each flat, which depends on whether we're doing cubic or affine.  Set up the simplex data for the faces and intersection polytopes in the semimagic part of the magilatin series.  This is for \begin{mapleinput}
\mapleinline{active}{1d}{simplexname[1]:="OABC": ehrgf[1]:= 1/((1-x)*(1-x\symbol{94}2)\symbol{94}3) : dimen[1]:=3:
simplexname[2]:="OEE2": ehrgf[2]:= 1/((1-x)*(1-x\symbol{94}4)\symbol{94}2) : dimen[2]:=2:
simplexname[3]:="OAE2": ehrgf[3]:= 1/((1-x)*(1-x\symbol{94}2)*(1-x\symbol{94}4)) : dimen[3]:=2:
simplexname[4]:="ADE2": ehrgf[4]:= 1/((1-x\symbol{94}2)*(1-x\symbol{94}3)*(1-x\symbol{94}4)) : dimen[4]:=2:
simplexname[5]:="DE1E2": ehrgf[5]:= 1/((1-x\symbol{94}3)*(1-x\symbol{94}4)\symbol{94}2) : dimen[5]:=2:
simplexname[6]:="OCE": ehrgf[6]:= 1/((1-x)*(1-x\symbol{94}2)*(1-x\symbol{94}4)) : dimen[6]:=2:
simplexname[7]:="BDE1": ehrgf[7]:= 1/((1-x\symbol{94}2)*(1-x\symbol{94}3)*(1-x\symbol{94}4)) : dimen[7]:=2:
simplexname[8]:="ABD": ehrgf[8]:= 1/((1-x\symbol{94}2)\symbol{94}2*(1-x\symbol{94}3)) : dimen[8]:=2:
simplexname[9]:="FG1": ehrgf[9]:= 1/((1-x\symbol{94}5)*(1-x\symbol{94}8)) : dimen[9]:=1:
simplexname[10]:="EF": ehrgf[10]:= 1/((1-x\symbol{94}4)*(1-x\symbol{94}5)) : dimen[10]:=1:
simplexname[11]:="OG": ehrgf[11]:= 1/((1-x)*(1-x\symbol{94}6)) : dimen[11]:=1:
simplexname[12]:="FG": ehrgf[12]:= 1/((1-x\symbol{94}5)*(1-x\symbol{94}6)) : dimen[12]:=1:
simplexname[13]:="AF": ehrgf[13]:= 1/((1-x\symbol{94}2)*(1-x\symbol{94}5)) : dimen[13]:=1:
simplexname[14]:="DG": ehrgf[14]:= 1/((1-x\symbol{94}3)*(1-x\symbol{94}6)) : dimen[14]:=1:
simplexname[15]:="DG2": ehrgf[15]:= 1/((1-x\symbol{94}3)*(1-x\symbol{94}8)) : dimen[15]:=1:
simplexname[16]:="DE": ehrgf[16]:= 1/((1-x\symbol{94}3)*(1-x\symbol{94}4)) : dimen[16]:=1:
simplexname[17]:="H": ehrgf[17] := 1/(1-x\symbol{94}7) : dimen[17]:=0:
# for n from 1 to 17 do print(simplexname[n], dimen[n], ehrgf[n]); od;}{}
\end{mapleinput}The closed E.g.f. is converted to the open E.g.f.  The first step is to compute the Mobius function of the intersection poset.\begin{mapleinput}
\mapleinline{active}{1d}{for n from 1 to 17 do 
  mu[n]:=(-1)\symbol{94}(dimen[1]-dimen[n]):
od:
mu[14]:=2*mu[14]:
for n from 1 to 17 do  
  openehrgf[n]:=simplify(-(-1)\symbol{94}dimen[n]*subs(x=1/x,ehrgf[n])):
od:}{}
\end{mapleinput}Set up basic g.f.'s.\begin{mapleinput}
\mapleinline{active}{1d}{for n from 1 to 17 do 
  rsgfterm[n]:=openehrgf[n]: 
od:
rsgfsum:=sum(mu[nn]*rsgfterm[nn],nn=1..17):
rsgf:=simplify(rsgfsum):
sgf:=simplify(RtoLfactor*rsgf):}{}
\end{mapleinput}The additional faces and intersection polytopes involved in the magilatin computation.  They depend on whether we're cubic or affine.  These are for . \begin{mapleinput}
\mapleinline{active}{1d}{mlsimplexname[1]:="OAB": mlehrgf[1]:= 1 / ((1-x)*(1-x\symbol{94}2)\symbol{94}2) : mldimen[1]:=2:
mlsimplexname[2]:="OE": mlehrgf[2]:= 1 / ((1-x)*(1-x\symbol{94}4)) : mldimen[2]:=1:
mlsimplexname[3]:="OAC": mlehrgf[3]:= 1 / ((1-x)*(1-x\symbol{94}2)\symbol{94}2) : mldimen[3]:=2:
mlsimplexname[4]:="AD": mlehrgf[4]:= 1 / ((1-x\symbol{94}3)*(1-x\symbol{94}2)) : mldimen[4]:=1:
mlsimplexname[5]:="DE1": mlehrgf[5]:= 1 / ((1-x\symbol{94}3)*(1-x\symbol{94}4)) : mldimen[5]:=1:
mlsimplexname[6]:="OBC": mlehrgf[6]:= 1 / ((1-x)*(1-x\symbol{94}2)\symbol{94}2) : mldimen[6]:=2:
mlsimplexname[7]:="OE2": mlehrgf[7]:= 1 / ((1-x)*(1-x\symbol{94}4)) : mldimen[7]:=1:
mlsimplexname[8]:="BD": mlehrgf[8]:= 1 / ((1-x\symbol{94}2)*(1-x\symbol{94}3)) : mldimen[8]:=1:
mlsimplexname[9]:="DE2": mlehrgf[9]:= 1 / ((1-x\symbol{94}3)*(1-x\symbol{94}4)) : mldimen[9]:=1:
mlsimplexname[10]:="F": mlehrgf[10]:= 1/(1-x\symbol{94}5) : mldimen[10]:=0:
mlsimplexname[11]:="OB": mlehrgf[11]:= 1/((1-x)*(1-x\symbol{94}2)) : mldimen[11]:=1:
# for n from 1 to 11 do print(mlsimplexname[n], mldimen[n], mlehrgf[n]); od;}{}
\end{mapleinput}Now a general computation.  First, open Ehrhart g.f.'s.\begin{mapleinput}
\mapleinline{active}{1d}{for n from 1 to 11 do  
  openmlehrgf[n]:=simplify(-(-1)\symbol{94}mldimen[n]*subs(x=1/x,mlehrgf[n])):
od:}{}
\end{mapleinput}(-1)\symbol{94}3 n_OAB(1/x) equals mlehrgf[1]+mlehrgf[2], and hence n_OAB(x) is, \begin{mapleinput}
\mapleinline{active}{1d}{mlnnew[1] := openmlehrgf[1]-openmlehrgf[2]:}{}
\end{mapleinput}(-1)\symbol{94}3 n_OAC(1/x) equals mlehrgf[3]+mlehrgf[4]+mlehrgf[5].  Hence n_OAC(x) equals\begin{mapleinput}
\mapleinline{active}{1d}{mlnnew[2] := openmlehrgf[3]-openmlehrgf[4]-openmlehrgf[5]:}{}
\end{mapleinput}(-1)\symbol{94}3 n_OBC(1/x) equals mlehrgf[6]+mlehrgf[7]+mlehrgf[8]+mlehrgf[9]+mlehrgf[10].  So n_OBC(x) equals\begin{mapleinput}
\mapleinline{active}{1d}{mlnnew[3] := openmlehrgf[6]-openmlehrgf[7]-openmlehrgf[8]-openmlehrgf[9]+openmlehrgf[10]:}{}
\end{mapleinput}Finally, OB gives mlehrgf[11], so that n_OB(x) is\begin{mapleinput}
\mapleinline{active}{1d}{mlnnew[4] := openmlehrgf[11]:}{}
\end{mapleinput}To compute r, we need rs=n from semimagic, which equals rgf:\begin{mapleinput}
\mapleinline{active}{1d}{Rgfsum:=72*rsgfsum+36*(mlnnew[1]+mlnnew[2]+mlnnew[3])+12*mlnnew[4]:
Rgf:=simplify(Rgfsum);}{}
\end{mapleinput}Hence L, the g.f. of the number of magilatin squares, equals\begin{mapleinput}
\mapleinline{active}{1d}{Lgf:=simplify(RtoLfactor*Rgf);}{}
\end{mapleinput}\begin{mapleinput}
\mapleinline{active}{1d}{rgfsum:=rsgfsum+mlnnew[1]+mlnnew[2]+mlnnew[3]+mlnnew[4]:
rgf:=simplify(rgfsum);}{}
\end{mapleinput}The g.f. of the total number of symmetry types, l_ml ("lgf"):\begin{mapleinput}
\mapleinline{active}{1d}{lgf:=simplify(RtoLfactor*rgf); }{}
\end{mapleinput}Expressing the rational function with standard denominator gives an orders-of-magnitude speedup in the series expansion.  Standard denominator (1-x\symbol{94}p)\symbol{94}{d+1}.  \begin{mapleinput}
\mapleinline{active}{1d}{pdenom:=(1-x\symbol{94}p):
standenom:=pdenom\symbol{94}(d+1);}{}
\end{mapleinput}G.f. as rational function with standard denominator.\begin{mapleinput}
\mapleinline{active}{1d}{Lgfstandnum:=simplify(numer(Lgf)*simplify(standenom/denom(Lgf))):
Lgf:=Lgfstandnum/standenom;}{}
\end{mapleinput}G.f. as rational function with standard denominator.\begin{mapleinput}
\mapleinline{active}{1d}{Rgfstandnum:=simplify(numer(Rgf)*standenom/denom(Rgf)):
Rgf:=Rgfstandnum/standenom;}{}
\end{mapleinput}G.f. as rational function with standard denominator.\begin{mapleinput}
\mapleinline{active}{1d}{lgfstandnum:=simplify(numer(lgf)*simplify(standenom/denom(lgf))): 
lgf:=lgfstandnum/standenom;}{}
\end{mapleinput}G.f. as rational function with standard denominator.\begin{mapleinput}
\mapleinline{active}{1d}{rgfstandnum:=simplify(numer(rgf)*standenom/denom(rgf)):
rgf:=rgfstandnum/standenom;}{}
\end{mapleinput}Expand the series to find the first few values of the number of squares.\begin{mapleinput}
\mapleinline{active}{1d}{Lseries:=series(Lgf,x=0,enddegree+1):
print("Series computed.");}{}
\end{mapleinput}Expand the series to find the first few values of the number of reduced squares.\begin{mapleinput}
\mapleinline{active}{1d}{Rseries:=series(Rgf,x=0,enddegree+1):
print("Series computed.");}{}
\end{mapleinput}Expand the series to find the first few values of the number of symmetry types.\begin{mapleinput}
\mapleinline{active}{1d}{lseries:=series(lgf,x=0,enddegree+1):
print("Series computed.");}{}
\end{mapleinput}\begin{mapleinput}
\mapleinline{active}{1d}{rseries:=series(rgf,x=0,enddegree+1):
print("Series computed.");}{}
\end{mapleinput}Generate the labelled sequence of magilatin square numbers of all four kinds.  The first step is to compute the degree of the first non-zero term.\begin{mapleinput}
\mapleinline{active}{1d}{Lgfdegree:=ldegree(numer(Lgf),x);
Rgfdegree:=ldegree(numer(Rgf),x);
lgfdegree:=ldegree(numer(lgf),x);
rgfdegree:=ldegree(numer(rgf),x);}{}
\end{mapleinput}List the coefficients of each series, i.e., the terms of the counting sequences.  With large "enddegree" the output is huge so it's more convenient to run each sequence's output separately and copy it from the worksheet.\begin{mapleinput}
\mapleinline{active}{1d}{for n from Lgfdegree to enddegree do
 co:=coeff(Lseries,x,n):
# printf("%d  %d \n",n,co);
od:
print("Sequence complete.",n,co);}{}
\end{mapleinput}\begin{mapleinput}
\mapleinline{active}{1d}{for n from Rgfdegree to enddegree do
 co:=coeff(Rseries,x,n):
# printf("%d  %d \n",n,co);
od:
print("Sequence complete.",n,co);}{}
\end{mapleinput}\begin{mapleinput}
\mapleinline{active}{1d}{for n from lgfdegree to enddegree do
 co:=coeff(lseries,x,n):
# printf("%d  %d \n",n,co);
od:
print("Sequence complete.",n,co);}{}
\end{mapleinput}\begin{mapleinput}
\mapleinline{active}{1d}{for n from rgfdegree to enddegree do
 co:=coeff(rseries,x,n):
# printf("%d  %d \n",n,co);
od:
print("Sequence complete.",n,co);}{}
\end{mapleinput}\end{document}
%% End of Maple 9.00 Output