function data(r,s,t)
    x:=PolynomialRing(FieldOfFractions(Parent(r))).1;
    H1:=x^2 + r*x + t;
    lambda1:=4*s;
    G1:=(-s*t + s - 1)*x^3 + (-3*s*t + 3*s*r)*x^2 + (-3*s*t*r + 3*s*r^2)*x - s*t^2 + s*r^3 + t;
    H2:=x^2 + x + r;
    lambda2:=4*s*t;
    G2:=(-s*t + s + 1)*x^3 + (-3*s*t + 3*s*r)*x^2 + (-3*s*t*r + 3*s*r^2)*x - s*t^2 + s*r^3 - t;
    H3:=s*x^2 + (-s*t + 2*s*r - 1)*x + s*r^2;
    lambda3:=4*t/(s^2*t^2 + 2*s*t + 1);
    G3:=(s^2*t^2 - s^2*t + 2*s*t + s + 1)/(s*t + 1)*x^3 + (3*s^2*t^2 - 3*s^2*t*r + 3*s*t + 3*s*r)/(s*t + 1)*x^2 + (3*s^2*t^2*r - 3*s^2*t*r^2 + 3*s*t*r + 3*s*r^2)/(s*t + 1)*x + (s^2*t^3 - s^2*t*r^3 + 2*s*t^2 + s*r^3 + t)/(s*t + 1);
    H4:=(s*t*r - s*t - s*r^2 + s*r + r)*x^2 + (s*t^2 - s*t*r - s*t - s*r^3 + 2*s*r^2 + t)*x + s*t^2 - s*t*r^2 - s*t*r + s*r^3 + t;
    lambda4:=4*s*t/(s^2*t^4 - 6*s^2*t^3*r + 2*s^2*t^3 + 2*s^2*t^2*r^3 + 9*s^2*t^2*r^2 - 6*s^2*t^2*r + s^2*t^2 - 6*s^2*t*r^4 + 2*s^2*t*r^3 + s^2*r^6 + 2*s*t^3 - 6*s*t^2*r + 2*s*t^2 + 2*s*t*r^3 + t^2);
    G4:=(s^2*t^3 - 3*s^2*t^2*r - s^2*t*r^3 + 6*s^2*t*r^2 - 3*s^2*t*r + s^2*t - s^2*r^3 + 2*s*t^2 - 3*s*t*r + 2*s*t - s*r^3 + t)/(s*t^2 - 3*s*t*r + s*t + s*r^3 + t)*x^3 + (3*s^2*t^3 - 6*s^2*t^2*r^2 - 3*s^2*t^2 + 9*s^2*t*r^3 - 3*s^2*t*r^2 + 3*s^2*t*r - 3*s^2*r^4 + 3*s*t^2 - 6*s*t*r^2 + 3*s*t*r)/(s*t^2 - 3*s*t*r + s*t + s*r^3 + t)*x^2 + (-3*s^2*t^3*r + 6*s^2*t^3 - 9*s^2*t^2*r + 3*s^2*t*r^4 + 3*s^2*t*r^3 + 3*s^2*t*r^2 - 3*s^2*r^5 - 3*s*t^2*r + 6*s*t^2 - 3*s*t*r^2)/(s*t^2 - 3*s*t*r + s*t + s*r^3 + t)*x + (-s^2*t^4 + 3*s^2*t^3*r + s^2*t^3 - 6*s^2*t^2*r^2 + 3*s^2*t*r^4 + s^2*t*r^3 - s^2*r^6 - 2*s*t^3 + 3*s*t^2*r + s*t^2 - 2*s*t*r^3 - t^2)/(s*t^2 - 3*s*t*r + s*t + s*r^3 + t);
    delta1:=s;
    delta2:=t;
    delta3:=s*t + 1;
    delta4:=r^3 - 3*r*t + t^2 + t;
    delta5:=r^3*s - 3*r*s*t + s*t^2 + s*t + t;
    delta6:=r^3*s^2 - 3*r*s^2*t - 3*r*s + s^2*t^2 + s^2*t + 2*s*t + s + 1;
    delta7:=r^3*s^2*t + r^3*s - 3*r*s^2*t^2 - 3*r*s*t + s^2*t^3 + s^2*t^2 + 2*s*t^2 + t;
    Delta:=r^6*s^2 - 6*r^4*s^2*t - 3*r^4*s + 2*r^3*s^2*t^2 + 2*r^3*s^2*t + 3*r^3*s*t + r^3*s + r^3 + 9*r^2*s^2*t^2 + 6*r^2*s*t - 6*r*s^2*t^3 - 6*r*s^2*t^2 - 9*r*s*t^2 - 3*r*s*t - 3*r*t + s^2*t^4 + 2*s^2*t^3 + s^2*t^2 + 2*s*t^3 + 3*s*t^2 + t^2 + t;
    Hhat1:=(s^2*t^3 - 4*s^2*t^2*r + s^2*t^2 + 6*s^2*t*r^2 - 4*s^2*t*r + s^2*t - s^2*r^4 + 2*s*t^2 - 5*s*t*r + s*t + 3*s*r^2 - s*r + t - r)*x^2 + (s^2*t^3 - 3*s^2*t^2*r^2 + s^2*t^2*r - s^2*t^2 + 6*s^2*t*r^3 - 5*s^2*t*r^2 + 2*s^2*t*r - s^2*r^5 + 2*s*t^2 - 4*s*t*r^2 + s*t*r - s*t + 3*s*r^3 - s*r^2 + t - r^2)*x - s^2*t^3*r + s^2*t^3 + 2*s^2*t^2*r^2 - 2*s^2*t^2*r - s^2*t*r^3 + s^2*t*r^2 - 2*s*t^2*r + 2*s*t^2 + 2*s*t*r^2 - 2*s*t*r - t*r + t;
    lhat1:=(4*s^3*t^4 - 24*s^3*t^3*r + 8*s^3*t^3 + 8*s^3*t^2*r^3 + 36*s^3*t^2*r^2 - 24*s^3*t^2*r + 4*s^3*t^2 - 24*s^3*t*r^4 + 8*s^3*t*r^3 + 4*s^3*r^6 + 8*s^2*t^3 - 36*s^2*t^2*r + 12*s^2*t^2 + 12*s^2*t*r^3 + 24*s^2*t*r^2 - 12*s^2*t*r - 12*s^2*r^4 + 4*s^2*r^3 + 4*s*t^2 - 12*s*t*r + 4*s*t + 4*s*r^3)/(s^4*t^4 - 6*s^4*t^3*r + 2*s^4*t^3 + 2*s^4*t^2*r^3 + 9*s^4*t^2*r^2 - 6*s^4*t^2*r + s^4*t^2 - 6*s^4*t*r^4 + 2*s^4*t*r^3 + s^4*r^6 + 4*s^3*t^3 - 18*s^3*t^2*r + 6*s^3*t^2 + 4*s^3*t*r^3 + 18*s^3*t*r^2 - 12*s^3*t*r + 2*s^3*t - 6*s^3*r^4 + 2*s^3*r^3 + 6*s^2*t^2 - 18*s^2*t*r + 6*s^2*t + 2*s^2*r^3 + 9*s^2*r^2 - 6*s^2*r + s^2 + 4*s*t - 6*s*r + 2*s + 1);
    Ghat1:=(-s^5*t^7 + 9*s^5*t^6*r - 2*s^5*t^6 - s^5*t^5*r^3 - 33*s^5*t^5*r^2 + 15*s^5*t^5*r - 2*s^5*t^5 + 6*s^5*t^4*r^4 + 64*s^5*t^4*r^3 - 48*s^5*t^4*r^2 + 15*s^5*t^4*r - 2*s^5*t^4 + s^5*t^3*r^6 - 21*s^5*t^3*r^5 - 54*s^5*t^3*r^4 + 64*s^5*t^3*r^3 - 33*s^5*t^3*r^2 + 9*s^5*t^3*r - s^5*t^3 - 3*s^5*t^2*r^7 + 40*s^5*t^2*r^6 - 21*s^5*t^2*r^5 + 6*s^5*t^2*r^4 - s^5*t^2*r^3 + s^5*t*r^9 - 6*s^5*t*r^8 - 3*s^5*t*r^7 + s^5*t*r^6 + s^5*r^9 - 3*s^4*t^6 + 21*s^4*t^5*r - 3*s^4*t^5 - 2*s^4*t^4*r^3 - 54*s^4*t^4*r^2 + 9*s^4*t^4*r + 2*s^4*t^4 + 3*s^4*t^3*r^4 + 80*s^4*t^3*r^3 - 18*s^4*t^3*r^2 - 6*s^4*t^3*r + s^4*t^3 + 4*s^4*t^2*r^6 - 12*s^4*t^2*r^5 - 60*s^4*t^2*r^4 + 16*s^4*t^2*r^3 + 12*s^4*t^2*r^2 - 6*s^4*t^2*r + s^4*t^2 - 9*s^4*t*r^7 + 26*s^4*t*r^6 + 9*s^4*t*r^5 - 9*s^4*t*r^4 + s^4*t*r^3 + s^4*r^9 - 6*s^4*r^7 + 2*s^4*r^6 - 2*s^3*t^5 + 6*s^3*t^4*r + 4*s^3*t^4 + s^3*t^3*r^3 + 6*s^3*t^3*r^2 - 39*s^3*t^3*r + 12*s^3*t^3 - 18*s^3*t^2*r^4 + 7*s^3*t^2*r^3 + 60*s^3*t^2*r^2 - 33*s^3*t^2*r + 4*s^3*t^2 + 5*s^3*t*r^6 + 18*s^3*t*r^5 - 36*s^3*t*r^4 - 9*s^3*t*r^3 + 15*s^3*t*r^2 - 3*s^3*t*r - 6*s^3*r^7 + 4*s^3*r^6 + 9*s^3*r^5 - 6*s^3*r^4 + s^3*r^3 + 2*s^2*t^4 - 18*s^2*t^3*r + 10*s^2*t^3 + 5*s^2*t^2*r^3 + 42*s^2*t^2*r^2 - 45*s^2*t^2*r + 10*s^2*t^2 - 21*s^2*t*r^4 - 6*s^2*t*r^3 + 30*s^2*t*r^2 - 12*s^2*t*r + s^2*t + 2*s^2*r^6 + 9*s^2*r^5 - 12*s^2*r^4 + 3*s^2*r^3 + 3*s*t^3 - 15*s*t^2*r + 6*s*t^2 + 4*s*t*r^3 + 15*s*t*r^2 - 12*s*t*r + 2*s*t - 6*s*r^4 + 3*s*r^3 + t^2 - 3*t*r + t + r^3)/(s^2*t^2 - 3*s^2*t*r + s^2*t + s^2*r^3 + 2*s*t - 3*s*r + s + 1)*x^3 + (-3*s^5*t^7 + 6*s^5*t^6*r^2 + 18*s^5*t^6*r - 3*s^5*t^6 - 51*s^5*t^5*r^3 - 12*s^5*t^5*r^2 - 3*s^5*t^5*r + 3*s^5*t^5 + 12*s^5*t^4*r^5 + 129*s^5*t^4*r^4 - 72*s^5*t^4*r^3 + 57*s^5*t^4*r^2 - 24*s^5*t^4*r + 3*s^5*t^4 - 57*s^5*t^3*r^6 - 81*s^5*t^3*r^5 + 63*s^5*t^3*r^4 - 48*s^5*t^3*r^3 + 21*s^5*t^3*r^2 - 3*s^5*t^3*r + 6*s^5*t^2*r^8 + 60*s^5*t^2*r^7 - 6*s^5*t^2*r^6 + 6*s^5*t^2*r^5 - 3*s^5*t^2*r^4 - 9*s^5*t*r^9 - 15*s^5*t*r^8 + 3*s^5*t*r^7 + 3*s^5*r^10 - 12*s^4*t^6 + 18*s^4*t^5*r^2 + 66*s^4*t^5*r - 12*s^4*t^5 - 129*s^4*t^4*r^3 - 63*s^4*t^4*r^2 + 3*s^4*t^4*r + 9*s^4*t^4 + 30*s^4*t^3*r^5 + 228*s^4*t^3*r^4 - 27*s^4*t^3*r^3 + 33*s^4*t^3*r^2 - 27*s^4*t^3*r + 3*s^4*t^3 - 87*s^4*t^2*r^6 - 90*s^4*t^2*r^5 - 51*s^4*t^2*r^4 + 63*s^4*t^2*r^3 - 18*s^4*t^2*r^2 + 3*s^4*t^2*r + 6*s^4*t*r^8 + 39*s^4*t*r^7 + 57*s^4*t*r^6 - 33*s^4*t*r^5 + 3*s^4*t*r^4 - 18*s^4*r^8 + 6*s^4*r^7 - 18*s^3*t^5 + 18*s^3*t^4*r^2 + 90*s^3*t^4*r - 18*s^3*t^4 - 108*s^3*t^3*r^3 - 99*s^3*t^3*r^2 + 18*s^3*t^3*r + 9*s^3*t^3 + 24*s^3*t^2*r^5 + 117*s^3*t^2*r^4 + 90*s^3*t^2*r^3 - 51*s^3*t^2*r^2 - 30*s^3*t*r^6 - 54*s^3*t*r^5 - 45*s^3*t*r^4 + 48*s^3*t*r^3 - 9*s^3*t*r^2 + 6*s^3*r^7 + 27*s^3*r^6 - 18*s^3*r^5 + 3*s^3*r^4 - 12*s^2*t^4 + 6*s^2*t^3*r^2 + 54*s^2*t^3*r - 12*s^2*t^3 - 33*s^2*t^2*r^3 - 57*s^2*t^2*r^2 + 15*s^2*t^2*r + 3*s^2*t^2 + 6*s^2*t*r^5 + 21*s^2*t*r^4 + 45*s^2*t*r^3 - 27*s^2*t*r^2 + 3*s^2*t*r - 18*s^2*r^5 + 6*s^2*r^4 - 3*s*t^3 + 12*s*t^2*r - 3*s*t^2 - 3*s*t*r^3 - 9*s*t*r^2 + 3*s*t*r + 3*s*r^4)/(s^2*t^2 - 3*s^2*t*r + s^2*t + s^2*r^3 + 2*s*t - 3*s*r + s + 1)*x^2 + (3*s^5*t^7*r - 6*s^5*t^7 - 18*s^5*t^6*r^2 + 51*s^5*t^6*r - 12*s^5*t^6 + 3*s^5*t^5*r^4 + 12*s^5*t^5*r^3 - 129*s^5*t^5*r^2 + 57*s^5*t^5*r - 6*s^5*t^5 + 3*s^5*t^4*r^5 + 72*s^5*t^4*r^4 + 81*s^5*t^4*r^3 - 60*s^5*t^4*r^2 + 9*s^5*t^4*r - 3*s^5*t^3*r^7 - 57*s^5*t^3*r^6 - 63*s^5*t^3*r^5 + 6*s^5*t^3*r^4 + 15*s^5*t^3*r^3 - 3*s^5*t^3*r^2 + 24*s^5*t^2*r^8 + 48*s^5*t^2*r^7 - 6*s^5*t^2*r^6 - 3*s^5*t^2*r^5 - 3*s^5*t*r^10 - 21*s^5*t*r^9 + 3*s^5*t*r^8 + 3*s^5*r^11 + 12*s^4*t^6*r - 24*s^4*t^6 - 60*s^4*t^5*r^2 + 174*s^4*t^5*r - 42*s^4*t^5 + 15*s^4*t^4*r^4 - 3*s^4*t^4*r^3 - 321*s^4*t^4*r^2 + 141*s^4*t^4*r - 12*s^4*t^4 - 6*s^4*t^3*r^5 + 249*s^4*t^3*r^4 + 69*s^4*t^3*r^3 - 81*s^4*t^3*r^2 + 9*s^4*t^3*r - 3*s^4*t^2*r^7 - 102*s^4*t^2*r^6 - 159*s^4*t^2*r^5 + 99*s^4*t^2*r^4 - 18*s^4*t^2*r^3 + 3*s^4*t^2*r^2 + 21*s^4*t*r^8 + 87*s^4*t*r^7 - 39*s^4*t*r^6 + 3*s^4*t*r^5 - 18*s^4*r^9 + 6*s^4*r^8 + 18*s^3*t^5*r - 36*s^3*t^5 - 72*s^3*t^4*r^2 + 216*s^3*t^4*r - 54*s^3*t^4 + 24*s^3*t^3*r^4 - 57*s^3*t^3*r^3 - 252*s^3*t^3*r^2 + 111*s^3*t^3*r - 6*s^3*t^3 - 9*s^3*t^2*r^5 + 252*s^3*t^2*r^4 - 45*s^3*t^2*r^3 - 18*s^3*t^2*r^2 - 72*s^3*t*r^6 - 63*s^3*t*r^5 + 54*s^3*t*r^4 - 9*s^3*t*r^3 + 6*s^3*r^8 + 27*s^3*r^7 - 18*s^3*r^6 + 3*s^3*r^5 + 12*s^2*t^4*r - 24*s^2*t^4 - 36*s^2*t^3*r^2 + 114*s^2*t^3*r - 30*s^2*t^3 + 15*s^2*t^2*r^4 - 57*s^2*t^2*r^3 - 57*s^2*t^2*r^2 + 27*s^2*t^2*r + 3*s^2*t*r^5 + 75*s^2*t*r^4 - 33*s^2*t*r^3 + 3*s^2*t*r^2 - 18*s^2*r^6 + 6*s^2*r^5 + 3*s*t^3*r - 6*s*t^3 - 6*s*t^2*r^2 + 21*s*t^2*r - 6*s*t^2 + 3*s*t*r^4 - 15*s*t*r^3 + 3*s*t*r^2 + 3*s*r^5)/(s^2*t^2 - 3*s^2*t*r + s^2*t + s^2*r^3 + 2*s*t - 3*s*r + s + 1)*x + (s^5*t^8 - 9*s^5*t^7*r + s^5*t^7 + 2*s^5*t^6*r^3 + 33*s^5*t^6*r^2 - 6*s^5*t^6*r - s^5*t^6 - 15*s^5*t^5*r^4 - 64*s^5*t^5*r^3 + 21*s^5*t^5*r^2 + 3*s^5*t^5*r - s^5*t^5 + 2*s^5*t^4*r^6 + 48*s^5*t^4*r^5 + 54*s^5*t^4*r^4 - 40*s^5*t^4*r^3 + 6*s^5*t^4*r^2 - 15*s^5*t^3*r^7 - 64*s^5*t^3*r^6 + 21*s^5*t^3*r^5 + 3*s^5*t^3*r^4 - s^5*t^3*r^3 + 2*s^5*t^2*r^9 + 33*s^5*t^2*r^8 - 6*s^5*t^2*r^7 - s^5*t^2*r^6 - 9*s^5*t*r^10 + s^5*t*r^9 + s^5*r^12 + 5*s^4*t^7 - 39*s^4*t^6*r + 5*s^4*t^6 + 10*s^4*t^5*r^3 + 114*s^4*t^5*r^2 - 21*s^4*t^5*r - 4*s^4*t^5 - 57*s^4*t^4*r^4 - 161*s^4*t^4*r^3 + 51*s^4*t^4*r^2 + 9*s^4*t^4*r - 2*s^4*t^4 + 7*s^4*t^3*r^6 + 129*s^4*t^3*r^5 + 72*s^4*t^3*r^4 - 59*s^4*t^3*r^3 + 6*s^4*t^3*r^2 - 36*s^4*t^2*r^7 - 97*s^4*t^2*r^6 + 48*s^4*t^2*r^5 - 6*s^4*t^2*r^4 + s^4*t^2*r^3 + 4*s^4*t*r^9 + 39*s^4*t*r^8 - 15*s^4*t*r^7 + s^4*t*r^6 - 6*s^4*r^10 + 2*s^4*r^9 + 10*s^3*t^6 - 66*s^3*t^5*r + 10*s^3*t^5 + 19*s^3*t^4*r^3 + 144*s^3*t^4*r^2 - 27*s^3*t^4*r - 6*s^3*t^4 - 78*s^3*t^3*r^4 - 129*s^3*t^3*r^3 + 39*s^3*t^3*r^2 + 9*s^3*t^3*r - s^3*t^3 + 9*s^3*t^2*r^6 + 111*s^3*t^2*r^5 + 9*s^3*t^2*r^4 - 19*s^3*t^2*r^3 - 27*s^3*t*r^7 - 31*s^3*t*r^6 + 21*s^3*t*r^5 - 3*s^3*t*r^4 + 2*s^3*r^9 + 9*s^3*r^8 - 6*s^3*r^7 + s^3*r^6 + 10*s^2*t^5 - 54*s^2*t^4*r + 10*s^2*t^4 + 17*s^2*t^3*r^3 + 78*s^2*t^3*r^2 - 15*s^2*t^3*r - 4*s^2*t^3 - 45*s^2*t^2*r^4 - 31*s^2*t^2*r^3 + 9*s^2*t^2*r^2 + 3*s^2*t^2*r + 5*s^2*t*r^6 + 30*s^2*t*r^5 - 9*s^2*t*r^4 - 6*s^2*r^7 + 2*s^2*r^6 + 5*s*t^4 - 21*s*t^3*r + 5*s*t^3 + 7*s*t^2*r^3 + 15*s*t^2*r^2 - 3*s*t^2*r - s*t^2 - 9*s*t*r^4 + s*t*r^3 + s*r^6 + t^3 - 3*t^2*r + t^2 + t*r^3)/(s^2*t^2 - 3*s^2*t*r + s^2*t + s^2*r^3 + 2*s*t - 3*s*r + s + 1);
    Hhat2:=(-s^2*t^3 + s^2*t^2*r^2 + 2*s^2*t^2*r - 2*s^2*t*r^3 - s^2*t*r^2 + s^2*r^4 - 2*s*t^2 + 2*s*t*r^2 + 2*s*t*r - 2*s*r^3 - t + r^2)*x^2 + (s^2*t^3*r - 3*s^2*t^3 + s^2*t^2*r^2 + 6*s^2*t^2*r - s^2*t^2 - s^2*t*r^4 - 5*s^2*t*r^3 + 2*s^2*r^5 + 2*s*t^2*r - 4*s*t^2 + s*t*r^2 + 3*s*t*r - s*r^4 - s*r^3 + t*r - t)*x + s^2*t^4 - 4*s^2*t^3*r + s^2*t^2*r^3 + 6*s^2*t^2*r^2 - s^2*t^2*r - 4*s^2*t*r^4 + s^2*r^6 + 2*s*t^3 - 5*s*t^2*r + s*t*r^3 + 3*s*t*r^2 - s*r^4 + t^2 - t*r;
    lhat2:=(4*s^3*t^5 - 24*s^3*t^4*r + 8*s^3*t^4 + 8*s^3*t^3*r^3 + 36*s^3*t^3*r^2 - 24*s^3*t^3*r + 4*s^3*t^3 - 24*s^3*t^2*r^4 + 8*s^3*t^2*r^3 + 4*s^3*t*r^6 + 8*s^2*t^4 - 36*s^2*t^3*r + 12*s^2*t^3 + 12*s^2*t^2*r^3 + 24*s^2*t^2*r^2 - 12*s^2*t^2*r - 12*s^2*t*r^4 + 4*s^2*t*r^3 + 4*s*t^3 - 12*s*t^2*r + 4*s*t^2 + 4*s*t*r^3)/(s^4*t^6 - 6*s^4*t^5*r + 2*s^4*t^5 + 2*s^4*t^4*r^3 + 9*s^4*t^4*r^2 - 6*s^4*t^4*r + s^4*t^4 - 6*s^4*t^3*r^4 + 2*s^4*t^3*r^3 + s^4*t^2*r^6 + 4*s^3*t^5 - 18*s^3*t^4*r + 4*s^3*t^4 + 6*s^3*t^3*r^3 + 18*s^3*t^3*r^2 - 6*s^3*t^3*r - 12*s^3*t^2*r^4 + 2*s^3*t^2*r^3 + 2*s^3*t*r^6 + 6*s^2*t^4 - 18*s^2*t^3*r + 2*s^2*t^3 + 6*s^2*t^2*r^3 + 9*s^2*t^2*r^2 - 6*s^2*t*r^4 + s^2*r^6 + 4*s*t^3 - 6*s*t^2*r + 2*s*t*r^3 + t^2);
    Ghat2:=(-s^5*t^8 + 9*s^5*t^7*r - 2*s^5*t^7 - s^5*t^6*r^3 - 33*s^5*t^6*r^2 + 15*s^5*t^6*r - 2*s^5*t^6 + 6*s^5*t^5*r^4 + 64*s^5*t^5*r^3 - 48*s^5*t^5*r^2 + 15*s^5*t^5*r - 2*s^5*t^5 + s^5*t^4*r^6 - 21*s^5*t^4*r^5 - 54*s^5*t^4*r^4 + 64*s^5*t^4*r^3 - 33*s^5*t^4*r^2 + 9*s^5*t^4*r - s^5*t^4 - 3*s^5*t^3*r^7 + 40*s^5*t^3*r^6 - 21*s^5*t^3*r^5 + 6*s^5*t^3*r^4 - s^5*t^3*r^3 + s^5*t^2*r^9 - 6*s^5*t^2*r^8 - 3*s^5*t^2*r^7 + s^5*t^2*r^6 + s^5*t*r^9 - 5*s^4*t^7 + 39*s^4*t^6*r - 10*s^4*t^6 - 5*s^4*t^5*r^3 - 114*s^4*t^5*r^2 + 57*s^4*t^5*r - 7*s^4*t^5 + 21*s^4*t^4*r^4 + 161*s^4*t^4*r^3 - 129*s^4*t^4*r^2 + 36*s^4*t^4*r - 4*s^4*t^4 + 4*s^4*t^3*r^6 - 51*s^4*t^3*r^5 - 72*s^4*t^3*r^4 + 97*s^4*t^3*r^3 - 39*s^4*t^3*r^2 + 6*s^4*t^3*r - 9*s^4*t^2*r^7 + 59*s^4*t^2*r^6 - 48*s^4*t^2*r^5 + 15*s^4*t^2*r^4 - 2*s^4*t^2*r^3 + 2*s^4*t*r^9 - 6*s^4*t*r^8 + 6*s^4*t*r^7 - s^4*t*r^6 - s^4*r^9 - 10*s^3*t^6 + 66*s^3*t^5*r - 19*s^3*t^5 - 10*s^3*t^4*r^3 - 144*s^3*t^4*r^2 + 78*s^3*t^4*r - 9*s^3*t^4 + 27*s^3*t^3*r^4 + 129*s^3*t^3*r^3 - 111*s^3*t^3*r^2 + 27*s^3*t^3*r - 2*s^3*t^3 + 6*s^3*t^2*r^6 - 39*s^3*t^2*r^5 - 9*s^3*t^2*r^4 + 31*s^3*t^2*r^3 - 9*s^3*t^2*r^2 - 9*s^3*t*r^7 + 19*s^3*t*r^6 - 21*s^3*t*r^5 + 6*s^3*t*r^4 + s^3*r^9 + 3*s^3*r^7 - s^3*r^6 - 10*s^2*t^5 + 54*s^2*t^4*r - 17*s^2*t^4 - 10*s^2*t^3*r^3 - 78*s^2*t^3*r^2 + 45*s^2*t^3*r - 5*s^2*t^3 + 15*s^2*t^2*r^4 + 31*s^2*t^2*r^3 - 30*s^2*t^2*r^2 + 6*s^2*t^2*r + 4*s^2*t*r^6 - 9*s^2*t*r^5 + 9*s^2*t*r^4 - 2*s^2*t*r^3 - 3*s^2*r^7 - 5*s*t^4 + 21*s*t^3*r - 7*s*t^3 - 5*s*t^2*r^3 - 15*s*t^2*r^2 + 9*s*t^2*r - s*t^2 + 3*s*t*r^4 - s*t*r^3 + s*r^6 - t^3 + 3*t^2*r - t^2 - t*r^3)/(s^2*t^3 - 3*s^2*t^2*r + s^2*t^2 + s^2*t*r^3 + 2*s*t^2 - 3*s*t*r + s*r^3 + t)*x^3 + (-3*s^5*t^8 + 6*s^5*t^7*r^2 + 18*s^5*t^7*r - 3*s^5*t^7 - 51*s^5*t^6*r^3 - 12*s^5*t^6*r^2 - 3*s^5*t^6*r + 3*s^5*t^6 + 12*s^5*t^5*r^5 + 129*s^5*t^5*r^4 - 72*s^5*t^5*r^3 + 57*s^5*t^5*r^2 - 24*s^5*t^5*r + 3*s^5*t^5 - 57*s^5*t^4*r^6 - 81*s^5*t^4*r^5 + 63*s^5*t^4*r^4 - 48*s^5*t^4*r^3 + 21*s^5*t^4*r^2 - 3*s^5*t^4*r + 6*s^5*t^3*r^8 + 60*s^5*t^3*r^7 - 6*s^5*t^3*r^6 + 6*s^5*t^3*r^5 - 3*s^5*t^3*r^4 - 9*s^5*t^2*r^9 - 15*s^5*t^2*r^8 + 3*s^5*t^2*r^7 + 3*s^5*t*r^10 - 12*s^4*t^7 + 24*s^4*t^6*r^2 + 60*s^4*t^6*r - 15*s^4*t^6 - 174*s^4*t^5*r^3 + 3*s^4*t^5*r^2 + 6*s^4*t^5*r + 3*s^4*t^5 + 42*s^4*t^4*r^5 + 321*s^4*t^4*r^4 - 249*s^4*t^4*r^3 + 102*s^4*t^4*r^2 - 21*s^4*t^4*r - 141*s^4*t^3*r^6 - 69*s^4*t^3*r^5 + 159*s^4*t^3*r^4 - 87*s^4*t^3*r^3 + 18*s^4*t^3*r^2 + 12*s^4*t^2*r^8 + 81*s^4*t^2*r^7 - 99*s^4*t^2*r^6 + 39*s^4*t^2*r^5 - 6*s^4*t^2*r^4 - 9*s^4*t*r^9 + 18*s^4*t*r^8 - 3*s^4*t*r^7 - 3*s^4*r^10 - 18*s^3*t^6 + 36*s^3*t^5*r^2 + 72*s^3*t^5*r - 24*s^3*t^5 - 216*s^3*t^4*r^3 + 57*s^3*t^4*r^2 + 9*s^3*t^4*r + 54*s^3*t^3*r^5 + 252*s^3*t^3*r^4 - 252*s^3*t^3*r^3 + 72*s^3*t^3*r^2 - 6*s^3*t^3*r - 111*s^3*t^2*r^6 + 45*s^3*t^2*r^5 + 63*s^3*t^2*r^4 - 27*s^3*t^2*r^3 + 6*s^3*t*r^8 + 18*s^3*t*r^7 - 54*s^3*t*r^6 + 18*s^3*t*r^5 + 9*s^3*r^8 - 3*s^3*r^7 - 12*s^2*t^5 + 24*s^2*t^4*r^2 + 36*s^2*t^4*r - 15*s^2*t^4 - 114*s^2*t^3*r^3 + 57*s^2*t^3*r^2 - 3*s^2*t^3*r + 30*s^2*t^2*r^5 + 57*s^2*t^2*r^4 - 75*s^2*t^2*r^3 + 18*s^2*t^2*r^2 - 27*s^2*t*r^6 + 33*s^2*t*r^5 - 6*s^2*t*r^4 - 3*s^2*r^7 - 3*s*t^4 + 6*s*t^3*r^2 + 6*s*t^3*r - 3*s*t^3 - 21*s*t^2*r^3 + 15*s*t^2*r^2 - 3*s*t^2*r + 6*s*t*r^5 - 3*s*t*r^4)/(s^2*t^3 - 3*s^2*t^2*r + s^2*t^2 + s^2*t*r^3 + 2*s*t^2 - 3*s*t*r + s*r^3 + t)*x^2 + (3*s^5*t^8*r - 6*s^5*t^8 - 18*s^5*t^7*r^2 + 51*s^5*t^7*r - 12*s^5*t^7 + 3*s^5*t^6*r^4 + 12*s^5*t^6*r^3 - 129*s^5*t^6*r^2 + 57*s^5*t^6*r - 6*s^5*t^6 + 3*s^5*t^5*r^5 + 72*s^5*t^5*r^4 + 81*s^5*t^5*r^3 - 60*s^5*t^5*r^2 + 9*s^5*t^5*r - 3*s^5*t^4*r^7 - 57*s^5*t^4*r^6 - 63*s^5*t^4*r^5 + 6*s^5*t^4*r^4 + 15*s^5*t^4*r^3 - 3*s^5*t^4*r^2 + 24*s^5*t^3*r^8 + 48*s^5*t^3*r^7 - 6*s^5*t^3*r^6 - 3*s^5*t^3*r^5 - 3*s^5*t^2*r^10 - 21*s^5*t^2*r^9 + 3*s^5*t^2*r^8 + 3*s^5*t*r^11 + 12*s^4*t^7*r - 18*s^4*t^7 - 66*s^4*t^6*r^2 + 129*s^4*t^6*r - 30*s^4*t^6 + 12*s^4*t^5*r^4 + 63*s^4*t^5*r^3 - 228*s^4*t^5*r^2 + 87*s^4*t^5*r - 6*s^4*t^5 - 3*s^4*t^4*r^5 + 27*s^4*t^4*r^4 + 90*s^4*t^4*r^3 - 39*s^4*t^4*r^2 - 9*s^4*t^3*r^7 - 33*s^4*t^3*r^6 + 51*s^4*t^3*r^5 - 57*s^4*t^3*r^4 + 18*s^4*t^3*r^3 + 27*s^4*t^2*r^8 - 63*s^4*t^2*r^7 + 33*s^4*t^2*r^6 - 6*s^4*t^2*r^5 - 3*s^4*t*r^10 + 18*s^4*t*r^9 - 3*s^4*t*r^8 - 3*s^4*r^11 + 18*s^3*t^6*r - 18*s^3*t^6 - 90*s^3*t^5*r^2 + 108*s^3*t^5*r - 24*s^3*t^5 + 18*s^3*t^4*r^4 + 99*s^3*t^4*r^3 - 117*s^3*t^4*r^2 + 30*s^3*t^4*r - 18*s^3*t^3*r^5 - 90*s^3*t^3*r^4 + 54*s^3*t^3*r^3 - 6*s^3*t^3*r^2 - 9*s^3*t^2*r^7 + 51*s^3*t^2*r^6 + 45*s^3*t^2*r^5 - 27*s^3*t^2*r^4 - 48*s^3*t*r^7 + 18*s^3*t*r^6 + 9*s^3*r^9 - 3*s^3*r^8 + 12*s^2*t^5*r - 6*s^2*t^5 - 54*s^2*t^4*r^2 + 33*s^2*t^4*r - 6*s^2*t^4 + 12*s^2*t^3*r^4 + 57*s^2*t^3*r^3 - 21*s^2*t^3*r^2 - 15*s^2*t^2*r^5 - 45*s^2*t^2*r^4 + 18*s^2*t^2*r^3 - 3*s^2*t*r^7 + 27*s^2*t*r^6 - 6*s^2*t*r^5 - 3*s^2*r^8 + 3*s*t^4*r - 12*s*t^3*r^2 + 3*s*t^3*r + 3*s*t^2*r^4 + 9*s*t^2*r^3 - 3*s*t^2*r^2 - 3*s*t*r^5)/(s^2*t^3 - 3*s^2*t^2*r + s^2*t^2 + s^2*t*r^3 + 2*s*t^2 - 3*s*t*r + s*r^3 + t)*x + (s^5*t^9 - 9*s^5*t^8*r + s^5*t^8 + 2*s^5*t^7*r^3 + 33*s^5*t^7*r^2 - 6*s^5*t^7*r - s^5*t^7 - 15*s^5*t^6*r^4 - 64*s^5*t^6*r^3 + 21*s^5*t^6*r^2 + 3*s^5*t^6*r - s^5*t^6 + 2*s^5*t^5*r^6 + 48*s^5*t^5*r^5 + 54*s^5*t^5*r^4 - 40*s^5*t^5*r^3 + 6*s^5*t^5*r^2 - 15*s^5*t^4*r^7 - 64*s^5*t^4*r^6 + 21*s^5*t^4*r^5 + 3*s^5*t^4*r^4 - s^5*t^4*r^3 + 2*s^5*t^3*r^9 + 33*s^5*t^3*r^8 - 6*s^5*t^3*r^7 - s^5*t^3*r^6 - 9*s^5*t^2*r^10 + s^5*t^2*r^9 + s^5*t*r^12 + 3*s^4*t^8 - 21*s^4*t^7*r + 2*s^4*t^7 + 3*s^4*t^6*r^3 + 54*s^4*t^6*r^2 - 3*s^4*t^6*r - 4*s^4*t^6 - 9*s^4*t^5*r^4 - 80*s^4*t^5*r^3 + 12*s^4*t^5*r^2 + 9*s^4*t^5*r - s^4*t^5 - 2*s^4*t^4*r^6 + 18*s^4*t^4*r^5 + 60*s^4*t^4*r^4 - 26*s^4*t^4*r^3 + 6*s^4*t^3*r^7 - 16*s^4*t^3*r^6 - 9*s^4*t^3*r^5 + 6*s^4*t^3*r^4 - s^4*t^2*r^9 - 12*s^4*t^2*r^8 + 9*s^4*t^2*r^7 - 2*s^4*t^2*r^6 + 6*s^4*t*r^10 - s^4*t*r^9 - s^4*r^12 + 2*s^3*t^7 - 6*s^3*t^6*r - s^3*t^6 - 4*s^3*t^5*r^3 - 6*s^3*t^5*r^2 + 18*s^3*t^5*r - 5*s^3*t^5 + 39*s^3*t^4*r^4 - 7*s^3*t^4*r^3 - 18*s^3*t^4*r^2 + 6*s^3*t^4*r - 12*s^3*t^3*r^6 - 60*s^3*t^3*r^5 + 36*s^3*t^3*r^4 - 4*s^3*t^3*r^3 + 33*s^3*t^2*r^7 + 9*s^3*t^2*r^6 - 9*s^3*t^2*r^5 - 4*s^3*t*r^9 - 15*s^3*t*r^8 + 6*s^3*t*r^7 + 3*s^3*r^10 - s^3*r^9 - 2*s^2*t^6 + 18*s^2*t^5*r - 5*s^2*t^5 - 10*s^2*t^4*r^3 - 42*s^2*t^4*r^2 + 21*s^2*t^4*r - 2*s^2*t^4 + 45*s^2*t^3*r^4 + 6*s^2*t^3*r^3 - 9*s^2*t^3*r^2 - 10*s^2*t^2*r^6 - 30*s^2*t^2*r^5 + 12*s^2*t^2*r^4 + 12*s^2*t*r^7 - 3*s^2*t*r^6 - s^2*r^9 - 3*s*t^5 + 15*s*t^4*r - 4*s*t^4 - 6*s*t^3*r^3 - 15*s*t^3*r^2 + 6*s*t^3*r + 12*s*t^2*r^4 - 3*s*t^2*r^3 - 2*s*t*r^6 - t^4 + 3*t^3*r - t^3 - t^2*r^3)/(s^2*t^3 - 3*s^2*t^2*r + s^2*t^2 + s^2*t*r^3 + 2*s*t^2 - 3*s*t*r + s*r^3 + t);
    Hhat3:=(-s*t*r^2 + 2*s*t*r - s*t + s*r^4 - 2*s*r^3 + s*r^2)*x^2 + (s*t^3 - 6*s*t^2*r + 3*s*t^2 + 3*s*t*r^3 + 3*s*t*r^2 - 3*s*t*r - 3*s*r^4 + 2*s*r^3 + t^2 - 3*t*r + t + r^3)*x - s*t^2*r + s*t^2 + 2*s*t*r^3 - 2*s*t*r^2 - s*r^5 + s*r^4;
    lhat3:=(4*s^2*t^5 - 24*s^2*t^4*r + 8*s^2*t^4 + 8*s^2*t^3*r^3 + 36*s^2*t^3*r^2 - 24*s^2*t^3*r + 4*s^2*t^3 - 24*s^2*t^2*r^4 + 8*s^2*t^2*r^3 + 4*s^2*t*r^6 + 8*s*t^4 - 36*s*t^3*r + 12*s*t^3 + 12*s*t^2*r^3 + 24*s*t^2*r^2 - 12*s*t^2*r - 12*s*t*r^4 + 4*s*t*r^3 + 4*t^3 - 12*t^2*r + 4*t^2 + 4*t*r^3)/(t^4 - 6*t^3*r + 2*t^3 + 2*t^2*r^3 + 9*t^2*r^2 - 6*t^2*r + t^2 - 6*t*r^4 + 2*t*r^3 + r^6);
    Ghat3:=(-s^3*t^7 + 9*s^3*t^6*r - 2*s^3*t^6 - s^3*t^5*r^3 - 33*s^3*t^5*r^2 + 15*s^3*t^5*r - 2*s^3*t^5 + 6*s^3*t^4*r^4 + 64*s^3*t^4*r^3 - 48*s^3*t^4*r^2 + 15*s^3*t^4*r - 2*s^3*t^4 + s^3*t^3*r^6 - 21*s^3*t^3*r^5 - 54*s^3*t^3*r^4 + 64*s^3*t^3*r^3 - 33*s^3*t^3*r^2 + 9*s^3*t^3*r - s^3*t^3 - 3*s^3*t^2*r^7 + 40*s^3*t^2*r^6 - 21*s^3*t^2*r^5 + 6*s^3*t^2*r^4 - s^3*t^2*r^3 + s^3*t*r^9 - 6*s^3*t*r^8 - 3*s^3*t*r^7 + s^3*t*r^6 + s^3*r^9 - 3*s^2*t^6 + 24*s^2*t^5*r - 6*s^2*t^5 - 4*s^2*t^4*r^3 - 72*s^2*t^4*r^2 + 36*s^2*t^4*r - 5*s^2*t^4 + 21*s^2*t^3*r^4 + 97*s^2*t^3*r^3 - 81*s^2*t^3*r^2 + 27*s^2*t^3*r - 4*s^2*t^3 - 48*s^2*t^2*r^5 - 18*s^2*t^2*r^4 + 33*s^2*t^2*r^3 - 15*s^2*t^2*r^2 + 3*s^2*t^2*r + 6*s^2*t*r^7 + 19*s^2*t*r^6 - 9*s^2*t*r^5 + 3*s^2*t*r^4 - s^2*t*r^3 - s^2*r^9 - 3*s^2*r^7 + s^2*r^6 - 3*s*t^5 + 21*s*t^4*r - 6*s*t^4 - 5*s*t^3*r^3 - 48*s*t^3*r^2 + 27*s*t^3*r - 4*s*t^3 + 21*s*t^2*r^4 + 31*s*t^2*r^3 - 30*s*t^2*r^2 + 9*s*t^2*r - s*t^2 - 2*s*t*r^6 - 21*s*t*r^5 + 9*s*t*r^4 - s*t*r^3 + 3*s*r^7 - t^4 + 6*t^3*r - 2*t^3 - 2*t^2*r^3 - 9*t^2*r^2 + 6*t^2*r - t^2 + 6*t*r^4 - 2*t*r^3 - r^6)/(t^2 - 3*t*r + t + r^3)*x^3 + (-3*s^3*t^7 + 6*s^3*t^6*r^2 + 18*s^3*t^6*r - 3*s^3*t^6 - 51*s^3*t^5*r^3 - 12*s^3*t^5*r^2 - 3*s^3*t^5*r + 3*s^3*t^5 + 12*s^3*t^4*r^5 + 129*s^3*t^4*r^4 - 72*s^3*t^4*r^3 + 57*s^3*t^4*r^2 - 24*s^3*t^4*r + 3*s^3*t^4 - 57*s^3*t^3*r^6 - 81*s^3*t^3*r^5 + 63*s^3*t^3*r^4 - 48*s^3*t^3*r^3 + 21*s^3*t^3*r^2 - 3*s^3*t^3*r + 6*s^3*t^2*r^8 + 60*s^3*t^2*r^7 - 6*s^3*t^2*r^6 + 6*s^3*t^2*r^5 - 3*s^3*t^2*r^4 - 9*s^3*t*r^9 - 15*s^3*t*r^8 + 3*s^3*t*r^7 + 3*s^3*r^10 - 6*s^2*t^6 + 12*s^2*t^5*r^2 + 27*s^2*t^5*r - 3*s^2*t^5 - 81*s^2*t^4*r^3 + 6*s^2*t^4*r^2 - 24*s^2*t^4*r + 9*s^2*t^4 + 18*s^2*t^3*r^5 + 132*s^2*t^3*r^4 - 66*s^2*t^3*r^3 + 54*s^2*t^3*r^2 - 18*s^2*t^3*r - 45*s^2*t^2*r^6 - 66*s^2*t^2*r^5 + 36*s^2*t^2*r^4 - 24*s^2*t^2*r^3 + 9*s^2*t^2*r^2 + 36*s^2*t*r^7 + 3*s^2*t*r^5 - 3*s^2*t*r^4 - 9*s^2*r^8 + 3*s^2*r^7 - 3*s*t^5 + 6*s*t^4*r^2 + 9*s*t^4*r - 30*s*t^3*r^3 + 9*s*t^3*r^2 - 12*s*t^3*r + 3*s*t^3 + 6*s*t^2*r^5 + 30*s*t^2*r^4 - 15*s*t^2*r^3 + 12*s*t^2*r^2 - 3*s*t^2*r - 9*s*t*r^6 - 6*s*t*r^5 + 3*s*r^7)/(t^2 - 3*t*r + t + r^3)*x^2 + (3*s^3*t^7*r - 6*s^3*t^7 - 18*s^3*t^6*r^2 + 51*s^3*t^6*r - 12*s^3*t^6 + 3*s^3*t^5*r^4 + 12*s^3*t^5*r^3 - 129*s^3*t^5*r^2 + 57*s^3*t^5*r - 6*s^3*t^5 + 3*s^3*t^4*r^5 + 72*s^3*t^4*r^4 + 81*s^3*t^4*r^3 - 60*s^3*t^4*r^2 + 9*s^3*t^4*r - 3*s^3*t^3*r^7 - 57*s^3*t^3*r^6 - 63*s^3*t^3*r^5 + 6*s^3*t^3*r^4 + 15*s^3*t^3*r^3 - 3*s^3*t^3*r^2 + 24*s^3*t^2*r^8 + 48*s^3*t^2*r^7 - 6*s^3*t^2*r^6 - 3*s^3*t^2*r^5 - 3*s^3*t*r^10 - 21*s^3*t*r^9 + 3*s^3*t*r^8 + 3*s^3*r^11 + 6*s^2*t^6*r - 12*s^2*t^6 - 27*s^2*t^5*r^2 + 81*s^2*t^5*r - 18*s^2*t^5 + 3*s^2*t^4*r^4 - 6*s^2*t^4*r^3 - 132*s^2*t^4*r^2 + 45*s^2*t^4*r + 24*s^2*t^3*r^5 + 66*s^2*t^3*r^4 + 66*s^2*t^3*r^3 - 36*s^2*t^3*r^2 - 9*s^2*t^2*r^7 - 54*s^2*t^2*r^6 - 36*s^2*t^2*r^5 + 9*s^2*t^2*r^3 + 18*s^2*t*r^8 + 24*s^2*t*r^7 - 3*s^2*t*r^6 - 3*s^2*t*r^5 - 9*s^2*r^9 + 3*s^2*r^8 + 3*s*t^5*r - 6*s*t^5 - 9*s*t^4*r^2 + 30*s*t^4*r - 6*s*t^4 - 9*s*t^3*r^3 - 30*s*t^3*r^2 + 9*s*t^3*r + 12*s*t^2*r^5 + 15*s*t^2*r^4 + 6*s*t^2*r^3 - 3*s*t^2*r^2 - 3*s*t*r^7 - 12*s*t*r^6 + 3*s*r^8)/(t^2 - 3*t*r + t + r^3)*x + (s^3*t^8 - 9*s^3*t^7*r + s^3*t^7 + 2*s^3*t^6*r^3 + 33*s^3*t^6*r^2 - 6*s^3*t^6*r - s^3*t^6 - 15*s^3*t^5*r^4 - 64*s^3*t^5*r^3 + 21*s^3*t^5*r^2 + 3*s^3*t^5*r - s^3*t^5 + 2*s^3*t^4*r^6 + 48*s^3*t^4*r^5 + 54*s^3*t^4*r^4 - 40*s^3*t^4*r^3 + 6*s^3*t^4*r^2 - 15*s^3*t^3*r^7 - 64*s^3*t^3*r^6 + 21*s^3*t^3*r^5 + 3*s^3*t^3*r^4 - s^3*t^3*r^3 + 2*s^3*t^2*r^9 + 33*s^3*t^2*r^8 - 6*s^3*t^2*r^7 - s^3*t^2*r^6 - 9*s^3*t*r^10 + s^3*t*r^9 + s^3*r^12 + 3*s^2*t^7 - 24*s^2*t^6*r + 4*s^2*t^6 + 6*s^2*t^5*r^3 + 72*s^2*t^5*r^2 - 21*s^2*t^5*r - 36*s^2*t^4*r^4 - 97*s^2*t^4*r^3 + 48*s^2*t^4*r^2 - 6*s^2*t^4*r + s^2*t^4 + 5*s^2*t^3*r^6 + 81*s^2*t^3*r^5 + 18*s^2*t^3*r^4 - 19*s^2*t^3*r^3 - 27*s^2*t^2*r^7 - 33*s^2*t^2*r^6 + 9*s^2*t^2*r^5 + 3*s^2*t^2*r^4 + 4*s^2*t*r^9 + 15*s^2*t*r^8 - 3*s^2*t*r^7 - s^2*t*r^6 - 3*s^2*r^10 + s^2*r^9 + 3*s*t^6 - 21*s*t^5*r + 5*s*t^5 + 6*s*t^4*r^3 + 48*s*t^4*r^2 - 21*s*t^4*r + 2*s*t^4 - 27*s*t^3*r^4 - 31*s*t^3*r^3 + 21*s*t^3*r^2 - 3*s*t^3*r + 4*s*t^2*r^6 + 30*s*t^2*r^5 - 9*s*t^2*r^4 - 9*s*t*r^7 + s*t*r^6 + s*r^9 + t^5 - 6*t^4*r + 2*t^4 + 2*t^3*r^3 + 9*t^3*r^2 - 6*t^3*r + t^3 - 6*t^2*r^4 + 2*t^2*r^3 + t*r^6)/(t^2 - 3*t*r + t + r^3);
    Hhat4:=(-s*t*r + s*t + s*r^2 - s*r - r + 1)*x^2 + (-s*t^2 + s*t*r + s*t + s*r^3 - 2*s*r^2 - t + r)*x - s*t^2 + s*t*r^2 + s*t*r - s*r^3 - t + r^2;
    lhat4:=4*s^3*t^5 - 24*s^3*t^4*r + 8*s^3*t^4 + 8*s^3*t^3*r^3 + 36*s^3*t^3*r^2 - 24*s^3*t^3*r + 4*s^3*t^3 - 24*s^3*t^2*r^4 + 8*s^3*t^2*r^3 + 4*s^3*t*r^6 + 8*s^2*t^4 - 36*s^2*t^3*r + 12*s^2*t^3 + 12*s^2*t^2*r^3 + 24*s^2*t^2*r^2 - 12*s^2*t^2*r - 12*s^2*t*r^4 + 4*s^2*t*r^3 + 4*s*t^3 - 12*s*t^2*r + 4*s*t^2 + 4*s*t*r^3;
    Ghat4:=(-s^3*t^5 + 6*s^3*t^4*r - s^3*t^4 - 2*s^3*t^3*r^3 - 9*s^3*t^3*r^2 + s^3*t^3 + 6*s^3*t^2*r^4 + 9*s^3*t^2*r^2 - 6*s^3*t^2*r + s^3*t^2 - s^3*t*r^6 - 6*s^3*t*r^4 + 2*s^3*t*r^3 + s^3*r^6 - 3*s^2*t^4 + 15*s^2*t^3*r - 3*s^2*t^3 - 5*s^2*t^2*r^3 - 15*s^2*t^2*r^2 + 2*s^2*t^2 + 9*s^2*t*r^4 + 6*s^2*t*r^2 - 3*s^2*t*r - s^2*r^6 - 3*s^2*r^4 + s^2*r^3 - 3*s*t^3 + 12*s*t^2*r - 3*s*t^2 - 4*s*t*r^3 - 6*s*t*r^2 + s*t + 3*s*r^4 - t^2 + 3*t*r - t - r^3)*x^3 + (-3*s^3*t^5 + 21*s^3*t^4*r - 6*s^3*t^4 - 6*s^3*t^3*r^3 - 45*s^3*t^3*r^2 + 24*s^3*t^3*r - 3*s^3*t^3 + 24*s^3*t^2*r^4 + 21*s^3*t^2*r^3 - 18*s^3*t^2*r^2 + 3*s^3*t^2*r - 3*s^3*t*r^6 - 18*s^3*t*r^5 + 6*s^3*t*r^4 + 3*s^3*r^7 - 6*s^2*t^4 + 33*s^2*t^3*r - 9*s^2*t^3 - 9*s^2*t^2*r^3 - 45*s^2*t^2*r^2 + 18*s^2*t^2*r + 18*s^2*t*r^4 + 15*s^2*t*r^3 - 9*s^2*t*r^2 - 9*s^2*r^5 + 3*s^2*r^4 - 3*s*t^3 + 12*s*t^2*r - 3*s*t^2 - 3*s*t*r^3 - 9*s*t*r^2 + 3*s*t*r + 3*s*r^4)*x^2 + (-3*s^3*t^5*r + 21*s^3*t^4*r^2 - 6*s^3*t^4*r - 6*s^3*t^3*r^4 - 45*s^3*t^3*r^3 + 24*s^3*t^3*r^2 - 3*s^3*t^3*r + 24*s^3*t^2*r^5 + 21*s^3*t^2*r^4 - 18*s^3*t^2*r^3 + 3*s^3*t^2*r^2 - 3*s^3*t*r^7 - 18*s^3*t*r^6 + 6*s^3*t*r^5 + 3*s^3*r^8 - 6*s^2*t^4*r + 33*s^2*t^3*r^2 - 9*s^2*t^3*r - 9*s^2*t^2*r^4 - 45*s^2*t^2*r^3 + 18*s^2*t^2*r^2 + 18*s^2*t*r^5 + 15*s^2*t*r^4 - 9*s^2*t*r^3 - 9*s^2*r^6 + 3*s^2*r^5 - 3*s*t^3*r + 12*s*t^2*r^2 - 3*s*t^2*r - 3*s*t*r^4 - 9*s*t*r^3 + 3*s*t*r^2 + 3*s*r^5)*x - s^3*t^6 + 6*s^3*t^5*r - 2*s^3*t^5 - s^3*t^4*r^3 - 9*s^3*t^4*r^2 + 6*s^3*t^4*r - s^3*t^4 + s^3*t^2*r^6 + 9*s^3*t^2*r^5 - 6*s^3*t^2*r^4 + s^3*t^2*r^3 - 6*s^3*t*r^7 + 2*s^3*t*r^6 + s^3*r^9 - 3*s^2*t^5 + 15*s^2*t^4*r - 5*s^2*t^4 - 3*s^2*t^3*r^3 - 15*s^2*t^3*r^2 + 9*s^2*t^3*r - s^2*t^3 + 2*s^2*t*r^6 + 6*s^2*t*r^5 - 3*s^2*t*r^4 - 3*s^2*r^7 + s^2*r^6 - 3*s*t^4 + 12*s*t^3*r - 4*s*t^3 - 3*s*t^2*r^3 - 6*s*t^2*r^2 + 3*s*t^2*r + s*r^6 - t^3 + 3*t^2*r - t^2 - t*r^3;
    return [<G1,lambda1,H1>,<G2,lambda2,H2>,<G3,lambda3,H3>,<G4,lambda4,H4>,
     <Ghat1,lhat1,Hhat1>,<Ghat2,lhat2,Hhat2>,
     <Ghat3,lhat3,Hhat3>,<Ghat4,lhat4,Hhat4>];
