(* "Copyright (c) 1996 The Regents of the University of California.
     All rights reserved.
  
     Permission to use, copy, modify, and distribute this software 
     and its documentation for any purpose, without fee, and without 
     written agreement is hereby granted, provided that the above 
     copyright notice and the following two paragraphs appear in 
     all copies of this software.
  
     IN NO EVENT SHALL THE UNIVERSITY OF CALIFORNIA BE LIABLE TO 
     ANY PARTY FOR DIRECT, INDIRECT, SPECIAL, INCIDENTAL, OR 
     CONSEQUENTIAL DAMAGES ARISING OUT OF THE USE OF THIS SOFTWARE 
     AND ITS DOCUMENTATION, EVEN IF THE UNIVERSITY OF CALIFORNIA 
     HAS BEEN ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
  
     THE UNIVERSITY OF CALIFORNIA SPECIFICALLY DISCLAIMS ANY 
     WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED 
     WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR 
     PURPOSE.  THE SOFTWARE PROVIDED HEREUNDER IS ON AN "AS IS" 
     BASIS, AND THE UNIVERSITY OF CALIFORNIA HAS NO OBLIGATION 
     TO PROVIDE MAINTENANCE, SUPPORT, UPDATES, ENHANCEMENTS, OR 
     MODIFICATIONS."
 
     Author:              OMER EGECIOGLU (omer@cs.ucsb.edu)
     Version:             1
     Creation Date:       May 1996
     Filename:            DiophantineGF.m                           
                                                                     *)
DiophantineGF[m_, b_, c_] :=

  (* Input is of the form a, b, c where a is the rxs input 
     matrix, b and c are rx1 vectors, and the triple corresponds 
     to the diophantine systems ax = bn+c, x>=0.                     *)

