symbolic - Improving the performance of a package for working with rational functions

As Mathematica gets slow for large symbolic calculations, the cost of putting terms over a common denominator (Together), in particular, gets too high. It occurred to me that, if one has a small number of variables, it should be much faster if we represent our expressions by arrays of their coefficients.

I have implemented a small package that does this with an improvement of a factor of 10 in both speed and memory (in putting Together a sum of rational terms, more or less independent of the size of the sum), but it seems to me that I should be able to do better. I post my package here with the hope that someone might point out glaring performance bottlenecks.

First, the public initializations, with their usability comments.

AA::usage = "DATA STRUCTURE AA[num,den] that stores polynomials as their tensor coefficients of a basis, multiplication is the tensor product"
PolytoAA::usage = "PolytoAA[poly_,vars_List] returns an AA[num_SparseArray,1]"
AAtoPoly::usage = "AAtoPoly[aa_AA,vars_] gives back a rational function of the polynomials num/den"

The actual code begins here, PolProd takes care of multiplying polynomials

SetAttributes[PolProd, Orderless]
PolProd[i__Integer] := Times @ i;

PolProd[sa__SparseArray] := Outer[Times, sa];
  (*multiplying 2 polynomials is the tensor product in their vector space*)
PolProd[i__Integer, sa__SparseArray] := Times[i] Outer[Times, sa];

(*It would be cleaner and more object oriented to call two instances of polprod 
  on i and sa but it would force more pattern matching so it would be less efficient*)

This code performs the arithmetic operations on the AA[num,den] data structure

AA /: i_Integer*AA[num_, den_] := AA[i*num, den]
 (*absorb all integer into AA, this means that the sum rule has to do less 

   pattern matching and is more efficient, overall more efficient?*)
AA /: AA[a_, a_] := 1
AA /: AA[0, _] := 0
AA /: AA[num1_, den1_]*AA[num2_, den2_] := 
  AA @@ (PolProd[#1, #2]& @@@ {{num1, num2}, {den1, den2}})
AA /: AA[num1_, den_] + AA[num2_, den_] := 
  AA[num1 + num2, den]
  (*can't check for a SparseArray cos it might be an Integer, 
    by dimensional analysis we should get two arrays here though*)
AA /: AA[num1_, den1_] + AA[num2_, den2_] := 

  AA[PolProd[num1, den2] + PolProd[num2, den1], PolProd[den1, den2]]
  (*together rule, this a recursive implementation that doesn't treat the 
    sum of many AA all at once, this is intentional cause it minimises the 
    number of array outer products that are calculated and I believe it's 
    the more efficient but maybe not*)
AA /: aa_AA^n_Integer := 
  If[n < 0, Reverse @ #, #]& @ (PolProd[Sequence @@ ConstantArray[#, Abs @ n]]& /@ aa)
  (*this takes care of both efficiency of recursively pattern matching the 
    product rule on large n and dividing by AA*)
(*NOTICE THE WEIRD SYNTAX/BUG If[n<0, Reverse @ #, #]& @ THERE ARE TWO INTANCES OF @*)

I also post the the functions that go between AA and the polynomials for the sake of clarity and completeness

AAtoPoly[aa_AA,vars_] := #1/#2& @@ (aa /. l_SparseArray :> 
  Dot[l, Sequence @@ ConstantArray[vars, Depth[l] - 1 (*gives rank of the tensor*)]])
PolytoAA[expr_^n_, vars_List] := PolytoAA[expr, vars]^n
PolytoAA[expr : (_Times | _Plus), vars_List] := PolytoAA[#, vars]& /@ expr
PolytoAA[expr_,vars_List] /; PolynomialQ[expr, vars] := 
  AA[Last[CoefficientArrays[expr, vars]](*num*), 1 (*den*)]

BENCHMARKS

I define some random, small polynomials

vars = a /@ Range[10]
pols = Table[
   a[RandomInteger[{1, 10}]] + a[RandomInteger[{1, 10}]], {20}] // 
  DeleteDuplicates
% // Length

{a[1], a[2], a[3], a[4], a[5], a[6], a[7], a[8], a[9], a[10]}

{a[1] + a[7], a[4] + a[6], a[2] + a[4], a[5] + a[9], 2 a[8], a[5] + a[6], a[2] + a[9], a[4] + a[7], a[4] + a[9], a[7] + a[9], a[3] + a[7], a[2] + a[6], a[1] + a[9], a[2] + a[3], a[7] + a[10], 2 a[2], a[3] + a[9]}

17

First I put them Together using my code

aa = PolytoAA[#, vars] & /@ pols;
togaa = Sum[aa[[i]]/aa[[i + 1]], {i, Length@pols - 1}] // 
  AbsoluteTiming
% // ByteCount

{0.8960513,AA[SparseArray[<274432>,{10,10,10,10,10,10,10,10,10,10,10,10,10,10,10,10}],SparseArray[<16384>,{10,10,10,10,10,10,10,10,10,10,10,10,10,10,10,10}]]}

18614128

and now without my implementation

togpols = (Sum[pols[[i]]/pols[[i + 1]], {i, Length@pols - 1}] // 
    Together) // AbsoluteTiming//
% // ByteCount

{8.2054693,...}

34567480

I check that both give the same result

togpols2 = AAtoPoly[togaa // Last, vars](*//Simplify*);
Last[togpols] - togpols2 // Together

0