%% Created by Maple 9.00, Mac OS X
%% Source Worksheet: gf-sequences.latte.sc.mw - /Users/zaslav/Tarticles/SLS/Maple/
%% Generated: Wed Mar 10 23:02:46 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}
S:  semimagic squares (all positive values).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.
RtoSfactor:  the rational function that multiplies Rgf to Sgf and rgf to sgf.This is for .\begin{mapleinput}
\mapleinline{active}{1d}{d:=5; p:=60;
RtoSfactor:=x\symbol{94}2/(1-x)\symbol{94}2;}{}
\end{mapleinput}The number of terms desired of each sequence is "enddegree".  Comment:  The slow part of the program is the series expansion (the following four commands).\begin{mapleinput}
\mapleinline{active}{1d}{enddegree:=100;}{}
\end{mapleinput}We start by computing 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 series.  These are  data.\begin{mapleinput}
\mapleinline{active}{1d}{simplexname[1]:="OABC": ehrgf[1]:= 1/((1-x)\symbol{94}3*(1-x\symbol{94}2)) : dimen[1]:=3:
simplexname[2]:="OEE2": ehrgf[2]:= 1/((1-x)*(1-x\symbol{94}2)*(1-x\symbol{94}3)) : dimen[2]:=2:
simplexname[3]:="OAE2": ehrgf[3]:= 1/((1-x)*(1-x\symbol{94}2)\symbol{94}2) : dimen[3]:=2:
simplexname[4]:="ADE2": ehrgf[4]:= 1/((1-x\symbol{94}2)\symbol{94}3) : dimen[4]:=2:
simplexname[5]:="DE1E2": ehrgf[5]:= 1/((1-x\symbol{94}2)\symbol{94}2*(1-x\symbol{94}3)) : dimen[5]:=2:
simplexname[6]:="OCE": ehrgf[6]:= 1/((1-x)\symbol{94}2*(1-x\symbol{94}3)) : dimen[6]:=2:
simplexname[7]:="BDE1": ehrgf[7]:= 1/((1-x)*(1-x\symbol{94}2)*(1-x\symbol{94}3)) : dimen[7]:=2:
simplexname[8]:="ABD": ehrgf[8]:= 1/((1-x)*(1-x\symbol{94}2)\symbol{94}2) : dimen[8]:=2:
simplexname[9]:="FG1": ehrgf[9]:= 1/((1-x\symbol{94}3)*(1-x\symbol{94}5)) : dimen[9]:=1:
simplexname[10]:="EF": ehrgf[10]:= 1/((1-x\symbol{94}3)\symbol{94}2) : dimen[10]:=1:
simplexname[11]:="OG": ehrgf[11]:= 1/((1-x)*(1-x\symbol{94}4)) : dimen[11]:=1:
simplexname[12]:="FG": ehrgf[12]:= 1/((1-x\symbol{94}3)*(1-x\symbol{94}4)) : dimen[12]:=1:
simplexname[13]:="AF": ehrgf[13]:= 1/((1-x\symbol{94}2)*(1-x\symbol{94}3)) : dimen[13]:=1:
simplexname[14]:="DG": ehrgf[14]:= 1/((1-x\symbol{94}2)*(1-x\symbol{94}4)) : dimen[14]:=1:
simplexname[15]:="DG2": ehrgf[15]:= 1/((1-x\symbol{94}2)*(1-x\symbol{94}5)) : dimen[15]:=1:
simplexname[16]:="DE": ehrgf[16]:= 1/((1-x\symbol{94}2)*(1-x\symbol{94}3)) : dimen[16]:=1:
simplexname[17]:="H": ehrgf[17] := 1/(1-x\symbol{94}5) : dimen[17]:=0:}{}
\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}Get the g.f. of reduced, normalized squares.\begin{mapleinput}
\mapleinline{active}{1d}{for n from 1 to 17 do 
  rgfterm[n]:=openehrgf[n]: 
od:
rgfsum:=sum(mu[nn]*rgfterm[nn],nn=1..17):
rgf:=simplify(rgfsum):}{}
\end{mapleinput}Get the g.f. of reduced squares.\begin{mapleinput}
\mapleinline{active}{1d}{Rgfsum:=72*rgfsum:
Rgf:=simplify(Rgfsum):}{}
\end{mapleinput}Hence S, the g.f. of the number of semimagic squares, equals\begin{mapleinput}
\mapleinline{active}{1d}{Sgf:=simplify(RtoSfactor*Rgf):}{}
\end{mapleinput}The g.f. of the total number of symmetry types, l_ml ("lgf"):\begin{mapleinput}
\mapleinline{active}{1d}{sgf:=simplify(RtoSfactor*rgf): }{}
\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}{Sgfdegree:=ldegree(numer(Sgf),x);
Rgfdegree:=ldegree(numer(Rgf),x);
sgfdegree:=ldegree(numer(sgf),x);
rgfdegree:=ldegree(numer(rgf),x);}{}
\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}{Sgfstandnum:=simplify(numer(Sgf)*standenom/denom(Sgf)):
Sgf:=Sgfstandnum/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}{sgfstandnum:=simplify(numer(sgf)*standenom/denom(sgf)):
sgf:=sgfstandnum/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 as a series to find the first few values of the number of squares and symmetry types.\begin{mapleinput}
\mapleinline{active}{1d}{Sseries:=series(Sgf,x=0,enddegree+1):
print("Series computed.");}{}
\end{mapleinput}\begin{mapleinput}
\mapleinline{active}{1d}{sseries:=series(sgf,x=0,enddegree+1):
print("Series computed.");}{}
\end{mapleinput}Expand as a series to find the first few values of the number of reduced squares and reduced symmetry types.\begin{mapleinput}
\mapleinline{active}{1d}{Rseries:=series(Rgf,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}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 Sgfdegree to enddegree do
 co:=coeff(Sseries,x,n):
# printf("%d  %d \n",n,co);
od:
print("Coefficients 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("Coefficients complete.",n,co);}{}
\end{mapleinput}\begin{mapleinput}
\mapleinline{active}{1d}{for n from sgfdegree to enddegree do
 co:=coeff(sseries,x,n):
# printf("%d  %d \n",n,co);
od:
print("Coefficients 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("Coefficients complete.",n,co);}{}
\end{mapleinput}\end{document}
%% End of Maple 9.00 Output