\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\
\\
\\  This file contains the formulas to accompany "Rank frequencies for 
\\  quadratic twists of elliptic curves", by K. Rubin and A. Silverberg, 
\\  in a form suitable for input into PARI. 
\\
\\  After reading this file into PARI, typing "cor42" sets the values 
\\  of f, g, p1, and p2 so that p1 and p2 are independent points of 
\\  infinite order on the twist g(u) y^2 = f(x) as in Corollary 4.2.
\\  Analogously for cor44, thm46, thm51, thm52a, thm52b, thm53, and thm55.
\\
\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\

\\ Corollary 4.2

cor42() = \
    f = x^3 + a*x^2 + b*x; \
    g = -a*b*(b^2 + u^2)*(u^4 + (2*b^2 - a^2*b)*u^2 + b^4); \
    p1 = [-((b^2 + u^2)/(a*b)), 1/(a^2*b^2)]; \
    p2 = [-((b*(b^2 + u^2))/(a*u^2)), b/(a^2*u^3)]; \
    print("Corollary 4.2");


\\ Corollary 4.4

cor44() = \
    f = x^3 + b^2/(4*c)*x^2 + b*x + c; \
    g = -b*c*(2*u^6 + (18*c^2 - b^3)*u^4 + (54*c^4 + 2*b^3*c^2)*u^2 + 54*c^6 - \
          b^3*c^4); \
    p1 = [-((3*c^2 + u^2)/(2*b*c)), 1/(4*b^2*c^2)]; \
    p2 = [(b^3*c^6 - 54*c^8 + (-3*b^3*c^4 - 54*c^6)*u^2 + (3*b^3*c^2 - 18*c^4)*u^4 + \
          (-b^3 - 2*c^2)*u^6)/(4*b*c*u^2*(3*c^2 + u^2)^2), \
         (b^3*c^6 - 54*c^8 + (b^3*c^4 - 54*c^6)*u^2 + (-5*b^3*c^2 - 18*c^4)*u^4 + \
           (3*b^3 - 2*c^2)*u^6)/(8*b^2*c*u^3*(3*c^2 + u^2)^3)]; \
    print("Corollary 4.4");


\\ Theorem 4.6

thm46() = \
    f = x^3 + a*x + b; \
    g = -a*b*(b^2*(u^4 + u^2 + 1)^3 + a^3*u^4*(u^2 + 1)^2)*(u^2 + 1); \
    p1 = [-((b*(-1 + u^6))/(a*(-1 + u^4))), 1/(a^2*(1 + u^2)^2)]; \
    p2 = [-((b*(-1 + u^6))/(a*u^2*(-1 + u^4))), 1/(a^2*u^3*(1 + u^2)^2)]; \
    print("Theorem 4.6");


\\ Theorem 5.1

thm51() = \
    lambda = -2*a^2; \
    f = x*(x - 1)*(x - lambda); \
    d = lambda*(2*lambda - 1)*u^2 + 2 - lambda; \
    n = lambda^2*(lambda + 1)*(2*lambda - 1)^2*u^4 - \
        4*lambda^2*(lambda - 1)*(2*lambda - 1)*u^3 + \
        2*lambda*(lambda + 1)*(2*lambda^2 - 3*lambda + 2)*u^2 - \
        4*lambda*(lambda - 1)*(lambda - 2)*u + (lambda - 2)^2*(lambda + 1); \
    g = 2*n*(n - 2*d^2)*(n - 2*lambda*d^2); \
    p1 = [n/(2*d^2), 1/(4*d^3)]; \
    p2 = [(lambda^2*(d^2 - \
             4*lambda*u*(u - 1)*(lambda*(2*lambda - 1)*u + 2 - lambda)))/ \
         (lambda*(2*lambda - 1)*u^2 - 2*lambda*(2*lambda - 1)*u + lambda - 2)^2, \
        (a*lambda)/(lambda*(2*lambda - 1)*u^2 - 2*lambda*(2*lambda - 1)*u + \
            lambda - 2)^3]; \
    p3 = [(d^2 + 4*lambda*u*(u - 1)*(lambda*(2*lambda - 1)*u + 2 - lambda))/ \
         (lambda*(lambda*(2*lambda - 1)*u^2 - (2*lambda - 4)*u + lambda - 2)^2), \
        -(a/(lambda^2*(lambda*(2*lambda - 1)*u^2 - (2*lambda - 4)*u + lambda - 2)^3))]; \
    print("Theorem 5.1");


