/* ************************* MainTheorem *************************** */

/* The following script determines if the surface w^2 = z^3 + p^3x^6 + p^3y^6 over the rationals gives a counterexample to Weak Approximation. By varying the parameter LastPrime you can search for primes p that give rise to counterexamples. */

LastPrime := 100;
ListOfPrimes := [];
for i in [5..LastPrime] do
if IsPrime(i) then

p := i;
Q := Rationals();
P1<t1> := PolynomialRing(Q);
f1 := t1^2 - t1 + 1;
N1<t> := NumberField(f1);
P2<t2> := PolynomialRing(N1);
f2 := t2^3 - 2;
N2<s> := NumberField(f2);
P3<t3> := PolynomialRing(N2);
f3 := t3^2 - p;
N3<A> := NumberField(f3); 

R<w,z,y,x> := PolynomialRing(N3,[3,2,1,1]);
Q1<w,z,y,x> := FieldOfFractions(R);

WPSpace := Proj(R);

/* L below stands for \Gamma in the paper, e.g., L6' is the exceptional curve
\Gamma_6' in the paper. Also, SL6 means sigma applied to L6, for example,
and TL6 means tau applied to L6. */

/* The funtion g below vanishes on L8, L8', SL6 and SL6' */
/* The function g2 below vanishes on L8, L8', TL6, TL6'. */

// P<w,z,y,x,s,t,p,A> := PolynomialRing(Rationals(),8);

g := (z - A^2*((s^2*t + 2*s*t + 2*t)*x^2 + (2*s^2 + 3*s + 4)*x*y + (-s^2*t + s^2 - 
2*s*t + 2*s - 2*t + 2)*y^2))*(z - A^2*((s^2*t + 2*s*t + 2*t)*x^2 + (-2*s^2*t - 
3*s*t - 4*t)*x*y + (s^2*t + 2*s*t + 2*t)*y^2));

g2 := (z - A^2*((s^2*t + 2*s*t + 2*t)*x^2 + (2*s^2 + 3*s + 4)*x*y + (-s^2*t + 
s^2 - 2*s*t + 2*s - 2*t + 2)*y^2))*(z - A^2*((s^2*t - 2*s - 2*t + 2)*x^2 + 
(-2*s^2*t + 3*s + 4*t - 4)*x*y + (s^2*t - 2*s - 2*t + 2)*y^2));

// The following product vanishes on P' + (rhobar)P' cf section 7 of the paper.

hdenom := 

(z - A^2*((s^2*t + 2*s*t + 2*t)*x^2 + (2*s^2 + 3*s + 4)*x*y + (-s^2*t + s^2 - 2*s*t + 2*s - 2*t + 
    2)*y^2)) // this bit vanishes on L8 and L8'
*(z - A^2*((-s^2*t + s^2 - 2*s + 2*t)*x^2 + (-2*s^2*t + 3*s*t - 3*s + 4)*x*y + (-s^2 + 2*s*t - 2*t + 
    2)*y^2)) // this bit vanishes on s^2L8 and s^2L8'
*(z - A^2*((-s^2*t + s^2 - 2*s*t + 2*s - 2*t + 2)*x^2 + (2*s^2 + 3*s + 4)*x*y + (s^2*t + 2*s*t + 
    2*t)*y^2))  // this bit vanishes on tL8 and tL8'
*(z - A^2*((-s^2 - 2*s*t + 2*s + 2*t)*x^2 + (2*s^2*t - 2*s^2 - 3*s*t + 4)*x*y + (s^2*t - 2*s - 2*t + 
    2)*y^2)) // this bit vanishes on sL8 and sL8' 
*(z - A^2*((s^2*t - 2*s - 2*t + 2)*x^2 + (2*s^2*t - 2*s^2 - 3*s*t + 4)*x*y + (-s^2 - 2*s*t + 2*s + 
    2*t)*y^2))  // this bit vanishes on S2TL8 and S2TL8'
*(z - A^2*((-s^2 + 2*s*t - 2*t + 2)*x^2 + (-2*s^2*t + 3*s*t - 3*s + 4)*x*y + (-s^2*t + s^2 - 2*s + 
    2*t)*y^2)); // this bit vanishes on STL8 and STL8'

// S := Scheme(WPSpace, [w^2 - 2*x^6 - 1*y^6 - z^3, g]);
// PrimaryComponents(S);

Gset0 := ideal< R | w^2 - A^6*x^6 - A^6*y^6 - z^3 >;

Gset1 := ideal< R | z - A^2*((-s^2*t + s^2 - 2*s + 2*t)*x^2 + (2*s^2*t - 2*s^2 + 
3*s - 4*t)*x*y + (-s^2*t + s^2 - 2*s + 2*t)*y^2), w - A^3*((2*s^2*t - 4*s^2 + 2*s*t 
+ 2*s - 6*t + 3)*x^3 + (-5*s^2*t + 10*s^2 - 6*s*t - 6*s + 16*t - 8)*x^2*y + 
  (5*s^2*t - 10*s^2 + 6*s*t + 6*s - 16*t + 8)*x*y^2 + (-2*s^2*t + 4*s^2 - 2*s*t - 
  2*s + 6*t - 3)*y^3) >; // L6

Gset2 := ideal< R | z - A^2*((-s^2 - 2*s*t + 2*s + 2*t)*x^2 + (2*s^2*t - 2*s^2 - 
3*s*t + 4)*x*y + (s^2*t - 2*s - 2*t + 2)*y^2), w - A^3*((2*s^2*t + 2*s^2 + 2*s*t - 
4*s - 6*t + 3)*x^3 + (-5*s^2*t + 10*s^2 + 12*s*t - 6*s - 8*t - 8)*x^2*y + 
(-10*s^2*t + 5*s^2 + 6*s*t + 6*s + 8*t - 16)*x*y^2 + (-2*s^2*t - 2*s^2 - 2*s*t + 
4*s + 6*t - 3)*y^3) >; // SL8

Gset3 := ideal< R | z - A^2*((s^2*t + 2*s*t + 2*t)*x^2 + (2*s^2 + 3*s + 4)*x*y + 
(-s^2*t + s^2 - 2*s*t + 2*s - 2*t + 2)*y^2), w - A^3*((4*s^2*t - 2*s^2 + 4*s*t - 
2*s + 6*t - 3)*x^3 + (5*s^2*t + 5*s^2 + 6*s*t + 6*s + 8*t + 8)*x^2*y + (-5*s^2*t + 
10*s^2 - 6*s*t + 12*s - 8*t + 16)*x*y^2 + (-4*s^2*t + 2*s^2 - 4*s*t + 2*s - 6*t + 
3)*y^3) >; // L8'

Gset4 := ideal< R | z - A^2*((s^2*t + 2*s*t + 2*t)*x^2 + (-2*s^2*t - 3*s*t - 
4*t)*x*y + (s^2*t + 2*s*t + 2*t)*y^2), w - A^3*((4*s^2*t - 2*s^2 + 4*s*t - 2*s + 
6*t - 3)*x^3 + (-10*s^2*t + 5*s^2 - 12*s*t + 6*s - 16*t + 8)*x^2*y + (10*s^2*t - 
5*s^2 + 12*s*t - 6*s + 16*t - 8)*x*y^2 + (-4*s^2*t + 2*s^2 - 4*s*t + 2*s - 6*t + 
3)*y^3) >; // SL6'

Gset5 := ideal< R | z - A^2*((-s^2*t + s^2 - 2*s*t + 2*s - 2*t + 2)*x^2 + (2*s^2 
+ 3*s + 4)*x*y + (s^2*t + 2*s*t + 2*t)*y^2), w - A^3*((4*s^2*t - 2*s^2 + 4*s*t - 
  2*s + 6*t - 3)*x^3 + (5*s^2*t - 10*s^2 + 6*s*t - 12*s + 8*t - 16)*x^2*y + 
  (-5*s^2*t - 5*s^2 - 6*s*t - 6*s - 8*t - 8)*x*y^2 + (-4*s^2*t + 2*s^2 - 4*s*t + 
  2*s - 6*t + 3)*y^3) >; // TL8

Gset6 := ideal< R | z - A^2*((s^2*t - 2*s - 2*t + 2)*x^2 + (-2*s^2*t + 3*s + 4*t 
- 4)*x*y + (s^2*t - 2*s - 2*t + 2)*y^2), w - A^3*((2*s^2*t + 2*s^2 + 2*s*t - 4*s - 
  6*t + 3)*x^3 + (-5*s^2*t - 5*s^2 - 6*s*t + 12*s + 16*t - 8)*x^2*y + (5*s^2*t + 
  5*s^2 + 6*s*t - 12*s - 16*t + 8)*x*y^2 + (-2*s^2*t - 2*s^2 - 2*s*t + 4*s + 6*t - 
  3)*y^3) >; // TL6'

