=========================
/* This function computes the PERIODIC contribution to Hilbert series of an
isolated orbifold point 1/r(a,b..). It assumes that the univariate polynomial
ring Q[t] is defined and that the input list L = [a1,a2,.. ] consists of
coprime elements. */
  
function Pc(r,L)
   if &or[GreatestCommonDivisor([r,a]) ne 1 : a in L] // other mugtrap
      then error "Error: Not Coprime";
   end if;
A := (t^r-1) div (t-1);
B := t*&*[1-t^i : i in L];
h_throwaway, al_throwaway, be := XGCD(A,B);
  return r*t*be;
end function;

/* This function computes the GROWING contribution to Hilbert series of an
isolated orbifold point 1/r(a,b..). It assumes that the univariate polynomial
ring Q[t] is defined and that the input list L = [a1,a2,.. ] consists of
coprime elements. */

function Pd(r,L)
n := #L; f := (1-t)^(n+1)*Pc(r,L);
return
&+ [-Coefficient(f,k)*t^i where k is Min(i,n+1-i) : i in [1..n] ]
 / (1-t)^(n+1);
end function;

----------

Exa 1: for a curve C, k=2, A = integral half of KC + 8/17P,
the fractional part is
  1/17 times 1/(1-t^17) times Pc(17, [15]),
where Pc(17,[15]) is the above function:

-9*t^16 - t^15 - 10*t^14 - 2*t^13 - 11*t^12 - 3*t^11 - 12*t^10 - 4*t^9 - 13*t^8
    - 5*t^7 - 14*t^6 - 6*t^5 - 15*t^4 - 7*t^3 - 16*t^2 - 8*t

There is a unique thing I can add to Pc to give it Gorenstein symmetry without
changing P1, P2, namely 8*t/(1-t)^2 times trivial factor. To see this,
multiply by (1-t), to get (1-t)*Pc(17, [15]) 

9*t^17 - 8*t^16 + 9*t^15 - 8*t^14 + 9*t^13 - 8*t^12 + 9*t^11 - 8*t^10 + 9*t^9 -
    8*t^8 + 9*t^7 - 8*t^6 + 9*t^5 - 8*t^4 + 9*t^3 - 8*t^2 - 8*t

Adding 8t+8t^2 would give a symmetric numerator. Instead I add 8t/(1-t):
to give

Pd := 8*t/(1-t)^2;

1/17*(1/(1-t^17)*Pc(17,[15])+ Pd); 

(t^17 + t^15 + t^13 + t^11 + t^9 + t^7 + t^5 + t^3)/(t^18 - t^17 - t + 1)

---------

Exa 2: for a curve C, k = 3, A = integral 1/3 KC + 6/19P

Pc(19,[16]);

-13*t^18 - 7*t^17 - t^16 - 14*t^15 - 8*t^14 - 2*t^13 - 15*t^12 - 9*t^11 - 3*t^10
    - 16*t^9 - 10*t^8 - 4*t^7 - 17*t^6 - 11*t^5 - 5*t^4 - 18*t^3 - 12*t^2 - 6*t

To see Gor symmetry
(1-t)^2*Pc(19,[16]);

-13*t^20 + 19*t^19 - 19*t^17 + 19*t^16 - 19*t^14 + 19*t^13 - 19*t^11 + 19*t^10 -
    19*t^8 + 19*t^7 - 19*t^5 + 19*t^4 - 6*t

Pd := 6*t/(1-t)^2;

1/19*(1/(1-t^19)*Pc(19,[16])+ Pd); 

(t^19 + t^16 + t^13 + t^10 + t^7 + t^4)/(t^20 - t^19 - t + 1)


? Question: the Pc is now automatic enough to be a computer function. How did
I know to put in 8*t/(1-t)^2 in Exa 1 and 6*/(1-t)^2 in Exa 2?

---------

Now do a surface with 1/7(2,4) and A = KS

Pc(7,[2,4]);
 // -t^6 - 3*t^5 + t^4 - 3*t^3 - t^2

Pd := t^2/(1-t)^3; // just enough to kill the -t^2
1/7*(1/(1-t^7)*Pc(7,[2,4]) + Pd);

 (-t^6 + t^5 - t^4)/(t^9 - 2*t^8 + t^7 - t^2 + 2*t - 1)
 S!$1; s^4 + s^5 + 2*s^6 + 3*s^7 + 4*s^8 + 5*s^9 + 6*s^10 + ..

1/7(3,3)

Pc(7,[3,3]);
 3*t^6 + 2*t^5 - 3*t^4 + 2*t^3 + 3*t^2

Pd := -3*t^2/(1-t)^3;

 1/7*(1/(1-t^7)*Pc(7,[3,3]) + Pd);
 
 (t^7 + t^6 - t^5 + t^4 + t^3)/(t^9 - 2*t^8 + t^7 - t^2 + 2*t - 1)

Pc(7,[2,3]);
 // -2*t^6 - 3*t^5 - 3*t^4 - 2*t^3 - 4*t