\\ Theorem 5.2(a)

thm52a() = \
    lambda = (1 - a^2)/(a^2 + 2); \
    f = x*(x - 1)*(x - lambda); \
    g = -2*(1 + lambda)*(4*u^2 + 8*u^3 + 4*u^4 - 4*u*lambda - 16*u^2*lambda - \
          20*u^3*lambda - 4*u^4*lambda + lambda^2 + 12*u*lambda^2 + \
          26*u^2*lambda^2 + 8*u^3*lambda^2 + u^4*lambda^2 - 2*lambda^3 - \
          12*u*lambda^3 - 10*u^2*lambda^3 + 3*lambda^4 + 8*u*lambda^4 - \
          2*u^2*lambda^4 - 2*lambda^5 + lambda^6)* \
        (4*u^4 - 8*u^3*lambda - 4*u^4*lambda + lambda^2 + 4*u*lambda^2 + \
          14*u^2*lambda^2 + 20*u^3*lambda^2 + u^4*lambda^2 - 2*lambda^3 - \
          12*u*lambda^3 - 26*u^2*lambda^3 - 8*u^3*lambda^3 + 3*lambda^4 + \
          12*u*lambda^4 + 14*u^2*lambda^4 - 2*lambda^5 - 8*u*lambda^5 + lambda^6)* \
        (2*u^2 + 4*u^3 + 4*u^4 - 2*u*lambda - 8*u^2*lambda - 14*u^3*lambda - \
          4*u^4*lambda + lambda^2 + 8*u*lambda^2 + 20*u^2*lambda^2 + \
          14*u^3*lambda^2 + u^4*lambda^2 - 2*lambda^3 - 12*u*lambda^3 - \
          18*u^2*lambda^3 - 4*u^3*lambda^3 + 3*lambda^4 + 10*u*lambda^4 + \
          6*u^2*lambda^4 - 2*lambda^5 - 4*u*lambda^5 + lambda^6); \
    p1 = [(2*lambda*(2*u^2 + 4*u^3 + 4*u^4 - 2*u*lambda - 8*u^2*lambda - \
             14*u^3*lambda - 4*u^4*lambda + lambda^2 + 8*u*lambda^2 + \
             20*u^2*lambda^2 + 14*u^3*lambda^2 + u^4*lambda^2 - 2*lambda^3 - \ 
             12*u*lambda^3 - 18*u^2*lambda^3 - 4*u^3*lambda^3 + 3*lambda^4 + \
             10*u*lambda^4 + 6*u^2*lambda^4 - 2*lambda^5 - 4*u*lambda^5 + \
             lambda^6))/ \
         ((1 + lambda)*(2*u^2 + lambda - u^2*lambda - lambda^2 + lambda^3)^2), \ 
        -(((-1 + lambda)*lambda)/ \
           ((1 + lambda)^2*(2*u^2 + lambda - u^2*lambda - lambda^2 + lambda^3)^3))]; \
    p2 = [(lambda*(4*u^2 + 8*u^3 + 4*u^4 - 4*u*lambda - \
             16*u^2*lambda - 20*u^3*lambda - 4*u^4*lambda + lambda^2 + \
             12*u*lambda^2 + 26*u^2*lambda^2 + 8*u^3*lambda^2 + u^4*lambda^2 - \
             2*lambda^3 - 12*u*lambda^3 - 10*u^2*lambda^3 + 3*lambda^4 + \
             8*u*lambda^4 - 2*u^2*lambda^4 - 2*lambda^5 + lambda^6))/ \
         ((-1 + 2*lambda)*(-2*u^2 + lambda + 4*u*lambda + u^2*lambda - \
              lambda^2 - 2*u*lambda^2 + lambda^3)^2), \
        (a*(-1 + lambda)*lambda)/ \
         ((-1 + 2*lambda)^2*(-2*u^2 + lambda + 4*u*lambda + u^2*lambda - \
              lambda^2 - 2*u*lambda^2 + lambda^3)^3)]; \
    p3 = [(lambda*(4*u^4 - 8*u^3*lambda - 4*u^4*lambda + lambda^2 + \
             4*u*lambda^2 + 14*u^2*lambda^2 + 20*u^3*lambda^2 + u^4*lambda^2 - \
             2*lambda^3 - 12*u*lambda^3 - 26*u^2*lambda^3 - 8*u^3*lambda^3 + \
             3*lambda^4 + 12*u*lambda^4 + 14*u^2*lambda^4 - 2*lambda^5 - \
             8*u*lambda^5 + lambda^6))/ \
         ((-1 + 2*lambda)*(-2*u - 2*u^2 + lambda + 2*u*lambda + u^2*lambda - \
              lambda^2 - 2*u*lambda^2 + lambda^3)^2), \
        -((a*(-1 + lambda)*lambda)/ \
           ((-1 + 2*lambda)^2*(-2*u - 2*u^2 + lambda + 2*u*lambda + u^2*lambda - \ 
                lambda^2 - 2*u*lambda^2 + lambda^3)^3))]; \
    print("Theorem 5.2(a)");