/*
I1 := Gset1 meet Gset2 meet Gset3 meet Gset4; 
I2 := Gset1 meet Gset3 meet Gset5 meet Gset6;

gb1 := GroebnerBasis(Saturation(I1));
gb2 := GroebnerBasis(Saturation(I2));
// print gb;

// The third element of GroebnerBasis(Saturation(I1)) is the function f 
// below. We have a normalization constant in front of it for purposes of descent
// The divisor of the function f/g is L6 + SL8 - L8 - SL6

f := (1/(1/(6*p)*(-2*t + 1)*s*A))*(gb1[3]);
// print f;
print gb1[3];

f2 := (1/(-1/(3*p)*s*t*A + 1/(6*p)*s*A))*gb2[3];

// The divisor of the function f2/g2 is - L8 - TL6 + L6 + TL8

// S2f is what happens when we apply G1^2 (i.e. sigma^2) to f. Ditto for the other functions 
// below. S is G1 and T is G2 (i.e. T is tau)

S2f := (-t + 2)*s^2*A*w*y + ((2*t - 1)*s^2 + (-t + 2)*s)*A*w*x + z^2 + (p*t*s^2 +
(-2*p*t + 2*p)*s - 2*p)*z*y^2 + (-p*s^2 +
    3*p*t*s - 2*p*t + 2*p)*z*y*x + ((-1*p*t + 1*p)*s^2 - 2*p*s + 2*p*t)*z*x^2 +
(p^2*t*s^2 +
(-2*p^2*t + 2*p^2)*s)*y^4 + (-p^2*s^2 +
    3*p^2*t*s - 2*p^2*t + 2*p^2)*y^3*x + ((-p^2*t + p^2)*s^2 - 2*p^2*s +
2*p^2*t)*y^2*x^2 +
(1*p^2*t*s^2 - 2*p^2)*y*x^3 + (-1*p^2*s^2 -
    1*p^2*t*s)*x^4;

S2g := z^2 + (p*t*s^2 + (-2*p*t + 2*p)*s - 2*p)*z*y^2 + (2*p*s^2 - 3*p*t*s + 4*p*t 
- 4*p)*z*y*x + ((2*p*t - 2*p)*s^2 + 4*p*s - 
    4*p*t)*z*x^2 + (-8*p^2*t*s^2 + (10*p^2*t - 10*p^2)*s + 12*p^2)*y^4 + (14*p^2*s^2 
- 18*p^2*t*s + 22*p^2*t - 22*p^2)*y^3*x + ((17*p^2*t - 
    17*p^2)*s^2 + 22*p^2*s - 28*p^2*t)*y^2*x^2 + (-14*p^2*t*s^2 + (18*p^2*t - 
18*p^2)*s + 22*p^2)*y*x^3 + (8*p^2*s^2 - 10*p^2*t*s + 12*p^2*t - 12*p^2)*x^4;

Sf := (-t - 1)*s^2*A*w*y + ((-t + 2)*s^2 + (2*t - 1)*s)*A*w*x + z^2 + ((-p*t + 
p)*s^2 + 2*p*t*s - 2*p)*z*y^2 + (p*t*s^2 - 3*p*s
    - 2*p*t + 2*p)*z*y*x + (-p*s^2 + (-2*p*t + 2*p)*s + 2*p*t)*z*x^2 + ((-p^2*t + 
p^2)*s^2 
+ 2*p^2*t*s)*y^4 + (p^2*t*s^2 - 3*p^2*s -
    2*p^2*t + 2*p^2)*y^3*x + (-1*p^2*s^2 + (-2*p^2*t + 2*p^2)*s + 2*p^2*t)*y^2*x^2 + 
((-p^2*t + p^2)*s^2 - 2*p^2)*y*x^3 + (p^2*t*s^2 + p^2*s)*x^4;

Sg := z^2 + ((-p*t + p)*s^2 + 2*p*t*s - 2*p)*z*y^2 + (-2*p*t*s^2 + 3*p*s + (4*p*t 
- 4*p))*z*y*x + (2*p*s^2 + (4*p*t - 4*p)*s - 
    4*p*t)*z*x^2 + ((8*p^2*t - 8*p^2)*s^2 - 10*p^2*t*s + 12*p^2)*y^4 + 
(-14*p^2*t*s^2 + 18*p^2*s + (22*p^2*t - 22*p^2))*y^3*x + (17*p^2*s^2 + 
    (22*p^2*t - 22*p^2)*s - 28*p^2*t)*y^2*x^2 + ((14*p^2*t - 14*p^2)*s^2 - 
18*p^2*t*s + 22*p^2)*y*x^3 + (-8*p^2*t*s^2 + 10*p^2*s + (12*p^2*t - 12*p^2))*x^4;

Tf2 := (-2*t + 1)*s^2*A*w*y + (2*t - 1)*s^2*A*w*x + z^2 + (-p*s^2 - 2*p*s - 
2*p)*z*y^2 + (-p*s^2 - 2*p)*z*y*x + (-p*s^2 - 2*p*s -
    2*p)*z*x^2 + (-p^2*s^2 - 2*p^2*s)*y^4 + (-p^2*s^2 - 2*p^2)*y^3*x + (-p^2*s^2 - 
2*p^2*s - 2*p^2)*y^2*x^2 + (-p^2*s^2 - 2*p^2)*y*x^3 + 
    (-p^2*s^2 - 2*p^2*s)*x^4;

Tg2 := z^2 + (-p*s^2 + (-2*p*t + 2*p)*s - 4*p*t)*z*y^2 + (-2*p*t*s^2 - 6*p*s + 4*p*t 
- 4*p)*z*y*x + ((2*p*t - 2*p)*s^2 + 2*p*t*s - 
    2*p)*z*x^2 + (-2*p^2*t*s^2 - 2*p^2*s)*y^4 + ((2*p^2*t - 2*p^2)*s^2 - 
2*p^2)*y^3*x + (-p^2*s^2 + (-2*p^2*t + 2*p^2)*s - 4*p^2*t)*y^2*x^2 + 
    (-2*p^2*t*s^2 - 6*p^2*s + 4*p^2*t - 4*p^2)*y*x^3 + ((2*p^2*t - 2*p^2)*s^2 + 
2*p^2*t*s)*x^4;

Sf2 := (-t - 1)*s^2*A*w*y + (t + 1)*s^2*A*w*x + z^2 + ((-p*t + p)*s^2 + 2*p*t*s - 
2*p)*z*y^2 + ((-p*t + p)*s^2 - 2*p)*z*y*x + 
    ((-p*t + p)*s^2 + 2*p*t*s - 2*p)*z*x^2 + ((-p^2*t + p^2)*s^2 + 2*p^2*t*s)*y^4 + 
((-p^2*t + p^2)*s^2 - 2*p^2)*y^3*x + ((-p^2*t +
    p^2)*s^2 + 2*p^2*t*s - 2*p^2)*y^2*x^2 + ((-p^2*t + p^2)*s^2 - 2*p^2)*y*x^3 + 
((-p^2*t + p^2)*s^2 + 2*p^2*t*s)*x^4;

Sg2 := z^2 + ((-p*t + p)*s^2 + (-2*p*t + 2*p)*s + (4*p*t - 4*p))*z*y^2 + (-2*p*t*s^2 
+ 6*p*t*s - 4*p*t)*z*y*x + (2*p*s^2 - 2*p*s - 
    2*p)*z*x^2 + (-2*p^2*t*s^2 + 2*p^2*t*s)*y^4 + (2*p^2*s^2 - 2*p^2)*y^3*x + 
((-p^2*t + p^2)*s^2 + (-2*p^2*t + 2*p^2)*s + (4*p^2*t - 
    4*p^2))*y^2*x^2 + (-2*p^2*t*s^2 + 6*p^2*t*s - 4*p^2*t)*y*x^3 + (2*p^2*s^2 - 
2*p^2*s)*x^4;

TSf := (t - 2)*s^2*A*w*y + ((t + 1)*s^2 + (-2*t + 1)*s)*A*w*x + z^2 + (p*t*s^2 + 
(-2*p*t + 2*p)*s - 2*p)*z*y^2 + ((-p*t + p)*s^2
    - 3*p*s + 2*p*t)*z*y*x + (-p*s^2 + 2*p*t*s - 2*p*t + 2*p)*z*x^2 + (p^2*t*s^2 + 
(-2*p^2*t + 2*p^2)*s)*y^4 + ((-p^2*t + p^2)*s^2 -
    3*p^2*s + 2*p^2*t)*y^3*x + (-p^2*s^2 + 2*p^2*t*s - 2*p^2*t + 2*p^2)*y^2*x^2 + 
(p^2*t*s^2 - 2*p^2)*y*x^3 + ((-p^2*t + p^2)*s^2 + p^2*s)*x^4;

TSg := z^2 + (p*t*s^2 + (-2*p*t + 2*p)*s - 2*p)*z*y^2 + ((2*p*t - 2*p)*s^2 + 3*p*s - 
4*p*t)*z*y*x + (2*p*s^2 - 4*p*t*s + 4*p*t - 
    4*p)*z*x^2 + (-8*p^2*t*s^2 + (10*p^2*t - 10*p^2)*s + 12*p^2)*y^4 + ((14*p^2*t - 
14*p^2)*s^2 + 18*p^2*s - 22*p^2*t)*y^3*x + (17*p^2*s^2 - 
    22*p^2*t*s + 28*p^2*t - 28*p^2)*y^2*x^2 + (-14*p^2*t*s^2 + (18*p^2*t - 18*p^2)*s 
+ 22*p^2)*y*x^3 + 
((8*p^2*t - 8*p^2)*s^2 + 10*p^2*s - 
    12*p^2*t)*x^4;

Tf := (-2*t + 1)*s^2*A*w*y + ((t - 2)*s^2 + (t - 2)*s)*A*w*x + z^2 + (-p*s^2 - 2*p*s 
- 2*p)*z*y^2 + (p*t*s^2 + 3*p*t*s + 
    2*p*t)*z*y*x + ((-p*t + p)*s^2 + (-2*p*t + 2*p)*s - 2*p*t + 2*p)*z*x^2 + 
(-p^2*s^2 - 2*p^2*s)*y^4 + (p^2*t*s^2 + 3*p^2*t*s + 
    2*p^2*t)*y^3*x + ((-p^2*t + p^2)*s^2 + (-2*p^2*t + 2*p^2)*s - 2*p^2*t + 
2*p^2)*y^2*x^2 + (-p^2*s^2 - 2*p^2)*y*x^3 + (p^2*t*s^2 - 
    p^2*t*s)*x^4;

Tg := z^2 + (-p*s^2 - 2*p*s - 2*p)*z*y^2 + (-2*p*t*s^2 - 3*p*t*s - 4*p*t)*z*y*x + 
((2*p*t 
- 2*p)*s^2 + (4*p*t - 4*p)*s + (4*p*t - 
    4*p))*z*x^2 + (8*p^2*s^2 + 10*p^2*s + 12*p^2)*y^4 + (-14*p^2*t*s^2 - 18*p^2*t*s 
- 22*p^2*t)*y^3*x 
+ ((17*p^2*t - 17*p^2)*s^2 + (22*p^2*t - 
    22*p^2)*s + (28*p^2*t - 28*p^2))*y^2*x^2 + (14*p^2*s^2 + 18*p^2*s + 
22*p^2)*y*x^3 + (-8*p^2*t*s^2 - 
10*p^2*t*s - 12*p^2*t)*x^4;

// The following are extra functions needed to average the rational function 1 
// over the Galois Group

STf2 := (t + 1)*s^2*A*w*y + (-t - 1)*s^2*A*w*x + z^2 + ((-p*t + p)*s^2 + 
2*p*t*s - 2*p)*z*y^2 + ((-p*t + p)*s^2 - 2*p)*z*y*x + ((-p*t + p)*s^2 + 2*p*t*s - 
2*p)*z*x^2 + ((-p^2*t + p^2)*s^2 + 2*p^2*t*s)*y^4 + ((-p^2*t + p^2)*s^2 - 
2*p^2)*y^3*x + 
((-p^2*t + p^2)*s^2 + 2*p^2*t*s - 2*p^2)*y^2*x^2 + ((-p^2*t + p^2)*s^2 - 
2*p^2)*y*x^3 + ((-p^2*t + p^2)*s^2 + 2*p^2*t*s)*x^4;

STg2 := z^2 + ((-p*t + p)*s^2 - 2*p*s - 4*p*t)*z*y^2 + (2*p*s^2 + 6*p*t*s + (4*p*t 
- 4*p))*z*y*x + (-2*p*t*s^2 + (-2*p*t + 2*p)*s - 2*p)*z*x^2 + (2*p^2*s^2 + 
2*p^2*t*s)*y^4 
+ (-2*p^2*t*s^2 - 2*p^2)*y^3*x + ((-p^2*t + p^2)*s^2 - 2*p^2*s - 4*p^2*t)*y^2*x^2 + 
(2*p^2*s^2 + 6*p^2*t*s + (4*p^2*t - 4*p^2))*y*x^3 + (-2*p^2*t*s^2 + (-2*p^2*t + 
2*p^2)*s)*x^4;

// good
TSf2 := (t - 2)*s^2*A*w*y + (-t + 2)*s^2*A*w*x + z^2 + (p*t*s^2 + (-2*p*t + 
2*p)*s - 2*p)*z*y^2 + (p*t*s^2 - 2*p)*z*y*x + (p*t*s^2 + (-2*p*t + 2*p)*s - 
2*p)*z*x^2 + (p^2*t*s^2 + (-2*p^2*t + 2*p^2)*s)*y^4 + (p^2*t*s^2 - 2*p^2)*y^3*x + 
(p^2*t*s^2 + (-2*p^2*t + 2*p^2)*s - 2*p^2)*y^2*x^2 + (p^2*t*s^2 - 2*p^2)*y*x^3 + 
(p^2*t*s^2 + 
(-2*p^2*t + 2*p^2)*s)*x^4;

TSg2 := z^2 + (p*t*s^2 + 2*p*t*s - 4*p*t)*z*y^2 + ((2*p*t - 2*p)*s^2 + (-6*p*t + 
6*p)*s + (4*p*t - 4*p))*z*y*x + (2*p*s^2 - 2*p*s - 2*p)*z*x^2 + ((2*p^2*t - 
2*p^2)*s^2 + 
(-2*p^2*t + 2*p^2)*s)*y^4 + (2*p^2*s^2 - 2*p^2)*y^3*x + (p^2*t*s^2 + 2*p^2*t*s - 
4*p^2*t)*y^2*x^2 + ((2*p^2*t - 2*p^2)*s^2 + (-6*p^2*t + 6*p^2)*s + (4*p^2*t - 
4*p^2))*y*x^3 + (2*p^2*s^2 - 2*p^2*s)*x^4;

// The following is the numerator of the average of the rational section 1 
// over the Galois group

hAverage := 1 + S2f/S2g + Tf2/Tg2 + (S2f/S2g)*(Sf/Sg) + 
(S2f/S2g)*(Sf/Sg)*(STf2/STg2) + (S2f/S2g)*(TSf2/TSg2);
*/