(    
     s = Length[m[[1]]]; 
  (* s is the number of variables = the number of columns of 
     the nonaugmented a.                                             *)

     r = Length[Transpose[m][[1]]]; 
  (* r is the number of equations = the number of rows of 
     the nonaugmented a.                                             *)

     temp = Transpose[Append[Transpose[m], -b]]; 
  (* Add another column to m, originally rxs. This is the s+1 
      column, and consists of -b.                                    *)

     texps = Table[0, {s + 1}]; texps[[-1]] = 1; 
     a = Prepend[temp, texps]; 
  (* Add an initial row marking t's: 0 0 0 ... 0 1.                  *)

     lambda = Array[l, r]; 
  (* l[1], ... , l[r] are the auxilary variables.                    *)

     fudgefactor = Product[lambda[[i]]^(-c[[i]]), {i, 1, r}]; 

     IterationM = Max[Abs[c]];
  (* When power series involving mixed terms (l[i]s and possibly 
     t's) are expanded, the series can be truncated past power 
     IterationM.  IterationM is only used as an upper bound, not 
     the exact cutoff point.                                         *)
 
     gf = 0; 
     sign = +1; 
  
   CheckAbort[recurse[a, sign], gf= Infinity]; 

   (* Ordinarily, this statement would be "recurse[a, sign];" only
      to start the computation at the root of the tree.
      However, when infinite number of solutions are detected via 
      the statement "If[checkinfiniteQ[a, u, v], Abort[ ] ];" in 
      procedure recurse, all computation/recursion is aborted, and
      the control returns to assign gf= Infinity in this statement.  *)

     gf = Factor[gf]
  (* Final output line.                                              *)

)



checkinfiniteQ[a_, u_, v_ ]:=  
  (a[[1]][[u]] == a[[1]][[v]] == 0) && 
     (!zerocolumnQ[a, u]) && (zerocolumnQ[ addcolumn[a, u, v], u ])
  (* checkinfiniteQ: returns True if the u-th and the v-th
     columns of a have opposite signs (and nonzero). The t 
     terms should not be there, this means the entries in the 
     first row must be zero.  The rows that are processed by 
     the called functions are in question are 2 through r+1.        *) 

prunetreeQ[a_] := 
  Module[{i}, prun = False;
    For[i = 2, i <= r+1, i++,
      If[ (Max[ a[[i]] ] < 0 && c[[i-1]] > 0) || 
           (Min[ a[[i]] ] > 0 && c[[i-1]] > 0),  prun = True]
       ]; prun ]
   (* prunetreeQ: returns True if the recursion can be stopped 
     at a.                                                           *)

truncateexp[a_, j_] := 
  Module[{i}, minimum= IterationM;
    For[i = 2, i <= r + 1, i++, 
        If[a[[i]][[j]] != 0 && minimum > Floor[c[[i-1]]/a[[i]][[j]] ],
          minimum =  Floor[ c[[i-1]]/ a[[i]][[j]] ]
          ]
       ];minimum] 
   (* truncateexp: returns the exponent past which the term coming 
      from column j of the leaf matrix a can have no contribution.
      This exponent is the minimum of lambda_i/c_i. Returned in 
      minimum.                                                       *)

zerocolumnQ[a_, j_] := 
  Module[{i}, ans := True; For[i = 2, i <= r + 1, i++, 
  If[a[[i]][[j]] != 0, ans := False]]; ans];
  (* zerocolumnQ: checks to see if the entries in the j-th column 
     of the augmented matrix a are all zeros, excepting the first 
     row. Thus the row indices run from 2 to r+1.                    *)

cleanrowQ[a_, i_] := Max[a[[i]]] <= 0 || Min[a[[i]]] >= 0;
  (* cleanrowQ: returns True iff the i-th row of a is all nonnegative 
     or all nonpositive (i.e. "clean"). All zeros return True.       *)
 
uindex[a_, i_]:= Min[ Position [ a[[i]], Min[a[[i]]] ] ];
vindex[a_, i_]:= Min[ Position [ a[[i]], Max[a[[i]]] ] ];
  (* uindex: scan the i-th row. Return the first index of a minimum 
     element of the row,
     vindex: scan the i-th row. Return the first index of a maximum 
     element of the row.  The elements returned are distinct if 
     cleanrowQ[a,i] = False.                                         *)
 
leafnode[a_] := 
  Module[{i}, i = 2; 
  While[i <= r + 1 && cleanrowQ[a, i], i++]; dirtyrow = i];
  (* leafnode: returns the index of the first dirty row in the 
     augmented matrix a. The relevant rows of a examined are 2 
     through r+1. The function returns r+2 if all the rows are  
     clean, i.e. a is a leaf node of the recursion.                  *)
 
addcolumn[a_, u_, v_] := 
  (at = Transpose[a]; at[[u]] = at[[u]] + at[[v]]; Transpose[at]);
  (* addcolumn: returns a matrix which is obtained from a by adding 
     its v-th column to its u-the column.                            *)

zapcolumn[a_, v_] := Transpose[Drop[Transpose[a], {v, v}]];
  (* zapcolumn: returns the matrix obtained by omitting the v-th 
     column of a.                                                    *)

help:= ( Print[" "];
Print["DiophantineGF[a, b, c]: Calculates the rational generating function"];
Print["        and the formula for the number of solutions to the linear "];
Print["        Diophantine system ax= bn+c, where the right hand side is "];
Print["        parametrized by n = 0, 1, 2, ... "];
Print[" "];
Print["        Inputs:          a :  rxs integral matrix, "];
Print["                   b and c :  rx1 integral vectors."];
Print[" "];
Print["        The unknowns x_1, x_2, ..., x_s are nonnegative integers."];
Print["    "];
Print["        After the computation of the generating function, typing"];
Print["        \"formula\" will generate the formula for the coefficient of "];
Print["        t^n in terms of the binomial coefficients C[x, k] where "];
Print[" "];
Print["                C[x, k]  = Binomial[x, k] if x is a positive integer,"];
Print["                         = zero otherwise."];
Print[" "];
Print["        A power basis expansion is also given if no argument x"];
Print["        is fractional in the formula found."];
Print[" "];
Print["        After the formula has been calculated, \"formulaN[n]\""];
Print["        will evaluate it for a specific argument n."];
       )
example:= ( Print[" "];
Print["EXAMPLE I:"]; Print["---------- "]; 
Print["  "];
Print["   In[2]:= a = {{1, 1, 1, 1, 1}};  b = {2};  c = {-1};"]; 
Print["  "];
Print["   In[3]:= DiophantineGF[a, b, c]"]; Print[" "];
Print["                           2"];
Print["            t (5 + 10 t + t )"];
Print["   Out[3]= -----------------"];
Print["                       5"];
Print["                (1 - t)"];
Print[" "];
Print["   In[4]:= formula"];
Print[" "];
Print["   Binomial Formula : C[1 + n, 4] + 10 C[2 + n, 4] + 5 C[3 + n, 4]"];
Print[" "];
Print["                      n (1 + n) (1 + 2 n) (3 + 2 n)"];
Print["   Power Formula    : -----------------------------"];
Print["                                    6"];
Print[" "];
Print["   In[5]:= formulaN[3]"];
Print["   Out[5]= 126"];
Print[" "];
Print["-----------"]; 
Print["   Type \"example2\" for another example."]; 
          )
example2:= ( Print[" "];
Print["EXAMPLE II:  "]; Print["-----------"]; 
Print["  "];
Print["   In[2]:= a = {{1, 0, 0, 1, 0, 0}, {0, 1, 0, 0, 1, 0}, {0, 0, 1, 0, 0, 1}};"];
Print["   In[3]:= b = {1, 1, 1};  c = {-1, -1, -1};"];
Print["  "];
Print["   In[4]:= DiophantineGF[a, b, c]"];
Print["  "];
Print["                          2"];
Print["            t (1 + 4 t + t )"];
Print["   Out[4]=  ----------------"];
Print["                       4"];
Print["               (-1 + t)"];
Print["    "];
Print["   In[5]:= formula"];
Print[" "]; 
Print["   Binomial Formula : C[n, 3] + 4 C[1 + n, 3] + C[2 + n, 3]"];
Print[" "]; 
Print["                       3"];
Print["   Power Formula    : n"];
Print[" "]; 
Print["   In[6]:= formulaN[5]"];
Print["   Out[6]= 125"];
Print[" "]; 
Print["-----------"]; 
Print["   Type \"example3\" for another example."]; 
         )
example3:= ( Print[" "];
Print["EXAMPLE III:"]; Print["------------ "]; 
Print["  "];
Print["   In[2]:= a = {{1, -1, 2}, {1, 0, 3}};  b = {1, 2};  c = {1, -2};"]; 
Print["  "];
Print["   In[3]:= DiophantineGF[a, b, c]"]; Print[" "];
Print["             3       2            3 "];
Print["            t  (1 + t ) (1 + t - t )"];
Print["   Out[3]=  ------------------------"];
Print["                     2           2  "];
Print["             (-1 + t)  (1 + t + t ) "];
Print[" "];
Print["   In[4]:= Series[%, {t, 0, 9}]"]; Print[" "];
Print["            3      4      5      6      7      8      9      10"];
Print["   Out[4]= t  + 2 t  + 3 t  + 4 t  + 5 t  + 5 t  + 6 t  + O[t]"];
Print[" "];
Print["   In[5]:= formula"];
Print["    "];
Print["   Binomial Formula : "];
Print["    "];
Print["          -10 + n         -8 + n                        -7 + n"];
Print["   >   (C[-------, 0] - C[------, 0] + C[-8 + n, 1] - C[------, 0] - "];
Print["             3              3                             3"];
Print["    "];
Print["                             -6 + n                        -4 + n"];
Print["   >      2 C[-7 + n, 1] - C[------, 0] - C[-5 + n, 1] + C[------, 0] + "];
Print["                               3                             3"];
Print["    "];
Print["                           -3 + n"];
Print["   >      C[-4 + n, 1] + C[------, 0] + C[-3 + n, 1] + 2 C[-2 + n, 1]) / 3"];
Print["                             3"];
Print["    "];
Print["   In[6]:= formulaN[7]"];
Print["   Out[6]= 5"];
Print["   "];
Print["   In[7]:= formulaN[9]"];
Print["   Out[7]= 6"];
Print["   "];
Print["   In[8]:= formulaN[10^20 - 3]"];
Print["   Out[8]= 66666666666666666665"];
Print["   "];
Print["   In[9]:= formulaN[10^20 - 2]"];
Print["   Out[9]= 66666666666666666665"];
Print["   "];
Print["   In[10]:= formulaN[10^20 - 1]"];
Print["   Out[10]= 66666666666666666666"];


        )
recurse[a_, sign_] := 
  Module[{a1, a2, a3, u, v, dirtyrow}, 
  If[ !prunetreeQ[a],

    dirtyrow = leafnode[a]; 
    If[dirtyrow == r + 2, updategf[a, sign], 

            (u = uindex[a, dirtyrow]; 
             v = vindex[a, dirtyrow]; 

               If[checkinfiniteQ[a, u, v], Abort[ ] ];

                  a1 = addcolumn[a, u, v]; 
                  a2 = addcolumn[a, v, u]; 
                  a3 = zapcolumn[a1, v]; 

                  recurse[a1, sign]; 
                  recurse[a2, sign]; 
                  recurse[a3, -sign]; )]
    ] ]

updategf[a_, sign_] := 
(
     tpart = mixedpart = 1;
 
     For[j = 1, j <= Length[a[[1]]], j++,
       If[zerocolumnQ[a, j], tpart = tpart/(1 - t^a[[1]][[j]]),
 
          ls = Product[lambda[[k - 1]]^a[[k]][[j]], {k, 2, r + 1}];
          mixedmono = t^a[[1]][[j]]*ls;
 
          maxexp = truncateexp[a, j];
          mixedpart = mixedpart*Sum[mixedmono^k, {k, 0, maxexp}];
         ];
        ];
 
     mixedpart = mixedpart*fudgefactor;
          For[k = 1, k <= r, k++, mixedpart =
            Coefficient[mixedpart, lambda[[k]], 0]
             ];
     gf = gf + tpart*mixedpart*sign;
  (* If the form of how the gf is accumulated needs to be controlled,
     then this is the line where the changes can be made.            *)
);
  (* updategf: updates the gf for a leaf matrix a. The 
     contribution of a column with d in the first row and 
     zeros below is 1/(1-t^d). The product of such terms over 
     the whole matrix is tpart.  Otherwise raise l[i] to the 
     power of the i-th entry in the row, multiply these out, 
     call it ls. Note the indexing of l[] and the rows are 
     shifted by one, since a is augmented. If W = mixedpart 
     denotes this product multiplied further by t^d (which 
     may or may not be 1), then mixedmono is 
     1+ W + .. + W^maxexp. The product of these terms over the 
     whole matrix is  mixedpart. After that adjust mixedpart 
     by fudgefactor.  The final For loop discards all the terms 
     in mixedpart that contains some power of a l[], and leaves 
     a polynomial in t. This polynomial multiplied by tpart is 
     the contribution of this leaf matrix to the gf.                 *)


formula :=
  (* Compute the formula as a function of n for from the expansion
     of the generating function gf. C[x,k] stands for Binomial[x,k]
     if x is an integer, and zero otherwise.                         *)
 Module[
 {i,j, gfden, gfnum, pi, ei, temp, partnum, expn, maxexpn },
    gfden = Denominator[gf];
    gfnum = Numerator[gf];

    binomialformula = 0; powerformula = 0; maxexpn = 0;
    factors = Complement[FactorList[gfden], FactorList[1]];
 
  (* This trick is to fix an incompatibility between versions 2.2 and 
     3.0 of Mathematica. In V2.2 FactorList[(-1 + t)^2]={{-1 + t, 2}} 
     (what is used here) whereas in V3.0 FactorList[(-1 + t)^2] =
     {{1, 1}, {-1 + t, 2}}. The Complement function removes {1,1} if 
     it is produced in FactorList[gfden] if running V3.0.           *)
 
    nooffactors= Length[ factors ];

    partfrac = Apart[ 1/gfden ];

    For[i=1, i <= nooffactors, i++,
      pi =  factors[[i]][[1]];
      ei =  factors[[i]][[2]];
  (* i-th factor of the denominator and its exponent                *)

    temp = Select[partfrac, PolynomialQ[ Factor[# pi^ei], t]&];

  (* This function picks all terms in the partial fraction
     expansion of 1/gfden where the denominator is a divisor of 
     the i-th factor in gfden. This factor is pi^ei, where
     pi(t) is a cyclotomic polynomial                               *)

    partnum = gfnum * pi^ei * temp;

  (* partnum is the numerator of 1/pi^ei in the expansion 
     of the whole generating function. Now multiply by an
     appropriate factor to make sure that it is of the 
     form 1/(1-t^di)^ei.                                            *) 

    For[expn= 1, PolynomialRemainder[1-t^expn, pi, t] =!= 0, expn++ ];
    If[ expn > maxexpn, maxexpn = expn];
    
    partnum =  Factor[ partnum * ((1-t^expn)/pi)^ei ];
    partnum =  Expand[ partnum];
  (* This last partnum is the actual numerator polynomial. 
     The corresponding denominator is (1-t^expn)^ei .               *)

    numdegree = Exponent[ partnum, t];

    For[j=0, j <= numdegree, j++,
      binomialformula = binomialformula + 
      Coefficient[partnum, t, j]*C[ei + (n -j)/expn -1, ei-1];
      powerformula = powerformula + 
      Coefficient[partnum, t, j]*Binomial[ei + (n -j)/expn -1, ei-1]; ];

      ];
    binomialformula = Factor[binomialformula];
    Print[" "];Print["Binomial Formula : ", binomialformula];
    If[maxexpn == 1, powerformula = Factor[powerformula];
      Print[" "]; Print["Power Formula    : ", powerformula]; 
      ]
     ]

B[x_, y_] := If[ Positive[x] && IntegerQ[x], Binomial[x, y], 0];

formulaN[v_] :=
(substitute = binomialformula /. {n -> v, C -> B};
Simplify[substitute]);

Print["Loaded. Type \"help\" for instructions, \"example\" for examples."];