Pd := (4*t-12*t^2)/(1-t)^3;
 1/7*(1/(1-t^7)*Pc(7,[2,3]) + Pd);
 // (2*t^8 + t^7 + t^6 + t^5 + t^4 + 2*t^3)/(t^9 - 2*t^8 + t^7 - t^2 + 2*t - 1)

Pc(7,[4,5]);
 Pd:= (2*t-3*t^2)/(1-t)^3;
 1/7*(1/(1-t^7)*Pc(7,[4,5]) + Pd); // (t^8 - t^7 + t^6)/(t^9 - 2*t^8 + t^7 - t^2 + 2*t - 1)

Pd:= &+[-Coefficient((1-t)^3*Pc(7,[4,5]),i)*t^i : i in [1..2]]/(1-t)^3;

r := 13; L := [1,2,9];
f := Pc(r,L); n := #L;
Pd:= &+[-Coefficient((1-t)^(n+1)*f,i)*t^i : i in [1..n]]/(1-t)^(n+1);
1/r*(1/(1-t^r)*Pc(r,L) + Pd);

r := 13; L := [1,2,9];
1/r*(1/(1-t^r)*Pc(r,L) + Pd(r,L));

r := 7; L := [1,2,4];
Pa := (1 + t^2)/(1 - t)^2;
P := Pa + 1/r*(1/(1-t^r)*Pc(r,L) + Pd(r,L));

P*&*[1-t^i : i in [1,1,2,7]];  // 1 + t^11
 // so Y(22) in PP(1,1,2,7,11)

r := 9; L := [1,2,2,4]; // I'm looking for a CY 4-fold
a1 := 0; a2 := 0;
Pa := (1 + a1*t + a2*t^2 + a2*t^3 + a1*t^4 + t^5)/(1-t)^5;
P := Pa + 1/r*(1/(1-t^r)*Pc(r,L) + Pd(r,L));
P *&*[1-t^i : i in [1,1,1,2,9]];

a1 := 0; a2 := 0;
Pa := (1 + a1*t + a2*t^2 + a2*t^3 + a1*t^4 + t^5)/(1-t)^4;

g := (-8/3*t)/(1-t)^5;  

g := (-8/3*t+41/3*t^2)/(1-t)^5;

g := (-8/3*t+41/3*t^2+41/3*t^3-8/3*t^4)/(1-t)^5;

Now I want to make a CY 4-fold with 1/4(1,1,3,3) and 1/5(1,2,3,4)

r:=4;L:=[1,1,3,3];
Pc(r,L);  // -5/8*t^3 - t^2 - 5/8*t

Pc4 := 1/r*(1/(1-t^r)*Pc(r,L)+ 0);
 g:=5/8*t/(1-t)^5;
 Pc4 := 1/r*(1/(1-t^r)*Pc(r,L)+ g);

 g:=(5/8*t-17/8*t^2-17/8*t^3+5/8*t^4)/(1-t)^5;
 Pc4 := 1/r*(1/(1-t^r)*Pc(r,L)+ g); // (t^3 + t^4 + t^5) / denom

r:=5; L:=[1,2,3,4]; // -t^4 - t^3 - t^2 - t

g:= t/(1-t^5);

 Pc5 := 1/r*(1/(1-t^r)*Pc(r,L)+ 0);

 g:= (t+t^2+t^3+t^4)/(1-t^5);       
 Pc5 := 1/r*(1/(1-t^r)*Pc(r,L)+ g); // 0

Now:
a1 :=-4; a2:= 0;

===================================================

Want a CY 4-fold with 1/4(1,1,3,3) and 1/5(1,2,3,4)?

a1 := -3; a2 := 3;
Pa := (1 + a1*t + a2*t^2 + a2*t^3 + a1*t^4 + t^5)/(1-t)^5;

P := Pa + 1/4*(1/(1-t^4)*Pc(4,[1,1,3,3]) + Pd(4,[1,1,3,3]))
                + 1/5*(1/(1-t^5)*Pc(5,[1,2,3,4]) + Pd(5,[1,2,3,4]));

// that tests the EVEN case of function Pd(r,L)

Want an orbicurve with 1/8(1) and 1/15(2)?
Pa := (1 -2*t + t^2 + t^5 -2*t^6 + t^7)/(1-t)^2;
P := Pa + 1/8*(1/(1-t^8)*Pc(8,[1]) + Pd(8,[1]))
          + 1/15*(1/(1-t^15)*Pc(15,[8]) + Pd(15,[8]));

The curve doesn't work, but the plurigenera are right.

==============

// CY wwith 3 x 1/3(1,1,1) and 1/13(1,3,9)

P1:=1; P2:=1;
Pa := (1 + (P1-3)*t+ (P2-3*P1+3)*t^2 + (P1-3)*t^3 + t^4)/(1-t)^4;

P := Pa + 3*1/3*(1/(1-t^3)*Pc(3,[1,1,1]) + Pd(3,[1,1,1]))
            + 1/13*(1/(1-t^13)*Pc(13,[1,3,9]) + Pd(13,[1,3,9]));

P*&*[1-t^i : i in [1,1,3,9,13]];
 //  -t^27 + 1