/* Unfortunately, the numerator of hAverage had w^2 and w^3 terms. We had to reduce these terms by hand using the relation 
w^2 = p^3*x^6 + p^3*y^6 + z^3. The result is what we've called MonstrousNumerator below. */

/* MonstrousNumerator and MonstrousDenominator have a factor in common! I know... it's hard to believe. Here it is: */

gcd := 
z^4 + (p*s^2 + (2*p*t - 2*p)*s - 8*p*t)*z^3*y^2 + (-2*p*s^2 + (-3*p*t + 3*p)*s + 
    16*p*t)*z^3*y*x + (p*s^2 + (2*p*t - 2*p)*s - 8*p*t)*z^3*x^2 + (-6*p^2*t*s^2 + 12*p^2*s
    + (12*p^2*t - 12*p^2))*z^2*y^4 + (24*p^2*t*s^2 - 42*p^2*s - 54*p^2*t + 54*p^2)*z^2*y^3*x + 
    (-36*p^2*t*s^2 + 60*p^2*s + (84*p^2*t - 84*p^2))*z^2*y^2*x^2 + (24*p^2*t*s^2 - 42*p^2*s - 
    54*p^2*t + 54*p^2)*z^2*y*x^3 + (-6*p^2*t*s^2 + 12*p^2*s + (12*p^2*t - 12*p^2))*z^2*x^4 + 
    4*p^3*z*y^6 + ((-6*p^3*t + 6*p^3)*s^2 + 12*p^3*t*s - 24*p^3)*z*y^5*x + ((24*p^3*t - 
    24*p^3)*s^2 - 42*p^3*t*s + 72*p^3)*z*y^4*x^2 + ((-36*p^3*t + 36*p^3)*s^2 + 60*p^3*t*s - 
    102*p^3)*z*y^3*x^3 + ((24*p^3*t - 24*p^3)*s^2 - 42*p^3*t*s + 72*p^3)*z*y^2*x^4 + 
    ((-6*p^3*t + 6*p^3)*s^2 + 12*p^3*t*s - 24*p^3)*z*y*x^5 + 4*p^3*z*x^6 + (4*p^4*s^2 + 
    (8*p^4*t - 8*p^4)*s - 8*p^4*t)*y^8 + (-20*p^4*s^2 + (-36*p^4*t + 36*p^4)*s + 
    40*p^4*t)*y^7*x + (52*p^4*s^2 + (92*p^4*t - 92*p^4)*s - 104*p^4*t)*y^6*x^2 + 
    (-90*p^4*s^2 + (-156*p^4*t + 156*p^4)*s + 180*p^4*t)*y^5*x^3 + (108*p^4*s^2 + 
    (186*p^4*t - 186*p^4)*s - 216*p^4*t)*y^4*x^4 + (-90*p^4*s^2 + (-156*p^4*t + 
    156*p^4)*s + 180*p^4*t)*y^3*x^5 + (52*p^4*s^2 + (92*p^4*t - 92*p^4)*s - 
    104*p^4*t)*y^2*x^6 + (-20*p^4*s^2 + (-36*p^4*t + 36*p^4)*s + 40*p^4*t)*y*x^7 + 
    (4*p^4*s^2 + (8*p^4*t - 8*p^4)*s - 8*p^4*t)*x^8;

// numerator := MonstrousNumerator/gcd;
// denominator := MonstrousDenominator/gcd;

numerator := 
(-24*p^6*s^2*t + 24*p^6*s^2 + 24*p^6*s*t - 24*p^6*s)*x^12 + (12*p^6*s^2*t + 12*p^6*s^2 + 24*p^6*s*t - 48*p^6*t)*x^11*y 
+ (72*p^6*s^2*t - 
    36*p^6*s^2 - 36*p^6*s*t + 24*p^6*s - 48*p^6*t + 24*p^6)*x^10*y^2 + (-12*z*p^5*s^2 - 36*z*p^5*s*t + 12*z*p^5*s + 
48*z*p^5*t)*x^10 + 
    (-108*p^6*s^2*t + 72*p^6*s^2 + 12*p^6*s*t - 192*p^6*s + 144*p^6*t + 84*p^6)*x^9*y^3 + (-12*z*p^5*s^2*t - 
48*z*p^5*s^2 - 12*z*p^5*s*t + 
    84*z*p^5*s + 24*z*p^5*t - 24*z*p^5)*x^9*y + (-24*w*A*p^4*s^2*t + 24*w*A*p^4*s^2 + 24*w*A*p^4*s*t - 
24*w*A*p^4*s)*x^9 + (-120*p^6*s^2*t 
    + 120*p^6*s^2 + 12*p^6*s*t - 216*p^6*s + 192*p^6*t)*x^8*y^4 + (180*z*p^5*s^2*t - 198*z*p^5*s^2 - 84*z*p^5*s*t + 
432*z*p^5*s - 
    168*z*p^5*t - 168*z*p^5)*x^8*y^2 + (-12*w*A*p^4*s^2*t + 36*w*A*p^4*s^2 + 24*w*A*p^4*s*t - 48*w*A*p^4*s)*x^8*y + 