\\ Theorem 5.2(b)

thm52b() = \
    lambda = (a*(a - 2))/(a^2 + 1); \
    f = x*(x - 1)*(x - lambda); \
    g = (-2 + a)*a*(1 + 8*a + 34*a^2 + 88*a^3 + 139*a^4 + 124*a^5 + 26*a^6 - \
          12*a^7 + 41*a^8 - 12*a^9 + 4*a^10 - 4*u - 20*a*u - 44*a^2*u - \
          32*a^3*u + 80*a^4*u + 160*a^5*u + 140*a^6*u + 172*a^7*u + 36*a^8*u + \ 
          16*a^10*u + 6*u^2 + 12*a*u^2 - 2*a^2*u^2 - 20*a^3*u^2 - 68*a^4*u^2 +\ 
          20*a^5*u^2 + 114*a^6*u^2 - 4*a^7*u^2 + 230*a^8*u^2 - 56*a^9*u^2 + \ 
          56*a^10*u^2 - 4*u^3 + 4*a*u^3 - 12*a^3*u^3 + 32*a^4*u^3 - 44*a^5*u^3 + \ 
          56*a^6*u^3 - 20*a^7*u^3 + 36*a^8*u^3 + 24*a^9*u^3 + 8*a^10*u^3 + \
          16*a^11*u^3 + u^4 - 4*a*u^4 + 12*a^2*u^4 - 24*a^3*u^4 + 42*a^4*u^4 - \
          56*a^5*u^4 + 68*a^6*u^4 - 64*a^7*u^4 + 57*a^8*u^4 - 36*a^9*u^4 + \
          24*a^10*u^4 - 8*a^11*u^4 + 4*a^12*u^4)* \
        (1 + 8*a + 34*a^2 + 88*a^3 + 139*a^4 + 124*a^5 + 26*a^6 - 12*a^7 + \
          41*a^8 - 12*a^9 + 4*a^10 - 4*u - 12*a*u - 24*a^2*u - 4*a^3*u + \
          44*a^4*u + 40*a^5*u + 140*a^6*u + 52*a^7*u + 72*a^8*u + 28*a^9*u - \
          4*a^10*u + 8*a^11*u + 6*u^2 - 4*a*u^2 + 6*a^2*u^2 - 36*a^3*u^2 + \
          4*a^4*u^2 - 28*a^5*u^2 + 26*a^6*u^2 + 44*a^7*u^2 + 38*a^8*u^2 + \
          48*a^9*u^2 + 20*a^10*u^2 + 8*a^11*u^2 + 4*a^12*u^2 - 4*u^3 + \
          12*a*u^3 - 28*a^2*u^3 + 48*a^3*u^3 - 64*a^4*u^3 + 72*a^5*u^3 - \
          56*a^6*u^3 + 48*a^7*u^3 - 4*a^8*u^3 + 12*a^9*u^3 + 20*a^10*u^3 + \
          8*a^12*u^3 + u^4 - 4*a*u^4 + 12*a^2*u^4 - 24*a^3*u^4 + 42*a^4*u^4 - \
          56*a^5*u^4 + 68*a^6*u^4 - 64*a^7*u^4 + 57*a^8*u^4 - 36*a^9*u^4 + \
          24*a^10*u^4 - 8*a^11*u^4 + 4*a^12*u^4)* \
        (1 + 6*a + 20*a^2 + 36*a^3 + 31*a^4 + 22*a^5 + 56*a^6 + 184*a^7 + \
          117*a^8 - 118*a^9 + 110*a^10 - 32*a^11 + 8*a^12 - 4*u - 12*a*u - \
          28*a^2*u - 16*a^3*u + 20*a^4*u + 36*a^5*u + 184*a^6*u + 92*a^7*u + \
          212*a^8*u + 80*a^9*u + 68*a^10*u + 36*a^11*u - 4*a^12*u + 8*a^13*u + \ 
          6*u^2 + 18*a^2*u^2 - 24*a^3*u^2 + 8*a^4*u^2 - 44*a^5*u^2 + \
          106*a^6*u^2 - 64*a^7*u^2 + 282*a^8*u^2 - 92*a^9*u^2 + 228*a^10*u^2 - \ 
          48*a^11*u^2 + 60*a^12*u^2 + 4*a^14*u^2 - 4*u^3 + 12*a*u^3 - \
          32*a^2*u^3 + 60*a^3*u^3 - 92*a^4*u^3 + 120*a^5*u^3 - 120*a^6*u^3 + \
          120*a^7*u^3 - 60*a^8*u^3 + 60*a^9*u^3 + 16*a^10*u^3 + 12*a^11*u^3 + \
          28*a^12*u^3 + 8*a^14*u^3 + u^4 - 6*a*u^4 + 22*a^2*u^4 - 56*a^3*u^4 + \ 
          114*a^4*u^4 - 188*a^5*u^4 + 264*a^6*u^4 - 312*a^7*u^4 + 321*a^8*u^4 - \ 
          278*a^9*u^4 + 210*a^10*u^4 - 128*a^11*u^4 + 68*a^12*u^4 - 24*a^13*u^4 + \ 
          8*a^14*u^4); \
    p1 = [((1 + a^2)*(1 + 8*a + 34*a^2 + 88*a^3 + 139*a^4 + 124*a^5 + \
             26*a^6 - 12*a^7 + 41*a^8 - 12*a^9 + 4*a^10 - 4*u - 20*a*u - \ 
             44*a^2*u - 32*a^3*u + 80*a^4*u + 160*a^5*u + 140*a^6*u + \
             172*a^7*u + 36*a^8*u + 16*a^10*u + 6*u^2 + 12*a*u^2 - 2*a^2*u^2 - \
             20*a^3*u^2 - 68*a^4*u^2 + 20*a^5*u^2 + 114*a^6*u^2 - 4*a^7*u^2 + \
             230*a^8*u^2 - 56*a^9*u^2 + 56*a^10*u^2 - 4*u^3 + 4*a*u^3 - \
             12*a^3*u^3 + 32*a^4*u^3 - 44*a^5*u^3 + 56*a^6*u^3 - 20*a^7*u^3 + \
             36*a^8*u^3 + 24*a^9*u^3 + 8*a^10*u^3 + 16*a^11*u^3 + u^4 - \
             4*a*u^4 + 12*a^2*u^4 - 24*a^3*u^4 + 42*a^4*u^4 - 56*a^5*u^4 + \
             68*a^6*u^4 - 64*a^7*u^4 + 57*a^8*u^4 - 36*a^9*u^4 + 24*a^10*u^4 - \
             8*a^11*u^4 + 4*a^12*u^4))/ \
         ((-2 + a)*a*(-1 - 4*a - 9*a^2 - 8*a^3 + 3*a^4 - 2*a^5 + u^2 - \
              2*a*u^2 + 4*a^2*u^2 - 4*a^3*u^2 + 5*a^4*u^2 - 2*a^5*u^2 + \
              2*a^6*u^2)^2), (1 + 2*a)/ \
         ((-2 + a)^2*a^2*(-1 - 4*a - 9*a^2 - 8*a^3 + 3*a^4 - 2*a^5 + u^2 - \
              2*a*u^2 + 4*a^2*u^2 - 4*a^3*u^2 + 5*a^4*u^2 - 2*a^5*u^2 + \
              2*a^6*u^2)^3)]; \
    p2 = [((-2 + a)*a*(1 + a^2)*(1 + 8*a + 34*a^2 + 88*a^3 + 139*a^4 + 124*a^5 + \
             26*a^6 - 12*a^7 + 41*a^8 - 12*a^9 + 4*a^10 - 4*u - 20*a*u - \
             44*a^2*u - 32*a^3*u + 80*a^4*u + 160*a^5*u + 140*a^6*u + \
             172*a^7*u + 36*a^8*u + 16*a^10*u + 6*u^2 + 12*a*u^2 - 2*a^2*u^2 - \
             20*a^3*u^2 - 68*a^4*u^2 + 20*a^5*u^2 + 114*a^6*u^2 - 4*a^7*u^2 + \
             230*a^8*u^2 - 56*a^9*u^2 + 56*a^10*u^2 - 4*u^3 + 4*a*u^3 - \
             12*a^3*u^3 + 32*a^4*u^3 - 44*a^5*u^3 + 56*a^6*u^3 - 20*a^7*u^3 + \
             36*a^8*u^3 + 24*a^9*u^3 + 8*a^10*u^3 + 16*a^11*u^3 + u^4 - \
             4*a*u^4 + 12*a^2*u^4 - 24*a^3*u^4 + 42*a^4*u^4 - 56*a^5*u^4 + \
             68*a^6*u^4 - 64*a^7*u^4 + 57*a^8*u^4 - 36*a^9*u^4 + 24*a^10*u^4 - \
             8*a^11*u^4 + 4*a^12*u^4))/ \
         (-1 - 3*a - 4*a^2 + 5*a^3 + 20*a^4 + 3*a^5 - a^6 + 2*a^7 + 2*u + \
            8*a^3*u - 6*a^4*u + 16*a^5*u - 4*a^6*u + 8*a^7*u - u^2 + 3*a*u^2 - \
            5*a^2*u^2 + 6*a^3*u^2 - 5*a^4*u^2 + 3*a^5*u^2 + a^6*u^2 + 2*a^8*u^2)^2, \
            ((-2 + a)*a*(1 + 2*a))/ \
         (-1 - 3*a - 4*a^2 + 5*a^3 + 20*a^4 + 3*a^5 - a^6 + 2*a^7 + 2*u + \
            8*a^3*u - 6*a^4*u + 16*a^5*u - 4*a^6*u + 8*a^7*u - u^2 + 3*a*u^2 - \
            5*a^2*u^2 + 6*a^3*u^2 - 5*a^4*u^2 + 3*a^5*u^2 + a^6*u^2 + 2*a^8*u^2)^3]; \
    p3 = [((-2 + a)^2*a^2*(1 + 8*a + 34*a^2 + 88*a^3 + 139*a^4 + 124*a^5 + \
             26*a^6 - 12*a^7 + 41*a^8 - 12*a^9 + 4*a^10 - 4*u - 12*a*u - \
             24*a^2*u - 4*a^3*u + 44*a^4*u + 40*a^5*u + 140*a^6*u + 52*a^7*u + \
             72*a^8*u + 28*a^9*u - 4*a^10*u + 8*a^11*u + 6*u^2 - 4*a*u^2 + \
             6*a^2*u^2 - 36*a^3*u^2 + 4*a^4*u^2 - 28*a^5*u^2 + 26*a^6*u^2 + \
             44*a^7*u^2 + 38*a^8*u^2 + 48*a^9*u^2 + 20*a^10*u^2 + 8*a^11*u^2 + \
             4*a^12*u^2 - 4*u^3 + 12*a*u^3 - 28*a^2*u^3 + 48*a^3*u^3 - \
             64*a^4*u^3 + 72*a^5*u^3 - 56*a^6*u^3 + 48*a^7*u^3 - 4*a^8*u^3 + \
             12*a^9*u^3 + 20*a^10*u^3 + 8*a^12*u^3 + u^4 - 4*a*u^4 + \
             12*a^2*u^4 - 24*a^3*u^4 + 42*a^4*u^4 - 56*a^5*u^4 + 68*a^6*u^4 - \
             64*a^7*u^4 + 57*a^8*u^4 - 36*a^9*u^4 + 24*a^10*u^4 - 8*a^11*u^4 + \
             4*a^12*u^4))/ \
         ((1 + a^2)^2*(1 + 4*a + 9*a^2 + 8*a^3 - 3*a^4 + 2*a^5 - 2*u - 2*a*u - \
              4*a^2*u + 18*a^3*u + 4*a^4*u - 2*a^5*u + 2*a^6*u + u^2 - \
              2*a*u^2 + 4*a^2*u^2 - 4*a^3*u^2 + 5*a^4*u^2 - 2*a^5*u^2 + \
              2*a^6*u^2)^2), -(((-2 + a)*a*(1 + 2*a))/ \
           ((1 + a^2)^3*(1 + 4*a + 9*a^2 + 8*a^3 - 3*a^4 + 2*a^5 - 2*u - 2*a*u - \ 
                4*a^2*u + 18*a^3*u + 4*a^4*u - 2*a^5*u + 2*a^6*u + u^2 - \
                2*a*u^2 + 4*a^2*u^2 - 4*a^3*u^2 + 5*a^4*u^2 - 2*a^5*u^2 + \
                2*a^6*u^2)^3))]; \
    print("Theorem 5.2(b)");


