function AB(r,L)
K<ep> := CyclotomicField(r);
return  1/r * 1/(1-t^r) *
  &+[ Q! &+[(ep^(-i*j)-1)/&*[1-ep^(-a*j) : a in L]
    : j in [1..r-1] | &and[(a*j) mod r ne 0 : a in L]] * t^i : i in [1..r-1]];
end function;

function denom(L) return &*[1-t^i:i in L]; end function;

function Qorb(r,LL,k)
  if (k + &+[i : i in LL]) mod r ne 0
     then error "Error: Canonical weight not compatible";
  end if;
 n := #LL;
 h := Degree(GCD(1-t^r, denom(LL))); // the degree of GCD(A,B) -- should calc mor
 l := Floor((k+n+1)/2+h);
 // If l < 0 we need a kludge to avoid programming Laurent polynomials properly
 de := Maximum(0,Ceiling(-l/r));
 m := l + de*r;
 A:=(t^r-1) div (t-1);
 B:= &*[(1-t^i) div (1-t) : i in LL];
 H,al_throwaway,be:=XGCD(A,t^m*B);
 return t^m*be/(H*(1-t)^n*(1-t^r)*t^(de*r));
end function;

> X:=1/denom([2,5,15]); X;
-1/(t^22 - t^20 - t^17 + t^15 - t^7 + t^5 + t^2 - 1)
> X-Qorb(2,[1,1],-22)-Qorb(15,[2,5],-22)+Qorb(5,[2],-17)/(1-t^5);
(-t^3 - 1)/(t^16 - 2*t^15 + t^14 - t^11 + 2*t^10 - t^9)

(1/t^9 + 1/t^6)/denom([1,1,5]);

a:=3;b:=3*5;
1/denom([2,a,b])-Qorb(2,[1,1],k)-Qorb(a,[2,b],k)-Qorb(b,[2,a],k)
 -Qorb(3,[2],k+3)/(1-t^3) + (1/t^8 + 1/t^7)/denom([1,1,3]) where k is -2-a-b;

a:=3;b:=3*3;
1/denom([2,a,b])-Qorb(2,[1,1],k)-Qorb(a,[2,b],k)-Qorb(b,[2,a],k)
 + Qorb(3,[2],k+3)/(1-t^3) - (1+t)/t^5/denom([1,1,3]) where k is -2-a-b;

c:=1;a:=3;b:=3*(2*c+1);
1/denom([2,a,b])-Qorb(2,[1,1],k)-Qorb(a,[2,b],k)-Qorb(b,[2,a],k)
 - (-1)^c*( Qorb(3,[2],k+3)/(1-t^3) - (1+t)/t^5/denom([1,1,3])) where k is -2-a-b;

c:=2;a:=3;b:=3*(2*c+1);
1/denom([2,a,b])-Qorb(2,[1,1],k)-Qorb(a,[2,b],k)-Qorb(b,[2,a],k)
 -Qorb(3,[2],k+3)/(1-t^3) + (1/t^8 + 1/t^7)/denom([1,1,3]) where k is -2-a-b;

c:=3;a:=3;b:=3*(2*c+1);
1/denom([2,a,b])-Qorb(2,[1,1],k)-Qorb(a,[2,b],k)-Qorb(b,[2,a],k)
 + Qorb(3,[2],k+3)/(1-t^3) - (1+t)/t^(-(k div 2)-2)/denom([1,1,3]) where k is -2-a-b;

c:=4;a:=3;b:=3*(2*c+1);
1/denom([2,a,b])-Qorb(2,[1,1],k)-Qorb(a,[2,b],k)-Qorb(b,[2,a],k)
 - Qorb(3,[2],k+3)/(1-t^3) + (1+t)/t^(-(k div 2)-2)/denom([1,1,3]) where k is -2-a-b;

[1/denom([2,a,b])-Qorb(2,[1,1],k)-Qorb(a,[2,b],k)-Qorb(b,[2,a],k)
  -(-1)^c*(Qorb(3,[2],k+3)/(1-t^3) - (1+t)/t^(-(k div 2)-2)/denom([1,1,3]))
  where k is -2-a-b where a is 3 where b is 3*(2*c+1)
  : c in [1..50]];

&and[(1/denom([2,a,b])-Qorb(2,[1,1],k)-Qorb(a,[2,b],k)-Qorb(b,[2,a],k)
  -(-1)^c*(Qorb(3,[2],k+3)/(1-t^3) - (1+t)/t^(-(k div 2)-2)/denom([1,1,3]))) eq 0
  where k is -2-a-b where a is 3 where b is 3*(2*c+1)
  : c in [1..50]]; // true

&and[(1/denom([1,a*r,b*r]) - Qorb(a*r,[1,b*r],k) - Qorb(b*r,[1,a*r],k)
  where k is -1-a*r-b*r) eq 0
  : a in [2..15], b in [2..30], r in [4,7,8,9,10] | a ne b]; // true

&and[(1/denom([2,a,b])-Qorb(2,[1,1],k)-Qorb(a,[2,b],k)-Qorb(b,[2,a],k)
  -(-1)^c*(Qorb(5,[2],k+5)/(1-t^5) - (1+t^3)/t^(-(k div 2)-2)/denom([1,1,5]))) eq 0
  where k is -2-a-b where a is 5 where b is 5*(2*c+1)
  : c in [2..50]];

&and[(1/denom([2,a,b])-Qorb(2,[1,1],k)-Qorb(a,[2,b],k)-Qorb(b,[2,a],k)
  -(-1)^c*(Qorb(7,[2],k+7)/(1-t^7) - (1+t^2+t^3+t^5)/t^(-(k div 2)-2)/denom([1,1,7]))
  ) eq 0
  where k is -2-a-b where a is 7 where b is 7*(2*c+1)
  : c in [6..25]]; // true

&and[(1/denom([2,a,b])-Qorb(2,[1,1],k)-Qorb(a,[2,b],k)-Qorb(b,[2,a],k)
  -(-1)^c*(Qorb(11,[2],k+11)/(1-t^11) - (1+t^2+t^4+t^5+t^7+t^9)
   /t^(-(k div 2)-2)/denom([1,1,11]))
   ) eq 0
  where k is -2-a-b where a is 11 where b is 11*(2*c+1)
  : c in [101..115]]; // true

&and[(1/denom([2,a,b])-Qorb(2,[1,1],k)-Qorb(a,[2,b],k)-Qorb(b,[2,a],k)
  -(-1)^c*(Qorb(23,[2],k+23)/(1-t^23) - (1+t^2+t^4+t^6+t^8+t^10+t^11+t^13+t^15+t^17+t^19+t^21)
   /t^(-(k div 2)-2)/denom([1,1,23]))
    ) eq 0
  where k is -2-a-b where a is 23 where b is 23*(2*c+1)
  : c in [6..10]]; // true

1/denom([3,4,4])-Qorb(3,[4,4],k)
  -Qorb(4,[3],k+4)/(1-t^4) + (1+t)/t^3/denom([1,1,4])
  where k is -3-4-4;

1/denom([3,4,8])-Qorb(3,[4,8],k)-Qorb(8,[3,4],k)
  -Qorb(4,[3],k+4)/(1-t^4) + (1+t)/t^3/denom([1,1,4])
  where k is -3-4-8;