(78*z^2*p^4*s^2*t - 
    36*z^2*p^4*s^2 - 60*z^2*p^4*s*t + 36*z^2*p^4*s - 24*z^2*p^4*t - 12*z^2*p^4)*x^8 + (48*p^6*s^2*t + 48*p^6*s^2 - 
24*p^6*s*t - 24*p^6*s - 
    24*p^6*t - 132*p^6)*x^7*y^5 + (24*z*p^5*s^2*t - 102*z*p^5*s^2 - 12*z*p^5*s*t + 468*z*p^5*s - 120*z*p^5*t - 
180*z*p^5)*x^7*y^3 + 
    (72*w*A*p^4*s^2*t - 36*w*A*p^4*s^2 - 36*w*A*p^4*s*t - 48*w*A*p^4*t + 24*w*A*p^4)*x^7*y^2 + (-114*z^2*p^4*s^2*t + 
66*z^2*p^4*s^2 + 
    72*z^2*p^4*s*t - 294*z^2*p^4*s + 72*z^2*p^4*t + 192*z^2*p^4)*x^7*y + (24*z*w*A*p^3*s^2*t - 36*z*w*A*p^3*s^2 - 
36*z*w*A*p^3*s*t + 
    60*z*w*A*p^3*s)*x^7 + (96*p^6*s^2*t - 12*p^6*s^2 + 12*p^6*s*t - 216*p^6*t + 12*p^6)*x^6*y^6 + (-48*z*p^5*s^2*t - 
24*z*p^5*s^2 - 
    36*z*p^5*s*t + 384*z*p^5*s + 144*z*p^5*t - 192*z*p^5)*x^6*y^4 + (36*w*A*p^4*s^2*t - 48*w*A*p^4*s^2 - 
12*w*A*p^4*s*t + 72*w*A*p^4*s - 
    120*w*A*p^4*t + 60*w*A*p^4)*x^6*y^3 + (54*z^2*p^4*s^2*t - 54*z^2*p^4*s^2 + 24*z^2*p^4*s*t - 438*z^2*p^4*s - 
12*z^2*p^4*t + 
    348*z^2*p^4)*x^6*y^2 + (-60*z*w*A*p^3*s^2*t + 24*z*w*A*p^3*s^2 - 36*z*w*A*p^3*s*t + 36*z*w*A*p^3*s + 
144*z*w*A*p^3*t - 
    72*z*w*A*p^3)*x^6*y + (66*z^3*p^3*s^2*t - 48*z^3*p^3*s^2 - 12*z^3*p^3*s*t + 102*z^3*p^3*s - 72*z^3*p^3*t - 
66*z^3*p^3)*x^6 + 
    (84*p^6*s^2*t - 24*p^6*s^2 + 60*p^6*s*t - 72*p^6*s - 264*p^6*t + 216*p^6)*x^5*y^7 + (-252*z*p^5*s^2*t + 
126*z*p^5*s^2 + 72*z*p^5*s*t + 
    72*z*p^5*s + 432*z*p^5*t - 324*z*p^5)*x^5*y^5 + (12*w*A*p^4*s^2*t - 24*w*A*p^4*s^2 + 36*w*A*p^4*s*t + 
144*w*A*p^4*s - 96*w*A*p^4*t + 
    48*w*A*p^4)*x^5*y^4 + (36*z^2*p^4*s^2*t - 72*z^2*p^4*s^2 + 36*z^2*p^4*s*t - 504*z^2*p^4*s - 216*z^2*p^4*t + 
540*z^2*p^4)*x^5*y^3 + 
    (-72*z*w*A*p^3*s^2*t + 18*z*w*A*p^3*s^2 - 36*z*w*A*p^3*s*t - 72*z*w*A*p^3*s + 288*z*w*A*p^3*t - 
144*z*w*A*p^3)*x^5*y^2 + 
    (-96*z^3*p^3*s^2*t + 66*z^3*p^3*s^2 + 60*z^3*p^3*s*t + 144*z^3*p^3*s - 12*z^3*p^3*t - 126*z^3*p^3)*x^5*y + 
(90*z^2*w*A*p^2*s^2*t - 
    72*z^2*w*A*p^2*s^2 - 36*z^2*w*A*p^2*s*t + 36*z^2*w*A*p^2*s - 72*z^2*w*A*p^2*t + 36*z^2*w*A*p^2)*x^5 + 
(72*p^6*s^2*t - 36*p^6*s^2 + 
    36*p^6*s*t - 120*p^6*s - 120*p^6*t + 276*p^6)*x^4*y^8 + (-180*z*p^5*s^2*t + 78*z*p^5*s^2 + 36*z*p^5*s*t + 
84*z*p^5*s + 336*z*p^5*t - 
    468*z*p^5)*x^4*y^6 + (-24*w*A*p^4*s^2*t + 12*w*A*p^4*s^2 + 84*w*A*p^4*s*t + 48*w*A*p^4*s + 24*w*A*p^4*t - 
12*w*A*p^4)*x^4*y^5 + 
    (270*z^2*p^4*s^2*t - 162*z^2*p^4*s^2 - 108*z^2*p^4*s*t - 162*z^2*p^4*s - 396*z^2*p^4*t + 576*z^2*p^4)*x^4*y^4 + 
(-132*z*w*A*p^3*s^2*t +
    30*z*w*A*p^3*s^2 - 12*z*w*A*p^3*s*t - 228*z*w*A*p^3*s + 216*z*w*A*p^3*t - 108*z*w*A*p^3)*x^4*y^3 + 
(18*z^3*p^3*s^2*t + 18*z^3*p^3*s^2 +
    36*z^3*p^3*s*t + 366*z^3*p^3*s + 60*z^3*p^3*t - 408*z^3*p^3)*x^4*y^2 + (54*z^2*w*A*p^2*s^2*t - 18*z^2*w*A*p^2*s^2 
- 36*z^2*w*A*p^2*s*t 
    + 54*z^2*w*A*p^2*s - 144*z^2*w*A*p^2*t + 72*z^2*w*A*p^2)*x^4*y + (18*z^4*p^2*s^2*t - 48*z^4*p^2*s^2 - 60*z^4*p^2*s 
+ 12*z^4*p^2*t + 
    72*z^4*p^2)*x^4 + (-108*p^6*s^2*t + 72*p^6*s^2 + 48*p^6*s*t - 264*p^6*s + 144*p^6*t + 192*p^6)*x^3*y^9 + 
(60*z*p^5*s^2*t - 84*z*p^5*s^2
    - 48*z*p^5*s*t + 372*z*p^5*s - 120*z*p^5*t - 492*z*p^5)*x^3*y^7 + (-84*w*A*p^4*s^2*t + 72*w*A*p^4*s^2 + 
48*w*A*p^4*s*t - 72*w*A*p^4*s +
    96*w*A*p^4*t - 48*w*A*p^4)*x^3*y^6 + (18*z^2*p^4*s^2*t + 18*z^2*p^4*s^2 - 270*z^2*p^4*s - 72*z^2*p^4*t + 
576*z^2*p^4)*x^3*y^5 + 
    (60*z*w*A*p^3*s^2*t - 84*z*w*A*p^3*s^2 - 96*z*w*A*p^3*s*t - 60*z*w*A*p^3*s - 72*z*w*A*p^3*t + 
36*z*w*A*p^3)*x^3*y^4 + 
    (-198*z^3*p^3*s^2*t + 90*z^3*p^3*s^2 + 120*z^3*p^3*s*t + 186*z^3*p^3*s + 180*z^3*p^3*t - 474*z^3*p^3)*x^3*y^3 + 
(126*z^2*w*A*p^2*s^2*t 
    - 18*z^2*w*A*p^2*s^2 - 36*z^2*w*A*p^2*s*t + 162*z^2*w*A*p^2*s - 144*z^2*w*A*p^2*t + 72*z^2*w*A*p^2)*x^3*y^2 + 
(-12*z^4*p^2*s^2*t - 
    48*z^4*p^2*s^2 - 48*z^4*p^2*s*t - 114*z^4*p^2*s - 12*z^4*p^2*t + 210*z^4*p^2)*x^3*y + (-6*z^3*w*A*p*s^2*t - 
12*z^3*w*A*p*s^2 + 
    24*z^3*w*A*p*s*t - 42*z^3*w*A*p*s + 36*z^3*w*A*p*t - 18*z^3*w*A*p)*x^3 + (-120*p^6*s^2*t + 120*p^6*s^2 - 
24*p^6*s*t - 144*p^6*s + 
    192*p^6*t)*x^2*y^10 + (216*z*p^5*s^2*t - 216*z*p^5*s^2 - 48*z*p^5*s*t + 468*z*p^5*s - 384*z*p^5*t - 
168*z*p^5)*x^2*y^8 + 
    (-72*w*A*p^4*s^2*t + 120*w*A*p^4*s^2 - 24*w*A*p^4*s*t - 96*w*A*p^4*s + 96*w*A*p^4*t - 48*w*A*p^4)*x^2*y^7 + 
(-102*z^2*p^4*s^2*t + 
    108*z^2*p^4*s^2 + 12*z^2*p^4*s*t - 324*z^2*p^4*s + 264*z^2*p^4*t + 384*z^2*p^4)*x^2*y^6 + (96*z*w*A*p^3*s^2*t - 
120*z*w*A*p^3*s^2 + 
    48*z*w*A*p^3*s*t + 12*z*w*A*p^3*s - 240*z*w*A*p^3*t + 120*z*w*A*p^3)*x^2*y^5 + (6*z^3*p^3*s^2*t - 24*z^3*p^3*s^2 - 
24*z^3*p^3*s*t + 
    288*z^3*p^3*s + 12*z^3*p^3*t - 450*z^3*p^3)*x^2*y^4 + (-78*z^2*w*A*p^2*s^2*t + 120*z^2*w*A*p^2*s^2 + 
60*z^2*w*A*p^2*s*t - 
    12*z^2*w*A*p^2*s + 96*z^2*w*A*p^2*t - 48*z^2*w*A*p^2)*x^2*y^3 + (72*z^4*p^2*s^2*t - 90*z^4*p^2*s^2 - 
48*z^4*p^2*s*t - 72*z^4*p^2*s - 
    60*z^4*p^2*t + 318*z^4*p^2)*x^2*y^2 + (-66*z^3*w*A*p*s^2*t + 24*z^3*w*A*p*s*t - 84*z^3*w*A*p*s + 36*z^3*w*A*p*t - 
18*z^3*w*A*p)*x^2*y +
    (6*z^5*p*s^2*t + 12*z^5*p*s*t + 12*z^5*p*t - 30*z^5*p)*x^2 + (48*p^6*s^2*t + 48*p^6*s^2 - 96*p^6*s*t + 120*p^6*s + 
48*p^6*t - 
    168*p^6)*x*y^11 + (60*z*p^5*s^2*t - 120*z*p^5*s^2 + 96*z*p^5*s*t + 144*z*p^5*s - 192*z*p^5*t + 72*z*p^5)*x*y^9 + 