\\ Theorem 5.3

thm53() = \
    f = x*(x - b)*(x - a^2*b); \
    t = (a^2*b*(1 + 2*a^2 + a^4 - 4*a*u + 8*a^2*u - 8*a^3*u + 8*a^4*u - \
          4*a^5*u + 4*u^2 - 8*a*u^2 + 14*a^2*u^2 - 32*a^3*u^2 + \
          40*a^4*u^2 - 24*a^5*u^2 + 6*a^6*u^2 + 16*a^2*u^3 - \
          60*a^3*u^3 + 88*a^4*u^3 - 64*a^5*u^3 + 24*a^6*u^3 - \
          4*a^7*u^3 + 16*a^2*u^4 - 56*a^3*u^4 + 81*a^4*u^4 - \
          64*a^5*u^4 + 30*a^6*u^4 - 8*a^7*u^4 + a^8*u^4))/ \
        (-1 - a^2 - 4*a*u^2 + 7*a^2*u^2 - 4*a^3*u^2 + a^4*u^2)^2; \
    g = -4*b*u*(-a + (-1 + a)^2*u)*(-(-1 + a)*(1 + a^2) + a^2*(4 - 3*a + a^2)*u)* \
        (1 + a - 2*(-1 + a)*a*u + a*(4 - 3*a + a^2)*u^2)* \
        ((1 + a^2)^2 - 2*(-1 + a)^2*a*(1 + a^2)*u + \
          (-1 + a)^2*a*(1 + a)*(4 - 3*a + a^2)*u^2)* \
        ((1 + a^2)^2 - 4*(-1 + a)^2*a*(1 + a^2)*u + \
          2*(-1 + a)^2*(2 + 5*a^2 - 6*a^3 + 3*a^4)*u^2 - \
          4*(-1 + a)^3*a^2*(4 - 3*a + a^2)*u^3 + \
          (-1 + a)^2*a^2*(4 - 3*a + a^2)^2*u^4); \
    p1 = [(a^2*b*(1 + 2*a^2 + a^4 - 4*a*u + 8*a^2*u - 8*a^3*u + \
             8*a^4*u - 4*a^5*u + 4*u^2 - 8*a*u^2 + 14*a^2*u^2 - 32*a^3*u^2 + \
             40*a^4*u^2 - 24*a^5*u^2 + 6*a^6*u^2 + 16*a^2*u^3 - 60*a^3*u^3 + \
             88*a^4*u^3 - 64*a^5*u^3 + 24*a^6*u^3 - 4*a^7*u^3 + 16*a^2*u^4 - \
             56*a^3*u^4 + 81*a^4*u^4 - 64*a^5*u^4 + 30*a^6*u^4 - 8*a^7*u^4 + \
             a^8*u^4))/(-1 - a^2 - 4*a*u^2 + 7*a^2*u^2 - 4*a^3*u^2 + a^4*u^2)^2, \
        ((-1 + a)*a^2*b)/ \
         (-1 - a^2 - 4*a*u^2 + 7*a^2*u^2 - 4*a^3*u^2 + a^4*u^2)^3]; \
    p2 = [(a^2*b*(1 + 2*a^2 + a^4 - 4*a*u + 8*a^2*u - 8*a^3*u + \
             8*a^4*u - 4*a^5*u + 4*u^2 - 8*a*u^2 + 14*a^2*u^2 - 32*a^3*u^2 + \
             40*a^4*u^2 - 24*a^5*u^2 + 6*a^6*u^2 + 16*a^2*u^3 - 60*a^3*u^3 + \
             88*a^4*u^3 - 64*a^5*u^3 + 24*a^6*u^3 - 4*a^7*u^3 + 16*a^2*u^4 - \
             56*a^3*u^4 + 81*a^4*u^4 - 64*a^5*u^4 + 30*a^6*u^4 - 8*a^7*u^4 + \
             a^8*u^4))/ \
         (a + a^3 - 2*u + 4*a*u - 4*a^2*u + 4*a^3*u - 2*a^4*u - 4*a^2*u^2 + \
            7*a^3*u^2 - 4*a^4*u^2 + a^5*u^2)^2, \
        -(((-1 + a)*a^2*b)/ \
           (a + a^3 - 2*u + 4*a*u - 4*a^2*u + 4*a^3*u - 2*a^4*u - 4*a^2*u^2 + \
              7*a^3*u^2 - 4*a^4*u^2 + a^5*u^2)^3)]; \
    p3 = [-((16*a^2*b*u*(-a + u - 2*a*u + a^2*u)* \
             (1 - a + a^2 - a^3 + 4*a^2*u - 3*a^3*u + a^4*u)* \
             (1 + 2*a^2 + a^4 - 4*a*u + 8*a^2*u - 8*a^3*u + 8*a^4*u - 4*a^5*u + \
               4*u^2 - 8*a*u^2 + 14*a^2*u^2 - 32*a^3*u^2 + 40*a^4*u^2 - \
               24*a^5*u^2 + 6*a^6*u^2 + 16*a^2*u^3 - 60*a^3*u^3 + 88*a^4*u^3 - \
               64*a^5*u^3 + 24*a^6*u^3 - 4*a^7*u^3 + 16*a^2*u^4 - 56*a^3*u^4 + \
               81*a^4*u^4 - 64*a^5*u^4 + 30*a^6*u^4 - 8*a^7*u^4 + a^8*u^4))/ \
           ((1 + a + 2*a*u - 2*a^2*u + 4*a*u^2 - 3*a^2*u^2 + a^3*u^2)* \
             (-1 + a - a^2 + a^3 - 8*a^2*u + 6*a^3*u - 2*a^4*u + 4*a*u^2 - \
                11*a^2*u^2 + 11*a^3*u^2 - 5*a^4*u^2 + a^5*u^2)^2* \
             (1 + 2*a^2 + a^4 - 2*a*u + 4*a^2*u - 4*a^3*u + 4*a^4*u - 2*a^5*u + \
               4*a*u^2 - 7*a^2*u^2 + 6*a^4*u^2 - 4*a^5*u^2 + a^6*u^2))), \
        -((2*a^2*(1 + a)*b*(-1 + a - 2*a^2 + 2*a^3 - a^4 + a^5 + 8*a*u - \
               16*a^2*u + 20*a^3*u - 20*a^4*u + 12*a^5*u - 4*a^6*u - 8*u^2 + \
               24*a*u^2 - 38*a^2*u^2 + 78*a^3*u^2 - 96*a^4*u^2 + 72*a^5*u^2 - \
               30*a^6*u^2 + 6*a^7*u^2 - 32*a^2*u^3 + 120*a^3*u^3 - \
               192*a^4*u^3 + 172*a^5*u^3 - 92*a^6*u^3 + 28*a^7*u^3 - \
               4*a^8*u^3 - 16*a^2*u^4 + 72*a^3*u^4 - 137*a^4*u^4 + \
               145*a^5*u^4 - 94*a^6*u^4 + 38*a^7*u^4 - 9*a^8*u^4 + a^9*u^4)* \
             (-1 + a - 2*a^2 + 2*a^3 - a^4 + a^5 - 8*a^2*u + 12*a^3*u - \
               12*a^4*u + 12*a^5*u - 4*a^6*u - 6*a^2*u^2 + 14*a^3*u^2 - \
               48*a^4*u^2 + 56*a^5*u^2 - 30*a^6*u^2 + 6*a^7*u^2 + 32*a^3*u^3 - \
               104*a^4*u^3 + 132*a^5*u^3 - 84*a^6*u^3 + 28*a^7*u^3 - 4*a^8*u^3 - \
               16*a^2*u^4 + 72*a^3*u^4 - 137*a^4*u^4 + 145*a^5*u^4 - \
               94*a^6*u^4 + 38*a^7*u^4 - 9*a^8*u^4 + a^9*u^4))/ \
           ((1 + a + 2*a*u - 2*a^2*u + 4*a*u^2 - 3*a^2*u^2 + a^3*u^2)^2* \
             (-1 + a - a^2 + a^3 - 8*a^2*u + 6*a^3*u - 2*a^4*u + 4*a*u^2 - \
                11*a^2*u^2 + 11*a^3*u^2 - 5*a^4*u^2 + a^5*u^2)^3* \
             (1 + 2*a^2 + a^4 - 2*a*u + 4*a^2*u - 4*a^3*u + 4*a^4*u - 2*a^5*u + \
                4*a*u^2 - 7*a^2*u^2 + 6*a^4*u^2 - 4*a^5*u^2 + a^6*u^2)^2))]; \
    print("Theorem 5.3");


\\ Theorem 5.5

thm55() = \
    f = x^3 - x; \
    g = 6*(u^12 - 33*u^8 - 33*u^4 + 1); \
    p1 = [-((u^4 - 6*u^2 + 1)/(3*(u^2 + 1)^2)), 2/(9*(u^2 + 1)^3)]; \
    p2 = [-((u^4 + 6*u^2 + 1)/(3*(u^2 - 1)^2)), 2/(9*(u^2 - 1)^3)]; \
    p3 = [(u^4 + 1)/(6*u^2), 1/(36*u^3)]; \
    print("Theorem 5.5");

test(p) = (g*p[2]^2==subst(f,x,p[1]));