end function;

function CasselsMap(G,rtd)
  //given a polynomial G associated to a 3-torsion point,
  //compute the cassels map associated to y-rtd*G, as a polynomial in
  // h0,h1,g0,g1, describing the divisor class 
  // D= { h0+h1*x+h2*x^2=0, y=g0+g1*x }
  //This describes an appropriate cassels map on the hyperelliptic curve
  // y^2= (rtd)^2*( G^2+lambda* H^3)
  k:=FieldOfFractions(Parent(rtd));
  P<T,N,A,B>:=PolynomialRing(k,4);
  _<t>:=PolynomialRing(P);
  return Resultant(A*t+B - rtd*Evaluate(G,t),t^2-T*t+N);
end function;

function blockdiagonal(L,n)
  zero:=[0: i in L[1]];
  R:=[];
  for j in [1..n] do
     prefix:= &cat[zero: k in [1..j-1]];
     postfix:= &cat[zero: k in [j+1..n]];
     R cat:=[prefix cat l cat postfix : l in L];
   end for;
   return R;
end function;

function pad(L,n,e)
  return L cat [e : i in [1..n-#L]];
end function;
 
function makejaclift(C,basefield)
  F,G:=HyperellipticPolynomials(C);
  assert G eq 0;
  f0,f1,f2,f3,f4,f5,f6:=Explode(pad(Eltseq(F),7,0));
  x4:=PolynomialRing(basefield).1;
  function k4solve(k1in,k2in,k3in)
    error if k1in eq 0, "We're not sure about that";
    kummereq:=(-4*f0*f2 + f1^2)*k1in^4 - 4*f0*f3*k1in^3*k2in - 2*f1*f3*k1in^3*k3in - 4*f0*k1in^3*x4 - 
    4*f0*f4*k1in^2*k2in^2 + (4*f0*f5 - 4*f1*f4)*k1in^2*k2in*k3in - 2*f1*k1in^2*k2in*x4 + 
    (-4*f0*f6 + 2*f1*f5 - 4*f2*f4 + f3^2)*k1in^2*k3in^2 - 4*f2*k1in^2*k3in*x4 - 
    4*f0*f5*k1in*k2in^3 + (8*f0*f6 - 4*f1*f5)*k1in*k2in^2*k3in + (4*f1*f6 - 
    4*f2*f5)*k1in*k2in*k3in^2 - 2*f3*k1in*k2in*k3in*x4 - 2*f3*f5*k1in*k3in^3 - 4*f4*k1in*k3in^2*x4 -
    4*k1in*k3in*x4^2 - 4*f0*f6*k2in^4 - 4*f1*f6*k2in^3*k3in - 4*f2*f6*k2in^2*k3in^2 + 
    k2in^2*x4^2 - 4*f3*f6*k2in*k3in^3 - 2*f5*k2in*k3in^2*x4 + (-4*f4*f6 + f5^2)*k3in^4 - 
    4*f6*k3in^3*x4;
    k4sols:=Roots(kummereq);
    return [[k1in,k2in,k3in,r[1]]:r in k4sols];
  end function;
  function jacsolve(kumpt)
    k1,k2,k3,k4:=Explode(kumpt);
    error if k1 eq 0, "We're not sure about that";
    Asquared := k1^3*k4+f2*k1^4+f3*k1^3*k2+f4*k1^2*k2^2+3*f5*k1^2*k2*k3
    +f5*k1*k2^3-4*f5*k1^2*k2*k3+9*f6*k1^2*k3^2+6*f6*k1*k2^2*k3
    -24*f6*k1^2*k3^2+f6*k2^4-8*f6*k1*k2^2*k3+16*f6*k1^2*k3^2;
    
    AB := -1/2*f3*k1^3*k3-f5*k1*k2^2*k3+4*f5*k1^2*k3^2-7/2*f5*k1^2*k3^2
    -f6*k2^3*k3+4*f6*k1*k2*k3^2-3*f6*k1*k2*k3^2
    +1/2*f1*k1^4-f4*k1^2*k2*k3-1/2*k1^2*k2*k4;

    bl,A:=IsSquare(Asquared);
    if bl and (basefield cmpeq Rationals() or RelativePrecision(A) gt 3) then
      return [[k2/k1,k3/k1,A,AB/A]];
    else
      return [];
    end if;
  end function;
  function fulllift(T,N)
    try
      L:=&cat[jacsolve(kumpoint) : kumpoint in k4solve(1,T,N)];
    catch E
      L:=[];
    end try;
    return L;
  end function;
  return fulllift;
end function;

function isodescent(r,s,t)
  L:=data(r,s,t);
  G1,lambda1,H1:=Explode(L[1]);
  G2,lambda2,H2:=Explode(L[2]);
  C:=HyperellipticCurve(G1^2+lambda1*H1^3);
  J:=Jacobian(C);
  cm:=CasselsMap(G1,1);

  Ghat1,lhat1,Hhat1:=Explode(L[5]);
  Ghat2,lhat2,Hhat2:=Explode(L[6]);
  Chat:=HyperellipticCurve(-3*(Ghat1^2+lhat1*Hhat1^3));
  Jhat:=Jacobian(Chat);

  Hints:=[ [-t[2],t[1]] where t:=Eltseq(l[3]/LeadingCoefficient(l[3])) : l in L];

  dn:=LCM([Denominator(a): a in Eltseq(HyperellipticPolynomials(C))]);
  SC:={p[1]: p in Factorisation(2*3*Numerator(Discriminant(C)))};
  dn:=LCM([Denominator(a): a in Eltseq(HyperellipticPolynomials(Chat))]);
  SChat:={p[1]: p in Factorisation(2*3*Numerator(Discriminant(Chat)))};
  S:=SC meet SChat;

  Q3,toQ3:=pSelmerGroup(3,S);
  HJ:=AbelianGroup(Invariants(Q3) cat Invariants(Q3));

  Qz<zeta>:=NumberField(x^2+x+1);
  OQz:=IntegerRing(Qz);
  rtm3:=2*zeta+1;
  assert rtm3^2 eq -3;
  Qz3,toQz3:=pSelmerGroup(3,Support(&*S * OQz));
  imgN:=[Eltseq(toQ3(Norm(a@@toQz3))): a in OrderedGenerators(Qz3)];
  HJhat:=AbelianGroup(Invariants(Qz3) cat Invariants(Qz3));
  HJhatNorm:=hom<HJhat->HJ|blockdiagonal(imgN,2)>;

  SelJ:=HJ;
  SelJhat:=Kernel(HJhatNorm);



  //////////////////////////////////////////////////////////////////////

  for p in S do

    Qp:=pAdicField(p);
    QtoQp:=Bang(Rationals(),Qp);
    Qp3,toQp3:=pSelmerGroup(3,Qp);
    imgp:=[Eltseq(toQp3(QtoQp(g@@toQ3))):g in OrderedGenerators(Q3)];
    zero:=[0 : c in imgp[1]];
    HJp:=AbelianGroup(Invariants(Qp3) cat Invariants(Qp3));
    HJtoHJp:=hom<HJ->HJp | blockdiagonal(imgp,2)>;

    //we're in the kernel of the norm map, so one place should be enough
    //for our computations.

    P:=Rep(Support(p*OQz));
    QzP,QztoQzP:=Completion(Qz,P);
    Qp`DefaultPrecision:=100;
    QzP`DefaultPrecision:=100;
    rtm3p:=QztoQzP(rtm3);
    QP3,toQP3:=pSelmerGroup(3,QzP);
    imgp:=[Eltseq(toQP3(QztoQzP(g@@toQz3))):g in OrderedGenerators(Qz3)];
    HJhatp:=AbelianGroup(Invariants(QP3) cat Invariants(QP3));
    HJhattoHJhatp:=hom<HJhat->HJhatp | blockdiagonal(imgp,2)>;

    Jp:=BaseChange(J,Qp);
    Kp:=KummerSurface(Jp);
    Cp:=Curve(Jp);

    Jhatp:=BaseChange(Jhat,Qp);
    Khatp:=KummerSurface(Jhatp);
    Chatp:=Curve(Jhatp);

    gamma1:=CasselsMap(G1,Qp!1);
    gamma2:=CasselsMap(G2,Qp!1);
    function CasselsJ(D)
      if Type(D) eq PtHyp then
        x0:=D[1]/D[3];
        y0:=D[2]/D[3]^3;
        L1:=Eltseq(toQp3((y0-Evaluate(G1,x0))^2/-lambda1));
        L2:=Eltseq(toQp3((y0-Evaluate(G2,x0))^2/-lambda2));
      else
        L1:=Eltseq(toQp3(Evaluate(gamma1,D)/-lambda1));
        L2:=Eltseq(toQp3(Evaluate(gamma2,D)/-lambda2));
      end if;
      return HJp!(L1 cat L2);
    end function;
    gamhat1:=CasselsMap(Ghat1,rtm3p);
    gamhat2:=CasselsMap(Ghat2,rtm3p);
    function CasselsJhat(D)
      if Type(D) eq PtHyp then
        x0:=D[1]/D[3];
        y0:=D[2]/D[3]^3;
        L1:=Eltseq(toQp3((y0-rtm3p*Evaluate(Ghat1,x0))^2/3/lhat1));
        L2:=Eltseq(toQp3((y0-rtm3p*Evaluate(Ghat2,x0))^2/3/lhat2));
      else
        L1:=Eltseq(toQP3(Evaluate(gamhat1,D)/3/lhat1));
        L2:=Eltseq(toQP3(Evaluate(gamhat2,D)/3/lhat2));
      end if;
      return HJhatp!(L1 cat L2);
    end function;

    if p eq 3 then
      expected_size:=4;
    elif p mod 3 eq 1 then
      expected_size:=4;
    else
      expected_size:=2;
    end if;

    imgJp:=sub<HJp|>;
    imgJhatp:=sub<HJhatp|>;

    Jplift:=makejaclift(C,Qp);
    Jhatplift:=makejaclift(Chat,Qp);

    extra_tries:=5;
    firstround:=true;
    repeat
      Jpts:={HJp|};Jhatpts:={HJhatp|};
      repeat
        if firstround then
          firstround:=false;
          candJp:=&cat[Jplift(t[1],t[2]): t in Hints];
          candJhatp:=&cat[Jhatplift(t[1],t[2]): t in Hints];
        elif Random([1..4]) eq 1 then
          x0:=p^Random([-3..3])*Random([-p^4..p^4]);
          candJp:=Points(Cp,x0);
          candJhatp:=Points(Chatp,x0);
        else
          slope:=Random([-7..7]);
          T:=Random([-100*p^2..p^2*100])*p^(slope div 2);
          N:=Random([-100*p^2..p^2*100])*p^slope;
          candJp:=Jplift(T,N);
          candJhatp:=Jhatplift(T,N);
        end if;
        for D in candJp do
           try
             Include(~Jpts,CasselsJ(D));
           catch E
             print "ignoring error",E;
           end try;
        end for;
        for D in candJhatp do
           try
             Include(~Jhatpts,CasselsJhat(D));
           catch E
             print "ignoring error",E;
           end try;
        end for;
      until #Jpts gt 0 or #Jhatpts gt 0;
      for img in Jpts do
        if img in imgJp then
          print "found point on J; already in generated image";
        else
          print "new image on J";
          imgJp:=imgJp+sub<HJp|img>;
        end if;
      end for;  
      for img in Jhatpts do
        if img in imgJhatp then
          print "found point on Jhat; already in generated image";
        else
          print "new image on Jhat";
          imgJhatp:=imgJhatp+sub<HJhatp|img>;
        end if;
      end for;
      found_size:=Ngens(imgJp)+Ngens(imgJhatp);
      error if found_size gt expected_size, "Found too big image";
      if found_size eq expected_size then
        extra_tries -:=1;
      end if;
    until extra_tries lt 1;

    SelJ:=SelJ meet (Image(HJtoHJp) meet imgJp)@@HJtoHJp;
    SelJhat:=SelJhat meet (Image(HJhattoHJhatp) meet imgJhatp)@@HJhattoHJhatp;
  end for;
  print "rank SelJ =",Ngens(SelJ);
  print "rank SelJhat =",Ngens(SelJhat);  
  Jlift:=makejaclift(C,Rationals());
  gamma1:=CasselsMap(G1,1);
  gamma2:=CasselsMap(G2,1);
  P:=&cat[ Jlift(-h[2],h[1]) where h:=Eltseq(g/LeadingCoefficient(g)) where g:=l[3] : l in L[1..4]];
  t:=Ngens(sub<SelJ|[HJ!(Eltseq(toQ3(w[1]))cat Eltseq(toQ3(w[2]))): v in P | w[1]*w[2] ne 0
               where w:=[Evaluate(gamma1,v)/-lambda1,Evaluate(gamma2,v)/-lambda2] ]>);
  print "found torsion ",t," room in SelJ",Ngens(SelJ)-t;
  return Ngens(SelJ),Ngens(SelJhat),t;
end function;

function twodescent(r,s,t)
  L:=data(r,s,t);
  G1,lambda1,H1:=Explode(L[1]);
  C:=HyperellipticCurve(G1^2+lambda1*H1^3);
  Cm:=SimplifiedModel(ReducedMinimalWeierstrassModel(C));
  J:=Jacobian(Cm);
  T:=TwoSelmerGroup(J);
  t2:=TwoTorsionSubgroup(J);
  c:=IsOdd(J) select 1 else 0;
  r:=Ngens(T)-Ngens(t2);
  return r, r-c,C;
end function;

function data_is_valid(r,s,t)
  try
    L:=data(r,s,t);
    if forall{Degree(l[1]) eq 3 : l in L} and
       forall{Degree(l[3]) eq 2 : l in L} and
       forall{Degree(l[2]) ne 0 : l in L} then
       C:=HyperellipticCurve(L[1][1]^2+L[1][2]*L[1][3]^3);
       if Genus(C) eq 2 then
         return true;
       end if;
    end if;
    return false;
  catch E
    return false;
  end try;
end function;
    