(48*w*A*p^4*s^2*t + 
    48*w*A*p^4*s^2 - 96*w*A*p^4*s*t - 24*w*A*p^4*s + 48*w*A*p^4*t - 24*w*A*p^4)*x*y^8 + (-168*z^2*p^4*s^2*t + 
120*z^2*p^4*s^2 + 
    36*z^2*p^4*s*t - 276*z^2*p^4*s + 288*z^2*p^4*t + 84*z^2*p^4)*x*y^7 + (84*z*w*A*p^3*s^2*t - 120*z*w*A*p^3*s^2 + 
96*z*w*A*p^3*s*t + 
    48*z*w*A*p^3*s - 240*z*w*A*p^3*t + 120*z*w*A*p^3)*x*y^6 + (102*z^3*p^3*s^2*t - 60*z^3*p^3*s^2 - 96*z^3*p^3*s*t + 
282*z^3*p^3*s - 
    96*z^3*p^3*t - 312*z^3*p^3)*x*y^5 + (-96*z^2*w*A*p^2*s^2*t + 120*z^2*w*A*p^2*s^2 + 12*z^2*w*A*p^2*s*t - 
60*z^2*w*A*p^2*s + 
    168*z^2*w*A*p^2*t - 84*z^2*w*A*p^2)*x*y^4 + (-48*z^4*p^2*s^2*t - 12*z^4*p^2*s^2 + 60*z^4*p^2*s*t - 162*z^4*p^2*s + 
24*z^4*p^2*t + 
    288*z^4*p^2)*x*y^3 + (54*z^3*w*A*p*s^2*t - 72*z^3*w*A*p*s^2 - 36*z^3*w*A*p*s*t + 18*z^3*w*A*p*s - 48*z^3*w*A*p*t + 
24*z^3*w*A*p)*x*y^2 
    + (-6*z^5*p*s^2*t + 12*z^5*p*s^2 + 12*z^5*p*s - 96*z^5*p)*x*y + (6*z^4*w*A*s^2*t + 6*z^4*w*A*s)*x + (48*p^6*s^2*t 
- 48*p^6*s*t + 
    96*p^6*s - 96*p^6)*y^12 + (-48*z*p^5*s^2*t - 24*z*p^5*s^2 + 72*z*p^5*s*t - 48*z*p^5*s + 96*z*p^5)*y^10 + 
(48*w*A*p^4*s^2*t - 
    48*w*A*p^4*s*t)*y^9 + (-72*z^2*p^4*s^2*t + 36*z^2*p^4*s^2 - 12*z^2*p^4*s*t - 96*z^2*p^4*s + 96*z^2*p^4*t + 
24*z^2*p^4)*y^8 + 
    (-24*z*w*A*p^3*s^2 + 72*z*w*A*p^3*s*t - 96*z*w*A*p^3*t + 48*z*w*A*p^3)*y^7 + (102*z^3*p^3*s^2*t - 54*z^3*p^3*s^2 - 
84*z^3*p^3*s*t + 
    168*z^3*p^3*s - 72*z^3*p^3*t - 108*z^3*p^3)*y^6 + (-96*z^2*w*A*p^2*s^2*t + 84*z^2*w*A*p^2*s^2 + 36*z^2*w*A*p^2*s*t 
- 72*z^2*w*A*p^2*s +
    96*z^2*w*A*p^2*t - 48*z^2*w*A*p^2)*y^5 + (-30*z^4*p^2*s^2*t - 6*z^4*p^2*s^2 + 36*z^4*p^2*s*t - 84*z^4*p^2*s + 
36*z^4*p^2*t + 
    132*z^4*p^2)*y^4 + (18*z^3*w*A*p*s^2*t - 30*z^3*w*A*p*s^2 - 12*z^3*w*A*p*s*t - 24*z^3*w*A*p*t + 12*z^3*w*A*p)*y^3 
+ (-12*z^5*p*s*t + 
    12*z^5*p*s - 12*z^5*p*t - 30*z^5*p)*y^2 + (-6*z^4*w*A*s^2*t + 6*z^4*w*A*s^2)*y + 6*z^6;

denominator :=
(-16*p^6*s^2*t + 32*p^6*s^2 - 16*p^6*s*t - 16*p^6*s + 48*p^6*t - 32*p^6)*x^12 + (-8*p^6*s^2*t + 16*p^6*s^2 + 16*p^6*t 
- 32*p^6)*x^11*y + 
    (24*p^6*s^2*t - 24*p^6*s^2 + 16*p^6*s*t + 40*p^6*s - 64*p^6*t - 16*p^6)*x^10*y^2 + (16*z*p^5*s^2*t + 16*z*p^5*s^2 
- 56*z*p^5*s*t + 
    16*z*p^5*s + 48*z*p^5*t - 48*z*p^5)*x^10 + (-48*p^6*s^2*t + 72*p^6*s^2 + 16*p^6*s*t - 56*p^6*s + 56*p^6*t - 
40*p^6)*x^9*y^3 + 
    (8*z*p^5*s^2*t - 16*z*p^5*s^2 - 48*z*p^5*s*t + 24*z*p^5*s + 56*z*p^5*t - 16*z*p^5)*x^9*y + (-44*p^6*s^2*t + 
112*p^6*s^2 - 72*p^6*s + 
    88*p^6*t - 104*p^6)*x^8*y^4 + (72*z*p^5*s^2*t - 156*z*p^5*s^2 - 16*z*p^5*s*t + 176*z*p^5*s - 104*z*p^5*t + 
16*z*p^5)*x^8*y^2 + 
    (56*z^2*p^4*s^2*t - 28*z^2*p^4*s^2 - 16*z^2*p^4*s*t + 56*z^2*p^4*s - 72*z^2*p^4*t - 24*z^2*p^4)*x^8 + 
(24*p^6*s^2*t + 48*p^6*s^2 - 
    12*p^6*s - 40*p^6*t - 88*p^6)*x^7*y^5 + (-8*z*p^5*s^2*t - 152*z*p^5*s^2 + 24*z*p^5*s*t + 156*z*p^5*s - 48*z*p^5*t 
+ 72*z*p^5)*x^7*y^3 +
    (-20*z^2*p^4*s^2*t + 76*z^2*p^4*s^2 - 108*z^2*p^4*s + 40*z^2*p^4*t + 16*z^2*p^4)*x^7*y + (28*p^6*s^2*t + 
28*p^6*s^2 + 16*p^6*s*t - 
    8*p^6*s - 80*p^6*t - 36*p^6)*x^6*y^6 + (-4*z*p^5*s^2*t - 136*z*p^5*s^2 + 192*z*p^5*s + 24*z*p^5*t)*x^6*y^4 + 
(36*z^2*p^4*s^2*t + 
    72*z^2*p^4*s^2 - 20*z^2*p^4*s*t - 140*z^2*p^4*s - 4*z^2*p^4*t + 32*z^2*p^4)*x^6*y^2 + (20*z^3*p^3*s^2*t - 
40*z^3*p^3*s^2 + 
    20*z^3*p^3*s*t + 56*z^3*p^3*s - 60*z^3*p^3*t - 4*z^3*p^3)*x^6 + (16*p^6*s^2*t + 4*p^6*s^2 + 48*p^6*s*t - 24*p^6*s 
- 80*p^6*t + 
    16*p^6)*x^5*y^7 + (-96*z*p^5*s^2*t + 12*z*p^5*s^2 - 48*z*p^5*s*t + 60*z*p^5*s + 240*z*p^5*t - 120*z*p^5)*x^5*y^5 + 
(-12*z^2*p^4*s^2*t +
    78*z^2*p^4*s^2 + 12*z^2*p^4*s*t - 240*z^2*p^4*s - 24*z^2*p^4*t + 156*z^2*p^4)*x^5*y^3 + (-44*z^3*p^3*s^2*t - 
2*z^3*p^3*s^2 + 
    36*z^3*p^3*s*t + 36*z^3*p^3*s + 16*z^3*p^3*t - 32*z^3*p^3)*x^5*y + (24*p^6*s^2*t - 24*p^6*s^2 + 40*p^6*s*t - 
8*p^6*s - 88*p^6*t + 
    32*p^6)*x^4*y^8 + (-68*z*p^5*s^2*t + 76*z*p^5*s^2 - 104*z*p^5*s*t + 76*z*p^5*s + 216*z*p^5*t - 204*z*p^5)*x^4*y^6 
+ (96*z^2*p^4*s^2*t -
    84*z^2*p^4*s^2 + 72*z^2*p^4*s*t - 36*z^2*p^4*s - 276*z^2*p^4*t + 228*z^2*p^4)*x^4*y^4 + (24*z^3*p^3*s^2*t - 
60*z^3*p^3*s^2 + 
    16*z^3*p^3*s*t + 166*z^3*p^3*s - 28*z^3*p^3*t - 124*z^3*p^3)*x^4*y^2 + (-2*z^4*p^2*s^2*t - 2*z^4*p^2*s^2 - 
2*z^4*p^2*s*t - 20*z^4*p^2*s
    + 12*z^4*p^2*t + 24*z^4*p^2)*x^4 + (-48*p^6*s^2*t + 72*p^6*s^2 - 8*p^6*s*t - 80*p^6*s + 80*p^6*t - 16*p^6)*x^3*y^9 
+ (32*z*p^5*s^2*t - 
    28*z*p^5*s^2 - 24*z*p^5*s*t + 120*z*p^5*s - 16*z*p^5*t - 88*z*p^5)*x^3*y^7 + (-12*z^2*p^4*s^2*t - 12*z^2*p^4*s^2 + 
84*z^2*p^4*s*t - 
    132*z^2*p^4*s - 48*z^2*p^4*t + 168*z^2*p^4)*x^3*y^5 + (-84*z^3*p^3*s^2*t + 54*z^3*p^3*s^2 - 20*z^3*p^3*s*t + 
16*z^3*p^3*s + 
    176*z^3*p^3*t - 172*z^3*p^3)*x^3*y^3 + (-10*z^4*p^2*s^2*t + 2*z^4*p^2*s^2 - 12*z^4*p^2*s*t - 48*z^4*p^2*s + 
20*z^4*p^2*t + 
    56*z^4*p^2)*x^3*y + (-32*p^6*s^2*t + 88*p^6*s^2 - 48*p^6*s*t - 48*p^6*s + 112*p^6*t - 80*p^6)*x^2*y^10 + 
(96*z*p^5*s^2*t - 
    168*z*p^5*s^2 + 104*z*p^5*s*t + 152*z*p^5*s - 272*z*p^5*t + 64*z*p^5)*x^2*y^8 + (-40*z^2*p^4*s^2*t + 
56*z^2*p^4*s^2 - 28*z^2*p^4*s*t - 
    64*z^2*p^4*s + 96*z^2*p^4*t - 36*z^2*p^4)*x^2*y^6 + (-8*z^3*p^3*s^2*t - 14*z^3*p^3*s^2 - 36*z^3*p^3*s*t + 
108*z^3*p^3*s + 16*z^3*p^3*t 
    - 104*z^3*p^3)*x^2*y^4 + (36*z^4*p^2*s^2*t - 48*z^4*p^2*s^2 + 2*z^4*p^2*s*t + 14*z^4*p^2*s - 68*z^4*p^2*t + 
88*z^4*p^2)*x^2*y^2 + 
    (2*z^5*p*s^2*t - z^5*p*s^2 + 2*z^5*p*s*t + 2*z^5*p*s - 6*z^5*p)*x^2 + (48*p^6*s^2*t - 48*p^6*s*t + 48*p^6*s - 
16*p^6*t - 64*p^6)*x*y^11
    + (-8*z*p^5*s^2*t - 80*z*p^5*s^2 + 120*z*p^5*s*t - 144*z*p^5*t + 120*z*p^5)*x*y^9 + (-80*z^2*p^4*s^2*t + 
124*z^2*p^4*s^2 - 
    72*z^2*p^4*s*t - 72*z^2*p^4*s + 208*z^2*p^4*t - 104*z^2*p^4)*x*y^7 + (36*z^3*p^3*s^2*t - 48*z^3*p^3*s^2 - 
12*z^3*p^3*s*t + 84*z^3*p^3*s
    - 40*z^3*p^3*t - 16*z^3*p^3)*x*y^5 + (-14*z^4*p^2*s^2*t + 4*z^4*p^2*s^2 + 12*z^4*p^2*s*t - 54*z^4*p^2*s + 
24*z^4*p^2*t + 
    54*z^4*p^2)*x*y^3 + (-2*z^5*p*s^2*t + 4*z^5*p*s^2 + 4*z^5*p*t - 20*z^5*p)*x*y + (32*p^6*s^2*t - 16*p^6*s^2 - 
16*p^6*s*t + 32*p^6*s - 
    32*p^6*t - 16*p^6)*y^12 + (-40*z*p^5*s^2*t + 8*z*p^5*s^2 + 48*z*p^5*s*t - 48*z*p^5*s + 48*z*p^5)*y^10 + 
(-24*z^2*p^4*s^2*t + 
    48*z^2*p^4*s^2 - 32*z^2*p^4*s*t - 8*z^2*p^4*s + 80*z^2*p^4*t - 64*z^2*p^4)*y^8 + (56*z^3*p^3*s^2*t - 
64*z^3*p^3*s^2 - 4*z^3*p^3*s*t + 
    80*z^3*p^3*s - 80*z^3*p^3*t)*y^6 + (-10*z^4*p^2*s^2*t + 2*z^4*p^2*s^2 - 24*z^4*p^2*s + 12*z^4*p^2*t + 
24*z^4*p^2)*y^4 + (-2*z^5*p*s*t +
    4*z^5*p*s - 4*z^5*p*t - 4*z^5*p)*y^2 + z^6;

// print numerator, denominator;
/* The numerator (the one after we've divided by the gcd) above can be expressed as a sum
    
p1 + s*p2 + s^2*p3 + t*p4 + t*s*p5 + t*s^2*p6
    
where the p_i are in Q(sqrt{p})[x,y,z,w]. We find the p_i below. */
    
P4<t4> := PolynomialRing(N2);
f4 := t4^2 - p;
N4<A> := NumberField(f4);

NewPolys<a1,b1,a2,b2,a3,b3,a4,b4,a5,b5,a6,b6> := PolynomialRing(N4,12);
    
T<w,z,y,x> := PolynomialRing(NewPolys,4);

p1 :=
12*p*A*w*z^3*y^3 +
    24*p*A*w*z^3*y^2*x - 18*p*A*w*z^3*y*x^2 - 18*p*A*w*z^3*x^3 - 48*p^2*A*w*z^2*y^5 -
    84*p^2*A*w*z^2*y^4*x - 48*p^2*A*w*z^2*y^3*x^2 + 72*p^2*A*w*z^2*y^2*x^3 + 72*p^2*A*w*z^2*y*x^4 +
    36*p^2*A*w*z^2*x^5 + 48*p^3*A*w*z*y^7 + 120*p^3*A*w*z*y^6*x + 120*p^3*A*w*z*y^5*x^2 +
    36*p^3*A*w*z*y^4*x^3 - 108*p^3*A*w*z*y^3*x^4 - 144*p^3*A*w*z*y^2*x^5 - 72*p^3*A*w*z*y*x^6 -
    24*p^4*A*w*y^8*x - 48*p^4*A*w*y^7*x^2 - 48*p^4*A*w*y^6*x^3 - 12*p^4*A*w*y^5*x^4 +
    48*p^4*A*w*y^4*x^5 + 60*p^4*A*w*y^3*x^6 + 24*p^4*A*w*y^2*x^7 + 6*z^6 - 30*p*z^5*y^2 -
    96*p*z^5*y*x - 30*p*z^5*x^2 + 132*p^2*z^4*y^4 + 288*p^2*z^4*y^3*x + 318*p^2*z^4*y^2*x^2 +
    210*p^2*z^4*y*x^3 + 72*p^2*z^4*x^4 - 108*p^3*z^3*y^6 - 312*p^3*z^3*y^5*x - 450*p^3*z^3*y^4*x^2 -
    474*p^3*z^3*y^3*x^3 - 408*p^3*z^3*y^2*x^4 - 126*p^3*z^3*y*x^5 - 66*p^3*z^3*x^6 + 24*p^4*z^2*y^8 +
    84*p^4*z^2*y^7*x + 384*p^4*z^2*y^6*x^2 + 576*p^4*z^2*y^5*x^3 + 576*p^4*z^2*y^4*x^4 +
    540*p^4*z^2*y^3*x^5 + 348*p^4*z^2*y^2*x^6 + 192*p^4*z^2*y*x^7 - 12*p^4*z^2*x^8 + 96*p^5*z*y^10
    + 72*p^5*z*y^9*x - 168*p^5*z*y^8*x^2 - 492*p^5*z*y^7*x^3 - 468*p^5*z*y^6*x^4 -
    324*p^5*z*y^5*x^5 - 192*p^5*z*y^4*x^6 - 180*p^5*z*y^3*x^7 - 168*p^5*z*y^2*x^8 - 24*p^5*z*y*x^9
    - 96*p^6*y^12 - 168*p^6*y^11*x + 192*p^6*y^9*x^3 + 276*p^6*y^8*x^4 + 216*p^6*y^7*x^5 +
    12*p^6*y^6*x^6 - 132*p^6*y^5*x^7 + 84*p^6*y^3*x^9 + 24*p^6*y^2*x^10;

p2 :=
6*A*w*z^4*x + 18*p*A*w*z^3*y^2*x - 84*p*A*w*z^3*y*x^2 - 42*p*A*w*z^3*x^3 - 72*p^2*A*w*z^2*y^5 -
    60*p^2*A*w*z^2*y^4*x - 12*p^2*A*w*z^2*y^3*x^2 + 162*p^2*A*w*z^2*y^2*x^3 + 54*p^2*A*w*z^2*y*x^4 +
    36*p^2*A*w*z^2*x^5 + 48*p^3*A*w*z*y^6*x + 12*p^3*A*w*z*y^5*x^2 - 60*p^3*A*w*z*y^4*x^3 -
    228*p^3*A*w*z*y^3*x^4 - 72*p^3*A*w*z*y^2*x^5 + 36*p^3*A*w*z*y*x^6 + 60*p^3*A*w*z*x^7 -  
    24*p^4*A*w*y^8*x - 96*p^4*A*w*y^7*x^2 - 72*p^4*A*w*y^6*x^3 + 48*p^4*A*w*y^5*x^4 +
    144*p^4*A*w*y^4*x^5 + 72*p^4*A*w*y^3*x^6 - 48*p^4*A*w*y*x^8 - 24*p^4*A*w*x^9 + 12*p*z^5*y^2 +
    12*p*z^5*y*x - 84*p^2*z^4*y^4 - 162*p^2*z^4*y^3*x - 72*p^2*z^4*y^2*x^2 - 114*p^2*z^4*y*x^3 -
    60*p^2*z^4*x^4 + 168*p^3*z^3*y^6 + 282*p^3*z^3*y^5*x + 288*p^3*z^3*y^4*x^2 + 186*p^3*z^3*y^3*x^3 
+
    366*p^3*z^3*y^2*x^4 + 144*p^3*z^3*y*x^5 + 102*p^3*z^3*x^6 - 96*p^4*z^2*y^8 - 276*p^4*z^2*y^7*x -
    324*p^4*z^2*y^6*x^2 - 270*p^4*z^2*y^5*x^3 - 162*p^4*z^2*y^4*x^4 - 504*p^4*z^2*y^3*x^5 -
    438*p^4*z^2*y^2*x^6 - 294*p^4*z^2*y*x^7 + 36*p^4*z^2*x^8 - 48*p^5*z*y^10 + 144*p^5*z*y^9*x +
    468*p^5*z*y^8*x^2 + 372*p^5*z*y^7*x^3 + 84*p^5*z*y^6*x^4 + 72*p^5*z*y^5*x^5 +
    384*p^5*z*y^4*x^6 + 468*p^5*z*y^3*x^7 + 432*p^5*z*y^2*x^8 + 84*p^5*z*y*x^9 + 12*p^5*z*x^10
    + 96*p^6*y^12 + 120*p^6*y^11*x - 144*p^6*y^10*x^2 - 264*p^6*y^9*x^3 - 120*p^6*y^8*x^4 -
    72*p^6*y^7*x^5 - 24*p^6*y^5*x^7 - 216*p^6*y^4*x^8 - 192*p^6*y^3*x^9 + 24*p^6*y^2*x^10 -
    24*p^6*x^12;

p3 := 
24*p^6*x^12 + 12*p^6*x^11*y - 36*p^6*x^10*y^2 - 12*z*p^5*x^10 + 72*p^6*x^9*y^3 - 48*z*p^5*x^9*y + 
24*w*A*p^4*x^9 + 120*p^6*x^8*y^4 - 198*z*p^5*x^8*y^2 + 36*w*A*p^4*x^8*y - 36*z^2*p^4*x^8 + 48*p^6*x^7*y^5 - 
102*z*p^5*x^7*y^3 - 36*w*A*p^4*x^7*y^2 + 66*z^2*p^4*x^7*y - 36*z*w*A*p^3*x^7 - 12*p^6*x^6*y^6 - 24*z*p^5*x^6*y^4 - 
48*w*A*p^4*x^6*y^3 -
    54*z^2*p^4*x^6*y^2 + 24*z*w*A*p^3*x^6*y - 48*z^3*p^3*x^6 - 24*p^6*x^5*y^7 + 126*z*p^5*x^5*y^5 - 
24*w*A*p^4*x^5*y^4 
    - 72*z^2*p^4*x^5*y^3 + 18*z*w*A*p^3*x^5*y^2 + 66*z^3*p^3*x^5*y - 72*z^2*w*A*p^2*x^5 - 
36*p^6*x^4*y^8 + 
    78*z*p^5*x^4*y^6 + 12*w*A*p^4*x^4*y^5 - 162*z^2*p^4*x^4*y^4 + 30*z*w*A*p^3*x^4*y^3 + 
18*z^3*p^3*x^4*y^2 - 
    18*z^2*w*A*p^2*x^4*y - 48*z^4*p^2*x^4 + 72*p^6*x^3*y^9 - 84*z*p^5*x^3*y^7 + 72*w*A*p^4*x^3*y^6 + 
18*z^2*p^4*x^3*y^5
    - 84*z*w*A*p^3*x^3*y^4 + 90*z^3*p^3*x^3*y^3 - 18*z^2*w*A*p^2*x^3*y^2 - 48*z^4*p^2*x^3*y - 
12*z^3*w*A*p*x^3 + 
    120*p^6*x^2*y^10 - 216*z*p^5*x^2*y^8 + 120*w*A*p^4*x^2*y^7 + 108*z^2*p^4*x^2*y^6 - 
120*z*w*A*p^3*x^2*y^5 - 
    24*z^3*p^3*x^2*y^4 + 120*z^2*w*A*p^2*x^2*y^3 - 90*z^4*p^2*x^2*y^2 + 48*p^6*x*y^11 - 
120*z*p^5*x*y^9 + 
    48*w*A*p^4*x*y^8 + 120*z^2*p^4*x*y^7 - 120*z*w*A*p^3*x*y^6 - 60*z^3*p^3*x*y^5 + 
120*z^2*w*A*p^2*x*y^4 - 
    12*z^4*p^2*x*y^3 - 72*z^3*w*A*p*x*y^2 + 12*z^5*p*x*y - 24*z*p^5*y^10 + 36*z^2*p^4*y^8 - 
24*z*w*A*p^3*y^7 - 
    54*z^3*p^3*y^6 + 84*z^2*w*A*p^2*y^5 - 6*z^4*p^2*y^4 - 30*z^3*w*A*p*y^3 + 6*z^4*w*A*y;

p4 :=
-24*p*A*w*z^3*y^3 - 48*p*A*w*z^3*y^2*x + 36*p*A*w*z^3*y*x^2 + 36*p*A*w*z^3*x^3 + 96*p^2*A*w*z^2*y^5
    + 168*p^2*A*w*z^2*y^4*x + 96*p^2*A*w*z^2*y^3*x^2 - 144*p^2*A*w*z^2*y^2*x^3 - 144*p^2*A*w*z^2*y*x^4 -
    72*p^2*A*w*z^2*x^5 - 96*p^3*A*w*z*y^7 - 240*p^3*A*w*z*y^6*x - 240*p^3*A*w*z*y^5*x^2 -
    72*p^3*A*w*z*y^4*x^3 + 216*p^3*A*w*z*y^3*x^4 + 288*p^3*A*w*z*y^2*x^5 + 144*p^3*A*w*z*y*x^6 +
    48*p^4*A*w*y^8*x + 96*p^4*A*w*y^7*x^2 + 96*p^4*A*w*y^6*x^3 + 24*p^4*A*w*y^5*x^4 -
    96*p^4*A*w*y^4*x^5 - 120*p^4*A*w*y^3*x^6 - 48*p^4*A*w*y^2*x^7 - 12*p*z^5*y^2 + 12*p*z^5*x^2 +
    36*p^2*z^4*y^4 + 24*p^2*z^4*y^3*x - 60*p^2*z^4*y^2*x^2 - 12*p^2*z^4*y*x^3 + 12*p^2*z^4*x^4 -
    72*p^3*z^3*y^6 - 96*p^3*z^3*y^5*x + 12*p^3*z^3*y^4*x^2 + 180*p^3*z^3*y^3*x^3 + 60*p^3*z^3*y^2*x^4 -
    12*p^3*z^3*y*x^5 - 72*p^3*z^3*x^6 + 96*p^4*z^2*y^8 + 288*p^4*z^2*y^7*x + 264*p^4*z^2*y^6*x^2 -
    72*p^4*z^2*y^5*x^3 - 396*p^4*z^2*y^4*x^4 - 216*p^4*z^2*y^3*x^5 - 12*p^4*z^2*y^2*x^6 +
    72*p^4*z^2*y*x^7 - 24*p^4*z^2*x^8 - 192*p^5*z*y^9*x - 384*p^5*z*y^8*x^2 - 120*p^5*z*y^7*x^3 +
    336*p^5*z*y^6*x^4 + 432*p^5*z*y^5*x^5 + 144*p^5*z*y^4*x^6 - 120*p^5*z*y^3*x^7 -
    168*p^5*z*y^2*x^8 + 24*p^5*z*y*x^9 + 48*p^5*z*x^10 + 48*p^6*y^11*x + 192*p^6*y^10*x^2 +
    144*p^6*y^9*x^3 - 120*p^6*y^8*x^4 - 264*p^6*y^7*x^5 - 216*p^6*y^6*x^6 - 24*p^6*y^5*x^7 + 
    192*p^6*y^4*x^8 + 144*p^6*y^3*x^9 - 48*p^6*y^2*x^10 - 48*p^6*y*x^11;

p5 :=
24*p^6*x^12 + 24*p^6*x^11*y - 36*p^6*x^10*y^2 - 36*z*p^5*x^10 + 12*p^6*x^9*y^3 - 12*z*p^5*x^9*y + 
24*w*A*p^4*x^9 + 
    12*p^6*x^8*y^4 - 84*z*p^5*x^8*y^2 + 24*w*A*p^4*x^8*y - 60*z^2*p^4*x^8 - 24*p^6*x^7*y^5 - 
12*z*p^5*x^7*y^3 - 
    36*w*A*p^4*x^7*y^2 + 72*z^2*p^4*x^7*y - 36*z*w*A*p^3*x^7 + 12*p^6*x^6*y^6 - 36*z*p^5*x^6*y^4 - 
12*w*A*p^4*x^6*y^3 +
    24*z^2*p^4*x^6*y^2 - 36*z*w*A*p^3*x^6*y - 12*z^3*p^3*x^6 + 60*p^6*x^5*y^7 + 72*z*p^5*x^5*y^5 + 
36*w*A*p^4*x^5*y^4 +
    36*z^2*p^4*x^5*y^3 - 36*z*w*A*p^3*x^5*y^2 + 60*z^3*p^3*x^5*y - 36*z^2*w*A*p^2*x^5 + 
36*p^6*x^4*y^8 + 
    36*z*p^5*x^4*y^6 + 84*w*A*p^4*x^4*y^5 - 108*z^2*p^4*x^4*y^4 - 12*z*w*A*p^3*x^4*y^3 + 
36*z^3*p^3*x^4*y^2 - 
    36*z^2*w*A*p^2*x^4*y + 48*p^6*x^3*y^9 - 48*z*p^5*x^3*y^7 + 48*w*A*p^4*x^3*y^6 - 
96*z*w*A*p^3*x^3*y^4 + 
    120*z^3*p^3*x^3*y^3 - 36*z^2*w*A*p^2*x^3*y^2 - 48*z^4*p^2*x^3*y + 24*z^3*w*A*p*x^3 - 
24*p^6*x^2*y^10 - 
    48*z*p^5*x^2*y^8 - 24*w*A*p^4*x^2*y^7 + 12*z^2*p^4*x^2*y^6 + 48*z*w*A*p^3*x^2*y^5 - 
24*z^3*p^3*x^2*y^4 + 
    60*z^2*w*A*p^2*x^2*y^3 - 48*z^4*p^2*x^2*y^2 + 24*z^3*w*A*p*x^2*y + 12*z^5*p*x^2 - 96*p^6*x*y^11 
+ 96*z*p^5*x*y^9 - 
    96*w*A*p^4*x*y^8 + 36*z^2*p^4*x*y^7 + 96*z*w*A*p^3*x*y^6 - 96*z^3*p^3*x*y^5 + 
12*z^2*w*A*p^2*x*y^4 + 
    60*z^4*p^2*x*y^3 - 36*z^3*w*A*p*x*y^2 - 48*p^6*y^12 + 72*z*p^5*y^10 - 48*w*A*p^4*y^9 - 
12*z^2*p^4*y^8 + 
    72*z*w*A*p^3*y^7 - 84*z^3*p^3*y^6 + 36*z^2*w*A*p^2*y^5 + 36*z^4*p^2*y^4 - 12*z^3*w*A*p*y^3 - 
12*z^5*p*y^2;

p6 :=
-24*p^6*x^12 + 12*p^6*x^11*y + 72*p^6*x^10*y^2 - 108*p^6*x^9*y^3 - 12*z*p^5*x^9*y - 24*w*A*p^4*x^9 - 
120*p^6*x^8*y^4 + 
    180*z*p^5*x^8*y^2 - 12*w*A*p^4*x^8*y + 78*z^2*p^4*x^8 + 48*p^6*x^7*y^5 + 24*z*p^5*x^7*y^3 + 
72*w*A*p^4*x^7*y^2 - 
    114*z^2*p^4*x^7*y + 24*z*w*A*p^3*x^7 + 96*p^6*x^6*y^6 - 48*z*p^5*x^6*y^4 + 36*w*A*p^4*x^6*y^3 + 
54*z^2*p^4*x^6*y^2 
    - 60*z*w*A*p^3*x^6*y + 66*z^3*p^3*x^6 + 84*p^6*x^5*y^7 - 252*z*p^5*x^5*y^5 + 12*w*A*p^4*x^5*y^4 
+ 
    36*z^2*p^4*x^5*y^3 - 72*z*w*A*p^3*x^5*y^2 - 96*z^3*p^3*x^5*y + 90*z^2*w*A*p^2*x^5 + 
72*p^6*x^4*y^8 - 
    180*z*p^5*x^4*y^6 - 24*w*A*p^4*x^4*y^5 + 270*z^2*p^4*x^4*y^4 - 132*z*w*A*p^3*x^4*y^3 + 
18*z^3*p^3*x^4*y^2 + 
    54*z^2*w*A*p^2*x^4*y + 18*z^4*p^2*x^4 - 108*p^6*x^3*y^9 + 60*z*p^5*x^3*y^7 - 84*w*A*p^4*x^3*y^6 
+ 
    18*z^2*p^4*x^3*y^5 + 60*z*w*A*p^3*x^3*y^4 - 198*z^3*p^3*x^3*y^3 + 126*z^2*w*A*p^2*x^3*y^2 - 
12*z^4*p^2*x^3*y - 
    6*z^3*w*A*p*x^3 - 120*p^6*x^2*y^10 + 216*z*p^5*x^2*y^8 - 72*w*A*p^4*x^2*y^7 - 
102*z^2*p^4*x^2*y^6 + 
    96*z*w*A*p^3*x^2*y^5 + 6*z^3*p^3*x^2*y^4 - 78*z^2*w*A*p^2*x^2*y^3 + 72*z^4*p^2*x^2*y^2 - 
66*z^3*w*A*p*x^2*y + 
    6*z^5*p*x^2 + 48*p^6*x*y^11 + 60*z*p^5*x*y^9 + 48*w*A*p^4*x*y^8 - 168*z^2*p^4*x*y^7 + 
84*z*w*A*p^3*x*y^6 + 
    102*z^3*p^3*x*y^5 - 96*z^2*w*A*p^2*x*y^4 - 48*z^4*p^2*x*y^3 + 54*z^3*w*A*p*x*y^2 - 6*z^5*p*x*y + 
6*z^4*w*A*x + 
    48*p^6*y^12 - 48*z*p^5*y^10 + 48*w*A*p^4*y^9 - 72*z^2*p^4*y^8 + 102*z^3*p^3*y^6 - 
96*z^2*w*A*p^2*y^5 - 
    30*z^4*p^2*y^4 + 18*z^3*w*A*p*y^3 - 6*z^4*w*A*y;

// print p1 + s*p2 + s^2*p3 + t*p4 + t*s*p5 + t*s^2*p6, numerator;

/* A little linear algebra reveals that 2*p1 + p4 belongs to Q[x,y,z,w], so this
polynomial vanishes on Z + (rhobar)Z (see the paper) */

testNumerator := 2*p1 + p4;
// print testNumerator;

// In general, for a fixed p, hdenom is as follows:

Newhdenom := z^6 - 6*p*z^5*y^2 - 24*p*z^5*y*x - 6*p*z^5*x^2 + 36*p^2*z^4*y^4 + 78*p^2*z^4*y^3*x +
132*p^2*z^4*y^2*x^2 +
   78*p^2*z^4*y*x^3 + 36*p^2*z^4*x^4 + 8*p^3*z^3*y^6 - 60*p^3*z^3*y^5*x - 168*p^3*z^3*y^4*x^2 -
276*p^3*z^3*y^3*x^3
    - 168*p^3*z^3*y^2*x^4 - 60*p^3*z^3*y*x^5 + 8*p^3*z^3*x^6 - 24*p^4*z^2*y^8 - 24*p^4*z^2*y^7*x +
    156*p^4*z^2*y^6*x^2 + 396*p^4*z^2*y^5*x^3 + 540*p^4*z^2*y^4*x^4 + 396*p^4*z^2*y^3*x^5 + 156*p^4*z^2*y^2*x^6
    - 24*p^4*z^2*y*x^7 - 24*p^4*z^2*x^8 + 24*p^5*z*y^9*x + 24*p^5*z*y^8*x^2 - 120*p^5*z*y^7*x^3 -
    324*p^5*z*y^6*x^4 - 432*p^5*z*y^5*x^5 - 324*p^5*z*y^4*x^6 - 120*p^5*z*y^3*x^7 + 24*p^5*z*y^2*x^8 +
    24*p^5*z*y*x^9 + 16*p^6*y^12 + 48*p^6*y^11*x + 48*p^6*y^10*x^2 + 48*p^6*y^9*x^3 + 120*p^6*y^8*x^4 +
    192*p^6*y^7*x^5 + 212*p^6*y^6*x^6 + 192*p^6*y^5*x^7 + 120*p^6*y^4*x^8 + 48*p^6*y^3*x^9 +
    48*p^6*y^2*x^10 + 48*p^6*y*x^11 + 16*p^6*x^12;

TheFunction := 12*z^6 + -72*p*z^5*y^2 + -192*p*z^5*y*x + -48*p*z^5*x^2 + 
300*p^2*z^4*y^4 + 
    600*p^2*z^4*y^3*x + 576*p^2*z^4*y^2*x^2 + 408*p^2*z^4*y*x^3 + 156*p^2*z^4*x^4 + 
    -288*p^3*z^3*y^6 + -720*p^3*z^3*y^5*x + -888*p^3*z^3*y^4*x^2 + -768*p^3*z^3*y^3*x^3
    + -756*p^3*z^3*y^2*x^4 + -264*p^3*z^3*y*x^5 + -204*p^3*z^3*x^6 + 144*p^4*z^2*y^8 + 
    456*p^4*z^2*y^7*x + 1032*p^4*z^2*y^6*x^2 + 1080*p^4*z^2*y^5*x^3 + 
    756*p^4*z^2*y^4*x^4 + 864*p^4*z^2*y^3*x^5 + 684*p^4*z^2*y^2*x^6 + 
    456*p^4*z^2*y*x^7 + -48*p^4*z^2*x^8 + 192*p^5*z*y^10 + -48*p^5*z*y^9*x + 
    -720*p^5*z*y^8*x^2 + -1104*p^5*z*y^7*x^3 + -600*p^5*z*y^6*x^4 + 
    -216*p^5*z*y^5*x^5 + -240*p^5*z*y^4*x^6 + -480*p^5*z*y^3*x^7 + 
    -504*p^5*z*y^2*x^8 + -24*p^5*z*y*x^9 + 48*p^5*z*x^10 + -192*p^6*y^12 + 
    -288*p^6*y^11*x + 192*p^6*y^10*x^2 + 528*p^6*y^9*x^3 + 432*p^6*y^8*x^4 + 
    168*p^6*y^7*x^5 + -192*p^6*y^6*x^6 + -288*p^6*y^5*x^7 + 192*p^6*y^4*x^8 + 
    312*p^6*y^3*x^9 + -48*p^6*y*x^11;

// print TheFunction;

if testNumerator ne TheFunction then
print "mmmmmm";
end if;

/* Now we check for differences in local Norm Residue symbols for the 
quaternion algebra (p,testNumerator/Newhdenom) at the points [w,z,y,x] = [0:-p:1:0]
and [0:-p:0:1]. If there is a discrepancy then we'll get a counter-example to weak approximation. */

if NormResidueSymbol(Rationals()!p,Rationals()!Evaluate(testNumerator/Newhdenom,[0,-p,1,0]),p) 
	ne NormResidueSymbol(Rationals()!p,Rationals()!Evaluate(testNumerator/Newhdenom,[0,-p,0,1]),p) then
	print p, "gives a counterexample to Weak Approximation!";
	Append(~ListOfPrimes,p);
end if;

if NormResidueSymbol(Rationals()!p,Rationals()!Evaluate(testNumerator/Newhdenom,[0,-p,1,0]),3)
        ne NormResidueSymbol(Rationals()!p,Rationals()!Evaluate(testNumerator/Newhdenom,[0,-p,0,1]),3) 
then
        print p, "gives a counterexample to Weak Approximation!";
        Append(~ListOfPrimes,p);
end if;

if NormResidueSymbol(Rationals()!p,Rationals()!Evaluate(testNumerator/Newhdenom,[0,-p,1,0]),2)
        ne NormResidueSymbol(Rationals()!p,Rationals()!Evaluate(testNumerator/Newhdenom,[0,-p,0,1]),2) 
then
        print p, "gives a counterexample to Weak Approximation!";
        Append(~ListOfPrimes,p);
end if;

end if;

end for;