# This Maple file gives the computational details for
# Examples 1 and 2 in "Rational Divisors in Rational Divisor Classes"
# by N. Bruin and E.V. Flynn.

# PART 1. For the curve y^2 = -x^6-x^5-2*x^4-2*x^3+x^2-2*x+2
#         we find a rational divisor in a given rational divisor
#         class. The following are adaptable, so that the users
#         can replace the first few lines with their own choices
#         of curves and simialar divisor classes.
#         This part relates to Example 1.

interface(prettyprint,prettyprint=false):

newarr := [
[C0dz, 
 -x^6-x^5-2*x^4-2*x^3+x^2-2*x+2,
 [[1/2+1/5*I*55^(1/2), 7/40+12/25*I*55^(1/2)], 
  [1/2-1/5*I*55^(1/2), 7/40-12/25*I*55^(1/2)]],
 -55, 
 [[1/2+1/5*I*55^(1/2), 7/40+12/25*I*55^(1/2)], 
  [1+I*5^(1/2), 2*I*55^(1/2)],
  [1-I*5^(1/2), 2*I*55^(1/2)]]]
]:

sqrtm55 := sqrt(-55):
sqrtm5 := sqrt(-5):

P1 := [1/2+1/5*sqrtm55, 7/40+12/25*sqrtm55] :
P2 := [1+sqrtm5, 2*sqrtm55]:
P3 := [1-sqrtm5, 2*sqrtm55]:

dofirst := proc(arr)
local arr1,arr2,narr1,narr2,M,pt1,x1,y1,found,WP,
      i,pt,count:
# The following takes a pair of arrays of points:
#  [arr1,arr2] =
# [ [[x1,y2],[x2,y2],..,[xn,yn]], [ whatever ] ]
# and looks to see if any of x2,..,xn are = x1.
# If so, and if any point has x1=xi and y1=-yi
#   then these 2 points are removed from arr1 and [[OO,OO],2]
#   is put at the beginning of arr2.
# Otherwise, all [xi,yi] are found s.t. [x1,y1]=[xi,yi] and
#   all of these are removed from arr1 and [[x1,y1],N] is
#   added at the end of arr2, where N is the multiplicity.
# N.B. pts. at infty are in the form: [infty,sqrt(f6)*infty**3]
# and [infty,-sqrt(f6)*infty**3] for even degree and
# [infty,0] for the odd degree situation.
 
  arr1 := arr[1]:
  arr2 := arr[2]:
  narr1 := arr1:
  narr2 := arr2:
  M := nops(arr1):
  if M = 0 then M := M
  else
    pt1 := arr1[1]:
    x1 := pt1[1]:
    y1 := pt1[2]:
    found := false:
    WP := false:
    i := 1:
    while (i<M) and (not found) do
      i := i+1:
      pt := arr1[i]:
      if (pt[1]=x1) and (pt[2]=-y1) then
        found := true: WP := true:
        narr1 := [op(narr1[2..(i-1)]), op(narr1[(i+1)..M])]:
        narr2 := [[[OO,OO],2], op(narr2)]:
      fi:
    od:
    i := 1:
    count := 1:
    if (not WP) then narr1 := [op(narr1[2..M])]
    fi:
    while (i<M) and (not WP) do
      i := i+1:
      pt := arr1[i]:
      if (pt[1]=x1) and (pt[2]=y1) then
        narr1 := [op(narr1[1..(i-count-1)]),op(narr1[(i-count+1)..(M-count)])]:
        count := count+1:
      fi:
    od:
    if (not WP) then narr2 := [op(narr2),[[x1,y1],count]]
    fi:
  fi:
  [narr1,narr2]
end:
 
convertt := proc(arr) local arrr;
# The following converts an array of points, arr:
# [ [x1,y1],..,[xn,yn] ]
# into an array of points with multiplicities in the style:
# [ [[x1,y1],N1],..,[[x2,y2],Nk] ]
# N.B. [[OO,OO],2]'s are put at the beginning of the outputted array.
# N.B. pts. at infty are in the form: [infty,sqrt(f6)*infty**3]
# and [infty,-sqrt(f6)*infty**3] for even degree and
# [infty,0] for the odd degree situation.
arrr := [arr,[]]:
while not arrr[1]=[] do
  arrr := dofirst(arrr)
od;
arrr[2]:
end:
 
deriv := proc(fnxy,N) local i,arr1,arr2,ffnxy;
# The following takes a polynomial fnxy in x,y
# and finds and expression (in terms of x,y) for
# diff(y(x),x$N).
ffnxy := subs(y=y(x),fnxy);
arr1 := []; arr2 := [];
for i from 0 to N do
  if i>0 then diff(ffnxy,x$i) else ffnxy fi;
  if i>0 then solve(%,diff(y(x),x$i)) else y(x) fi;
  if i>0 then subs(op(arr1),%) else % fi;
  if i>0 then arr1 := [diff(y(x),x$i)=% , op(arr1)]
                  else arr1 := [y(x)=y(x)] fi;
  arr2 := [%%,op(arr2)];
od;
arr2[1];
normal(subs(y(x)=y,%),expanded);
end;

readlib(unassign):

geteqn := proc(ptm,sextic) local pt,mult,x1,y1,cub,vsextic,
                     lhs,rhs,vcub,vx1,vy1,eqnset,i;
# The input is ptm = [[x1,y1],mult] and sextic (in x)
# The output is the set of linear eqns in a,b,c,d
# induced by the condition that the curve y**2 = sextic
# meets y = a*x**3 + b*x**2 + c*x + d at [x1,y1]
# with multiplicity mult.

pt := ptm[1]: mult := ptm[2]: x1 := pt[1]: y1 := pt[2]:
cub := a*x**3 + b*x**2 + c*x + d:
if x1=infty then
  vsextic := expand(x**6*subs(x=1/x,sextic)):
  vcub := expand(x**3*subs(x=1/x,cub)):
  vx1 := 0:
  vy1 := y1/(infty**3):
else
  vsextic := sextic:  vcub := cub:
  vx1 := x1: vy1 := y1:
fi:
eqnset := { subs(x=vx1,y=vy1,vcub) = vy1 }:
for i from 1 to (mult-1) do
  lhs := expand( subs(x=vx1,diff(vcub,x$i)) ):
  rhs := subs(x=vx1,y=vy1, deriv(y**2 - vsextic,i) ):
  radsimp(rhs,'ratdenom'):
  rhs := expand(%):
  eqnset := eqnset union {lhs=rhs}:
od:
eqnset:
end:

geteqns := proc(ptms,sextic) local i,eqnset:
# The input is ptms = [ptms[1],..,ptms[k]] and sextic (in x)
# where each ptms[i] = [[xi,yi],mult_i].
# The output is the set of linear eqns in a,b,c,d
# induced by the condition that the curve y**2 = sextic
# meets y = a*x**3 + b*x**2 + c*x + d at all of the
# [xi,yi] with multiplicity mult_i.

eqnset := {}:
for i from 1 to nops(ptms) do
  eqnset := eqnset union  geteqn(ptms[i] ,sextic) 
od:
eqnset:
end:

getcub := proc(ptms,sextic) global aa,bb,cc,dd,cubb;
# The input is ptms = [ptms[1],..,ptms[k]] and sextic (in x)
# where each ptms[i] = [[xi,yi],mult_i].
# The output is the cubic:
#       y = a*x**3 + b*x**2 + c*x + d
# which meets the curve y**2 = sextic at all the
# [xi,yi] with multiplicity mult_i.

cubb := aa*x**3 + bb*x**2 + cc*x + dd;
geteqns(ptms,sextic):
solve(%,{a,b,c,d}):
subs(a=aa,b=bb,c=cc,d=dd,%);
assign(%):
radsimp(aa,'ratdenom'):
aa := expand(%):
bb := radsimp(bb,'ratdenom'):
bb := expand(%):
cc := radsimp(cc,'ratdenom'):
cc := expand(%):
dd := radsimp(dd,'ratdenom'):
dd := expand(%):
cubb := cubb;
unassign('aa'): unassign('bb'): unassign('cc'): unassign('dd'):
cubb:
end;

deconv := proc(arrr,sextic) local ff6, arrrr;
# The following "deconverts" an array 
#    [ [[OO,OO],2],[[x1,y1],1],[[x2,y2],1] ] 
# into: [ [x1,y1],[x2,y2] ]
# and "deconverts" an array
#    [ [[OO,OO],2],[[x1,y1],2] ]
# into: [ [x1,y1],[x1,y1] ].

ff6 := sqrt(coeff(sextic,x,6)):
if arrr[2][2] = 1 then
  arrrr := [ arrr[2][1] , arrr[3][1] ]
elif arrr[2][1] = [OO,OO] then
  arrrr := [ [infty, ff6*infty**3] , [infty,-ff6*infty**3] ]:
else
  arrrr := [ arrr[2][1] , arrr[2][1] ]
fi:
arrrr:
end:

addd := proc(arr1,arr2,sext) local cub,varsextic,deg,i,varquartic,
                sextic,varquad,x1,y1,x2,y2,summ,arr,arrr;
global a,b,c,d,t;
# N.B. Global: a,b,c,d,t.
# The following inputs: arr1 = [ pt1 , pt2 ] , arr2 = [ pt3 , pt4 ]
# where each pt is of the form [x,y], with a pt. at
# infinity represented by: [infty,sqrtf6*infty**3]
# or [infty,-sqrtf6*infty**3]. The other input is a sextic in x.
# The output is the divisor: [[x1,y1],[x2,y2]] giving
# the sum of arr1 + arr2 on the Jacobian of the
# curve: y**2 = sextic. 

sextic := expand(sext):
arr := [op(arr1),op(arr2)];
arrr := convertt(arr):
if arrr[1] = [[OO,OO],2] then summ := deconv(arrr,sextic)
else
  cub := getcub(arrr,sextic):
  varsextic := expand((cub**2 - sextic)):
  radsimp(%,'ratdenom'):
  varsextic := expand(%):
  deg := degree(varsextic):
  if (not deg = 6) then varsextic := varsextic*(x-infty)**(6-deg) 
  fi:
  varsextic := expand(varsextic);
  varquartic := 1;
  for i from 1 to nops(arr) do 
    varquartic := varquartic*(x-arr[i][1])
  od:
  varquartic := expand(varquartic):
  radsimp(varquartic,'ratdenom'):
  varquartic := expand(%):
  varquad := simplify(varsextic/varquartic);
  [solve(%,x)]:
  x1 := %[1]: x2 := %%[2]:
  if x1 = infty then y1 := -coeff(cub,x,3)*infty**3
  else
    y1 := -subs(x=x1,cub):
    radsimp(%,'ratdenom'):
    y1 := expand(%):
  fi:
  if x2 = infty then y2 := -coeff(cub,x,3)*infty**3 
  else 
    y2 := -subs(x=x2,cub): 
    radsimp(%,'ratdenom'): 
    y2 := expand(%): 
  fi:
  summ := [[x1,y1],[x2,y2]]:
fi:
summ:
end:

#################################################################

# write function in singles of variables in terms of elementary
# symmetric polynomials...
#

gpolyutil[g3onesymmp] := proc(poly::polynom)
local i;
subs(seq(g[`i`]=(-1)^(1+`i`)*p[2-`i`],`i`=0..2),
     gpolyutil[onesymm](subs(seq(x[`i`]=w[`i`],`i`=1..3),poly),3));
end:

gpolyutil[checkonesym] := proc(expr,N)
local i, check, exprr;
check := true;
exprr := expand(expr);
i:=2;
while check = true and i < (N+1) do
  if subs({w[1]=w[i],w[i]=w[1]},exprr) <> exprr then check:=false fi;
  i := i+1 od;
check
end:

gpolyutil[onesymminit] := proc(N)
local
i,x,     # dummy variables
pol1,  # pol1 will be (x-w[1])*...*(x-w[n])
pol2;  # pol2 will be pol1 expanded, with coefficients of x collected.
global
ww, gg, subsgg;
pol1 := 1;
for i from 1 to N do pol1 := pol1 * (x-w[i]) od;
pol2 := collect(pol1,x);
for i from 0 to N-1 do gg[i] := expand(coeff(pol2,x,i)) od;
subsgg := g[0]=gg[0];
for i from 1 to (N-1) do subsgg := subsgg , g[i]=gg[i] od;
ww := [w[1]];
for i from 2 to N do ww := [ op(ww), w[i] ] od;
end:

gpolyutil[gethigh] := proc(expr, N)
local i, exprr , highlex , lead , t;
exprr := collect( expr , ww); # this collects coeffs of w[1],..,w[N]
lead[1] := lcoeff ( exprr , w[1] ,'t');
  if t=1 then highlex :=0
           elif nops(t)=1 then highlex :=[1] else highlex := [op(2,t)] fi;
for i from 2 to N do
               lead[i-1] := collect ( lead[i-1] , w[i] );
               lead[i] := lcoeff( lead[i-1] , w[i] , 't');
               if t=1 then highlex := [op(highlex),0]
               elif nops(t) = 1 then highlex := [op(highlex),1]
               else highlex :=[op(highlex),op(2,t)] fi od;
highlex := [op(highlex), lead[N] ];
end:

gpolyutil[conv] := proc(arr,N)
local i, monom;
monom := arr[N+1]*(g[0]*(-1)^N)^(arr[N]);
for i from 1 to (N-1) do
 monom := monom * (g[i]*(-1)^(N-i))^(arr[N-i] - arr[N-i+1]) od;
monom;
end:

gpolyutil[onesymm] := proc(expr,N)
local remains, gexpr, newg, newgg;
if gpolyutil[checkonesym](expr,N) = false then NOTSYM else
gpolyutil[onesymminit](N);
remains := expand(expr);
gexpr :=0;
while remains <> 0 do
  newg := gpolyutil[conv](gpolyutil[gethigh](remains,N),N);
  gexpr := gexpr + newg;
  newgg := subs(subsgg , newg);
  remains := expand( remains - newgg) od;
gexpr; fi
end:

gpolyutil[onesym] := proc(expr,vlist)
local i, newexpr, N;
N := nops(vlist);
newexpr := expr;
for i from 1 to N do newexpr := subs( vlist[i] = w[i] , newexpr) od;
gpolyutil[onesymm](newexpr,N)
end:

gpolyutil[onesymvar] := proc(expr,vlist)
local i, newexpr, N, subsgtoa;
N := nops(vlist);
subsgtoa := [];
for i from 0 to N do
subsgtoa := [op(subsgtoa),g[i]=(-1)^(i+N)*a[N-i]/a[0]]
od;
newexpr := expr;
for i from 1 to N do newexpr := subs( vlist[i] = w[i] , newexpr) od;
newexpr := gpolyutil[onesymm](newexpr,N);
newexpr := subs(subsgtoa, newexpr);
newexpr := normal(%,expanded);
while not type(newexpr, polynom) do
newexpr := a[0]*newexpr
od;
newexpr;
end:

with(gpolyutil);

#################################################################

i := 1:

sexticx[i] := newarr[i][2]:
sexticu[i] := subs(x=u,%):
ratdiv[i] := newarr[i][5]:
convertt( [ op(ratdiv[i]) ] ):
# The following line gets the cubic s.t. y=cubic passes through [u,v] and
# the 3 points in ratdiv.
cubicx[i] := simplify( getcub(%,sexticx[i]) ):
# The following will be a sextic in x whose roots include u and the x-coords
# of the 3 points in ratdiv. The remaining two roots are what we seek.
difference[i] := factor( numer( simplify( cubicx[i]^2 - sexticx[i] ) ) ):
outratdiv[i] := 1:
for j from 1 to 3 do
  if ratdiv[i][j][1] <> infty then
    outratdiv[i] := outratdiv[i]*(x - ratdiv[i][j][1])
  fi:
od:
difference[i] := simplify( difference[i] / outratdiv[i] ):
cubx[i] := %:
# At this point we have divided out by the known 3 roots and have the res cubic
cubx[i] := radsimp(cubx[i]/subs(x=0,cubx[i]),'ratdenom'):
cubx[i] := expand( numer( cubx[i] ) ):
# At this point, cubx[i] should be over Q(sqrtd). We just need to remove the
# factors not involving x:
relevantcubx := 1:
factorsofcubx := [op(factor(t*cubx[i]))]:
for j from 1 to nops( factorsofcubx ) do
  if has( factorsofcubx[j], x ) then
   relevantcubx := relevantcubx*factorsofcubx[j]
  fi
od:
# The following gives the cubic satidfied by the x-coords of
# the 3 remaining points
cubx[i] := relevantcubx;
# Now we find the quadratic which defines the 3 corresponding y-coords.
gx := g3*x^3+g2*x^2+g1*x+g0:
c1 := subs(x=x1, gx): c2 := subs(x=x2,gx): c3 := subs(x=x3, gx):
simplify(  gpolyutil[onesym]( simplify(
     -(-x1*c3+x3*c1+x2*c3+x1*c2-x2*c1-x3*c2)/
      (x1*x3^2-x2*x3^2+x2*x1^2-x3*x1^2+x3*x2^2-x1*x2^2) ),
       [x1,x2,x3] ));
simplify(  gpolyutil[onesym]( simplify(
       (-x1^2*c3+x1^2*c2+c1*x3^2-c2*x3^2+x2^2*c3-c1*x2^2)/
       (x1*x3^2-x2*x3^2+x2*x1^2-x3*x1^2+x3*x2^2-x1*x2^2) ),
       [x1,x2,x3] ));
simplify(  gpolyutil[onesym]( simplify(
       (-c1*x2*x3^2+x2^2*x3*c1+x1^2*x2*c3-x1^2*x3*c2+c2*x1*x3^2-x2^2*x1*c3)/
       (x1*x3^2-x2*x3^2+x2*x1^2-x3*x1^2+x3*x2^2-x1*x2^2) ),
       [x1,x2,x3] ));
quadgeneral := %%%*x^2 + %%*x + %;
# Gives: quadgeneral := (-g3*g[2]+g2)*x^2+(-g3*g[1]+g1)*x-g3*g[0]+g0
# y = quadgeneral is the quad through [x1,g(x1)], [x2,g(x2)], [x3,g(x3)]
# where g(x) = g3*x^3+g2*x^2+g1*x+g0 and g[0]=-x1*x2*x3, 
# g[1]=x1*x2 + x2*x3 + x3*x1, g[2] := -(x1+x2+x3).
check := [expand(subs(g[0]=-x1*x2*x3, g[1]=x1*x2+x2*x3+x3*x1, g[2]=-(x1+x2+x3), 
                 x=x1,quadgeneral)),
          expand(subs(g[0]=-x1*x2*x3, g[1]=x1*x2+x2*x3+x3*x1, g[2]=-(x1+x2+x3), 
                 x=x2,quadgeneral)),
          expand(subs(g[0]=-x1*x2*x3, g[1]=x1*x2+x2*x3+x3*x1, g[2]=-(x1+x2+x3), 
                 x=x3,quadgeneral))];
# Gives: check := [g3*x1^3+g2*x1^2+g1*x1+g0, g3*x2^3+g2*x2^2+g1*x2+g0]
# In the following, we want the defining quadratic which gives the
# y-coordinates of the linearly equiv divisor - N.B. we negate quadgeneral.
newquadx[i] := -simplify(subs( g[0] = coeff(cubx[i],x,0)/coeff(cubx[i],x,3),
                     g[1] = coeff(cubx[i],x,1)/coeff(cubx[i],x,3),
                     g[2] = coeff(cubx[i],x,2)/coeff(cubx[i],x,3),
                    g3 = coeff(cubicx[i],x,3), g2 = coeff(cubicx[i],x,2),
                    g1 = coeff(cubicx[i],x,1), g0 = coeff(cubicx[i],x,0),
                    -quadgeneral ));
# At this point, newquadx[i] should be over Q(sqrtd). We just need to remove the
# factors not involving x:
# In summary, so far, we see that the degree 3 divisor ratdiv[i] 
# is linearly equivalent to a new  degree 3 divisor
# [ [x1,y1], [x2,y2], [x3,y3] ], where x1,x2,x3 are the roots
# of cubx[i] and each yi = subs( x=xi, newquadx[i] ).
# We just need to let aa = a1 + a2*sqrtd and see what are the
# conditions for everything to be defined over Q.
# This is the same as the ratios of the coefficients of cubx[i] being in Q
# and the actual coefficients of newquadx[i] being in Q.

sqrtd := sqrt(newarr[i][4]):
eqn1 := expand(numer( radsimp( subs( aa = a1 + a2*sqrtd , 
              coeff( cubx[i],x,2 )/ coeff(cubx[i],x,3) ), `ratdenom`))):
eqn1arr := [op(eqn1)];
rateqn1 := 0:
for j from 1 to nops(eqn1arr) do if type(eqn1arr[j],ratpoly(rational,{a1,a2}))
then rateqn1 := rateqn1 + eqn1arr[j] fi od:
eqn1 := factor(simplify((eqn1 - rateqn1)/sqrtd));

eqn2 := expand(numer( radsimp( subs( aa = a1 + a2*sqrtd , 
              coeff( cubx[i],x,1 )/ coeff(cubx[i],x,3) ), `ratdenom`))):
eqn2arr := [op(eqn2)];
rateqn2 := 0:
for j from 1 to nops(eqn2arr) do if type(eqn2arr[j],ratpoly(rational,{a1,a2}))
then rateqn2 := rateqn2 + eqn2arr[j] fi od:
eqn2 := factor(simplify((eqn2 - rateqn2)/sqrtd));

eqn3 := expand(numer( radsimp( subs( aa = a1 + a2*sqrtd , 
              coeff( cubx[i],x,0 )/ coeff(cubx[i],x,3) ), `ratdenom`))):
eqn3arr := [op(eqn3)];
rateqn3 := 0:
for j from 1 to nops(eqn3arr) do if type(eqn3arr[j],ratpoly(rational,{a1,a2}))
then rateqn3 := rateqn3 + eqn3arr[j] fi od:
eqn3 := factor(simplify((eqn3 - rateqn3)/sqrtd));

eqn4 := expand(numer( radsimp( subs( aa = a1 + a2*sqrtd ,
               coeff(newquadx[i],x,2) ), `ratdenom`))):
eqn4arr := [op(eqn4)];
rateqn4 := 0:
for j from 1 to nops(eqn4arr) do if type(eqn4arr[j],ratpoly(rational,{a1,a2}))
then rateqn4 := rateqn4 + eqn4arr[j] fi od:
eqn4 := factor(simplify((eqn4 - rateqn4)/sqrtd));

eqn5 := expand(numer( radsimp( subs( aa = a1 + a2*sqrtd ,
               coeff(newquadx[i],x,1) ), `ratdenom`))):
eqn5arr := [op(eqn5)];
rateqn5 := 0:
for j from 1 to nops(eqn5arr) do if type(eqn5arr[j],ratpoly(rational,{a1,a2}))
then rateqn5 := rateqn5 + eqn5arr[j] fi od:
eqn5 := factor(simplify((eqn5 - rateqn5)/sqrtd));

eqn6 := expand(numer( radsimp( subs( aa = a1 + a2*sqrtd ,
               coeff(newquadx[i],x,0) ), `ratdenom`))):
eqn6arr := [op(eqn6)];
rateqn6 := 0:
for j from 1 to nops(eqn6arr) do if type(eqn6arr[j],ratpoly(rational,{a1,a2}))
then rateqn6 := rateqn6 + eqn6arr[j] fi od:
eqn6 := factor(simplify((eqn6 - rateqn6)/sqrtd));

eqn := gcd( gcd( gcd( gcd( gcd( eqn1, eqn2 ), eqn3 ), eqn4 ), eqn5 ), eqn6 );

projeqn := numer( subs(a1 = b1/b0, a2 = b2/b0, eqn) ):
solnset := []:
N := 30:
for j0 from -N to N do for j1 from -N to N do for j2 from -N to N do
  if igcd(j0,j1,j2)= 1 and subs(b0=j0,b1=j1,b2=j2,projeqn) = 0
  then print(j1/j0, j2/j0): solnset := [op(solnset), [j0,j1,j2]] fi 
od od od;

solnset;

if solnset <> [] then
  minsofar := max( abs(solnset[1][1]),abs(solnset[1][2]),abs(solnset[1][3])): 
  minsolnsofar := solnset[1]:
  for jj from 1 to nops(solnset) do
    if max( abs(solnset[jj][1]), abs(solnset[jj][2]), abs(solnset[jj][3]) )
           < minsofar
    then minsofar := 
         max( abs(solnset[jj][1]), abs(solnset[jj][2]), abs(solnset[jj][3]) ): 
         minsolnsofar := solnset[jj]:
    fi:
  od:
  minsolnsofar;
  aaa := minsolnsofar[2]/minsolnsofar[1]+minsolnsofar[3]/minsolnsofar[1]*sqrtd:

  finalcubicx[i] := factor( subs( aa = aaa, cubicx[i] ) ) ;
  finalcubx[i] := factor( subs( aa = aaa, cubx[i] ) ) ;
  finalnewquadx[i] := expand( factor( subs( aa = aaa, newquadx[i] ) ) );
fi:

ratdiv[i];
cubicx[i];
cubx[i];
eqn1; eqn2; eqn3;
newquadx[i];
eqn4; eqn5; eqn6;
eqn;
aaa;

# A summary of the above for Example 1 is as follows: 
# Recall, the curve is y^2 = -x^6-x^5-2*x^4-2*x^3+x^2-2*x+2
# and our Q-rational divisor class is:
# [ P1, P2, P3 ], where:
# P1 := [[1/2+1/5*sqrtm55, 7/40+12/25*sqrtm55];
# P2 := [1+sqrtm5, 2*sqrtm55];
# P3 := [1-sqrtm5, 2*sqrtm55];
# where sqrtm5 := sqrt(-5) and sqrtm55 := sqrt(-55).
# The general y = cubic in x, through P1,P2,P3 is:
# y = G(x) = a*x^3 - 5/2*x^2*a + 3/2*x^2 - 1/5*x^2*a*sqrtm55 
#            - 2/5*x^2*sqrtm55 + 7*x*a - 3*x + 2/5*x*a*sqrtm55
#            + 4/5*x*sqrtm55 - 3*a + 9 - 6/5*a*sqrtm55 - 2/5*sqrtm55;
# Replacing y by the right hand cubic in the curve, and
# removing the factor (x-x(P1))*(x-x(P2))*(x-x(P3)) leaves
# the residual cubic:
# C(x) := 10*x^3 + 35*x^2 + 2*x^2*sqrtd - 4*x*sqrtd - 50*x + 30+ 12*sqrtd
#         + 30*x^2*a - 8*x^2*a*sqrtd - 60*x*a + 16*x*a*sqrtd
#         - 8*a*sqrtd + 180*a + 10*a^2*x^3 - 2*x^2*a^2*sqrtd
#         - 25*x^2*a^2 + 70*x*a^2 + 4*x*a^2*sqrtd - 30*a^2 - 12*a^2*sqrtd;
# Let [[x1,y1],[x2,y2],[x3,y3]] be the linearly equivalent divisor.
# Then x1,x2,x3 are the roots of C(x). Furthermore, the yi = -G(xi),
# where G(x) is as above. We now find the quadratic y = H(x) which
# passes through [x1, -G(x1)], [x2, -G(x2)], [x3, -G(x3)] and find
# that the resulting H(x) is:
# H(x) := 1/10*(90+4*x^2*a^2*sqrtd - 30*x + 15*x^2 - 60*a + 8*x*sqrtd
#         - 4*x^2*sqrtd - 4*sqrtd + 4*a^2*sqrtd - 24*a*sqrtd
#         - 90*a^2 - 8*x*a^2*sqrtd - 15*x^2*a^2 + 30*x*a^2
#         - 4*x^2*a*sqrtd + 8*x*a*sqrtd - 60*x^2*a + 120*x*a)/(1+a^2);
# Now, our requirement for [[x1,y1],[x2,y2],[x3,y3]] to be Q-rational
# is the same as requiring that the ratios of the coefficients of C(x)
# are Q-rational and that the actual coefficients of H(x) are also. 
# Let a = a1 + a2*sqrtd, and let eqn1, eqn2, eqn3 be the equations
# obtained by setting the sqrtd parts to zero of:
# coeff-x^2/coeff-x^3, coeff-x/coeff-x^3, constant-term/coeff-x^3,
# respectively of C(x). Furthermore, let eqn4, eqn5, eqn6 be the
# equations obtained by setting the sqrtd parts to zero of:
# coeff-x^2, coeff-x, constant-term,respectively, of H(x).
# These are:
# -2*(a1^2+1+15*a2+55*a2^2)*(a1^2+4*a1-1+55*a2^2) = 0
# 2*(a1^2+1+15*a2+55*a2^2)*(a1^2+4*a1-1+55*a2^2) = 0
# -2*(a1^2+1+15*a2+55*a2^2)*(3*a1^2+2*a1-3+165*a2^2) = 0
# 4*(a1^2-a1-1+55*a2^2)*(a1^2+1+15*a2+55*a2^2) = 0
# -4*(a1^2-a1-1+55*a2^2)*(a1^2+1+15*a2+55*a2^2) = 0
# 2*(a1^2+1+15*a2+55*a2^2)*(a1^2-6*a1-1+55*a2^2) = 0
# which have a common factor:
# a1^2 + 1 + 15*a2 + 55*a2^2 = 0
# which is the genus 0 curve we are looking for.
# The smallest-height solution is: a1 = 1/7, a2 = -1/7,
# corresponding to a = 1/7 - 1/7*sqrtd. Substituting this
# into C(x) and H(x) gives:
# C(x) = -5/49*(5 + 2*sqrtd)*(2*x^3 + x^2 + 2*x + 2);
# H(x) = -1/2*x^2 + x - 1;
# Of course, we only need C(x) up to constants, so we can take:
# C(x) = 2*x^3 + x^2 + 2*x + 2;
# These C(x), H(x) define our Q-rational divisor class
# [[x1,y1],[x2,y2],[x3,y3]], where x1,x2,x3 are the roots
# of C(x) and yi = H(xi) for i = 1,2,3.
#   I have done all of this in Maple in a style amenable
# for mass production. That is, if a divisor curve and a divisor
# like [P1,P2,P3] is inputted, it returns C(x) and H(x).
# It should also be reasonably easy to use my procedures
# to obtain a case where the genus 0 curve has no Q-rational points
# (of course, I first should choose a curve which has no points
# in some Qp).

######################################################################

# PART 2. Maple script relevant to Example 2.

interface(prettyprint, prettyprint=false);

# We are supposing the we have the curve defined over Q:
# y^2 = k*F1*F2 = k*F1*F1'
#     = k*( a3*x^3 + a2*x^2 + a1*x + a0 )*( b3*x^3 + b2*x^2 + b1*x + b0 );
F1 := a3*x^3 + a2*x^2 + a1*x + a0;
F2 := b3*x^3 + b2*x^2 + b1*x + b0;

# where k is in Q, and F1,F2 are defined over Q(sqrt(d)) and are conjugate; so:
a3 := g3 + h3*sqrt(d); b3 := g3 - h3*sqrt(d);
a2 := g2 + h2*sqrt(d); b2 := g2 - h2*sqrt(d);
a1 := g1 + h1*sqrt(d); b1 := g1 - h1*sqrt(d);
a0 := g0 + h0*sqrt(d); b0 := g0 - h0*sqrt(d);
# where the gi,hi are in Q.
# Then [ [e1,0] + [e2,0] + [e3,0] ] is a Q-rat divisor class
# of degree 3 and the general cubic through it is:
# y = t*F1 = 
cubicx := expand( t*( a3*x^3 + a2*x^2 + a1*x + a0 ) );
# where:
t := t1 + t2*sqrt(d);
# Let s be the conjugate of t:
s :=  t1 - t2*sqrt(d);
cubicx := expand( cubicx );
# and t1,t2 are in Q.
# We want the residual cubic to have coefficient-ratios in Q.
# Substituting t*F1 for y in y^2 - k*F1*F2 and removing the known
# factor F1 gives: 

rescubx := expand( t^2*F1 - k*F2 );
conjrescubx := expand( s^2*F2 - k*F1 );
# The coefficient-ratios of rescubx are in Q iff the
# product (any coeff of rescubx)*(any coeff of conjrescubx) is in Q.

eqn1 := expand( coeff(rescubx , x^3)*coeff(conjrescubx , x^2) );
# This could also have been obtained by:
# eqn1 := ( t^2*(g3 + h3*sqrt(d)) - k*(g3 - h3*sqrt(d)) )*
#        ( s^2*(g2 - h2*sqrt(d)) - k*(g2 + h2*sqrt(d)) );
eqn1 := factor( expand( coeff(eqn1, sqrt(d)) + d*coeff(eqn1, d^(3/2))
               + d^2*coeff(eqn1, d^(5/2)) ) );

eqn2 := expand( coeff(rescubx , x^3)*coeff(conjrescubx , x) );
eqn2 := factor( expand( coeff(eqn2, sqrt(d)) + d*coeff(eqn2, d^(3/2))
               + d^2*coeff(eqn2, d^(5/2)) ) );

eqn3 := expand( coeff(rescubx , x^3)*coeff(conjrescubx , x, 0) );
eqn3 := factor( expand( coeff(eqn3, sqrt(d)) + d*coeff(eqn3, d^(3/2))
               + d^2*coeff(eqn3, d^(5/2)) ) );

eqn4 := expand( coeff(rescubx , x^2)*coeff(conjrescubx , x) );
eqn4 := factor( expand( coeff(eqn4, sqrt(d)) + d*coeff(eqn4, d^(3/2))
               + d^2*coeff(eqn4, d^(5/2)) ) );

eqn5 := expand( coeff(rescubx , x^2)*coeff(conjrescubx , x, 0) );
eqn5 := factor( expand( coeff(eqn5, sqrt(d)) + d*coeff(eqn5, d^(3/2))
               + d^2*coeff(eqn5, d^(5/2)) ) );

eqn6 := expand( coeff(rescubx , x)*coeff(conjrescubx , x, 0) );
eqn6 := factor( expand( coeff(eqn6, sqrt(d)) + d*coeff(eqn6, d^(3/2))
               + d^2*coeff(eqn6, d^(5/2)) ) );

quadgeneral := (-G3*G[2]+G2)*x^2+(-G3*G[1]+G1)*x-G3*G[0]+G0;
# y = quadgeneral is the quad through [x1,G(x1)], [x2,G(x2)], [x3,G(x3)]
# where G(x) = G3*x^3+G2*x^2+G1*x+G0 and G[0]=-x1*x2*x3, 
# G[1]=x1*x2 + x2*x3 + x3*x1, G[2] := -(x1+x2+x3).

check := [expand(subs(G[0]=-x1*x2*x3, G[1]=x1*x2+x2*x3+x3*x1, G[2]=-(x1+x2+x3), 
                 x=x1,quadgeneral)),
          expand(subs(G[0]=-x1*x2*x3, G[1]=x1*x2+x2*x3+x3*x1, G[2]=-(x1+x2+x3), 
                 x=x2,quadgeneral)),
          expand(subs(G[0]=-x1*x2*x3, G[1]=x1*x2+x2*x3+x3*x1, G[2]=-(x1+x2+x3), 
                 x=x3,quadgeneral))];
# Gives: check := [G3*x1^3+G2*x1^2+G1*x1+G0, G3*x2^3+G2*x2^2+G1*x2+G0]
# In the following, we want the defining quadratic which gives the
# y-coordinates of the linearly equiv divisor - N.B. we negate quadgeneral.
cubx := rescubx;
newquadx := -simplify(subs( G[0] = coeff(cubx,x,0)/coeff(cubx,x,3),
                     G[1] = coeff(cubx,x,1)/coeff(cubx,x,3),
                     G[2] = coeff(cubx,x,2)/coeff(cubx,x,3),
                    G3 = coeff(cubicx,x,3), G2 = coeff(cubicx,x,2),
                    G1 = coeff(cubicx,x,1), G0 = coeff(cubicx,x,0),
                    -quadgeneral ));

newquadx := expand(%);
newquadx := expand( numer( radsimp( newquadx, `ratdenom` ) ) );

eqn7 := expand( coeff(newquadx , x^2) );
eqn7 := factor( expand( coeff(eqn7, sqrt(d)) + d*coeff(eqn7, d^(3/2))
               + d^2*coeff(eqn7, d^(5/2)) ) );

eqn8 := expand( coeff(newquadx , x) );
eqn8 := factor( expand( coeff(eqn8, sqrt(d)) + d*coeff(eqn8, d^(3/2))
               + d^2*coeff(eqn8, d^(5/2)) ) );

eqn9 := expand( coeff(newquadx , x, 0) );
eqn9 := factor( expand( coeff(eqn9, sqrt(d)) + d*coeff(eqn9, d^(3/2))
               + d^2*coeff(eqn9, d^(5/2)) ) );

#   The above is summarised as follows, and includes an example in genus 2 of
# a rational divisor class which does not contain a rational divisor.
# We are supposing the we have the curve defined over Q:
# y^2 = k*F1*F2 = k*F1*F1'
#     = k*( a3*x^3 + a2*x^2 + a1*x + a0 )*( b3*x^3 + b2*x^2 + b1*x + b0 );
# F1 := a3*x^3 + a2*x^2 + a1*x + a0;
# F2 := b3*x^3 + b2*x^2 + b1*x + b0;
# where k is in Q, and F1,F2 are defined over Q(sqrt(d)) and are conjugate; so:
# a3 := g3 + h3*sqrt(d); b3 := g3 - h3*sqrt(d);
# a2 := g2 + h2*sqrt(d); b2 := g2 - h2*sqrt(d);
# a1 := g1 + h1*sqrt(d); b1 := g1 - h1*sqrt(d);
# a0 := g0 + h0*sqrt(d); b0 := g0 - h0*sqrt(d);
# where the gi,hi are in Q.
# Then [ [e1,0] + [e2,0] + [e3,0] ] is a Q-rat divisor class
# of degree 3 and the general cubic through it is:
# y = t*F1 = 
#     G(x) := expand( t*( a3*x^3 + a2*x^2 + a1*x + a0 ) );
# where:
# t := t1 + t2*sqrt(d);
# Let s be the conjugate of t:
# s :=  t1 - t2*sqrt(d);
# and t1,t2 are in Q.
# We want the residual cubic to have coefficient-ratios in Q.
# Substituting t*F1 for y in y^2 - k*F1*F2 and removing the known
# factor F1 gives the residual cubic: 
# rescubx := expand( t^2*F1 - k*F2 );
# So, if [ [x1,y1], [x2,y2], [x3,y3] ] is our linearly equiv divisor,
# then x1,x2,x3 are the roots of rescubx.
# conjrescubx := expand( s^2*F2 - k*F1 );
# The coefficient-ratios of rescubx are in Q iff the
# product (any coeff of rescubx)*(any coeff of conjrescubx) is in Q.
# This gives eqn1 , ... , eqn6 as:
# -(k+t2^2*d-t1^2)*(-k+t2^2*d-t1^2)*(g3*h2-g2*h3)
# -(k+t2^2*d-t1^2)*(-k+t2^2*d-t1^2)*(g3*h1-g1*h3)
# -(k+t2^2*d-t1^2)*(-k+t2^2*d-t1^2)*(g3*h0-g0*h3)
# -(k+t2^2*d-t1^2)*(-k+t2^2*d-t1^2)*(-g1*h2+g2*h1)
# -(k+t2^2*d-t1^2)*(-k+t2^2*d-t1^2)*(-g0*h2+g2*h0)
# -(k+t2^2*d-t1^2)*(-k+t2^2*d-t1^2)*(g1*h0-g0*h1)
#
# Furthermore, if we compute the y = quadratic-in-x through
# the points [x1, G(x1)], [x2, G(x2)], [x3, G(x3)], and set the
# coeff of sqrt(d) to 0 we get eqn7, eqn8, eqn9 as:
#
# 2*k*(g3*h2-g2*h3)*(h3*d*t2+g3*t1)*(k+t2^2*d-t1^2)
# 2*k*(g3*h1-g1*h3)*(h3*d*t2+g3*t1)*(k+t2^2*d-t1^2)
# 2*k*(g3*h0-g0*h3)*(h3*d*t2+g3*t1)*(k+t2^2*d-t1^2)
#
# Now, looking at the first 6 eqns, we see that we cannot have all of:
# g3*h2-g2*h3, g3*h1-g1*h3, g3*h0-g0*h3, -g1*h2+g2*h1, -g0*h2+g2*h0, g1*h0-g0*h1,
# since then the original curve would have zero discriminant.
# So, either N(t) = k or N(t) = -k.
#   Imagine that N(t) = -k. Then from the last 3 equations, either
#   g3*h2-g2*h3 = g3*h1-g1*h3 = g3*h0-g0*h3 = 0, which would again
#   make the discrim of the original curve 0 (both when g3=h4=0 and otherwise).
#   Then we would be forced into h3*d*t2+g3*t1 = 0 as well as N(t) = -k;
#   these force newquadx = 0 and rescubic to be:
#   2*t2^2*d*(d^(1/2)*g3-h3*d)*(g3*x^2*h2+g3*x*h1+g3*h0-x^2*h3*g2-g0*h3-x*h3*g1)
#   Now, the coeff-of-x^3 of rescubic is 0, and if the remaining part
#   g3*x^2*h2+g3*x*h1+g3*h0-x^2*h3*g2-g0*h3-x*h3*g1 = 0, then we would
#   again have zero discrim of the original curve. Otherwise, we must
#   have points at infinity occurring, but clearly f6 = d mod squares,
#   and our intersection cannot have both branches. We seem to
#   be forces to concluce that d is a square, when of course there
#   is indeed a rational divisor.
# Conclusion: 
# There exists a rational divisor linearly equivalent to our starting
# rat divisor class when either:
# d is a square (in which case our starting divisor is already Q-rational)
# or:
# N(t) = k for some t, in which case the lin equiv divisor is cut out by:
# t*F1 - t'*F1'.
# 4*d*N(t)^2*(t2*g3+t1*h3)*( (g3*h2-g2*h3)*x^2 + (g3*h1-g1*h3)*x + g3*h0-g0*h3 )
#
# So, just need to test these in each Qp.
#
# Ex A. k = 7, d = 2, F1 = any cubic in x over Q(sqrt(2)) [and not over Q]
#                          with nonzero discriminant,
#       and let C : y^2 = k*F1*F1'.  For t = t1 + t2*sqrt(d) the genus 0 
#       equation in t1,t2 is: t1^2 - d*t2^2 = k. Take t1=3, t2=1, t = 3 + sqrt(2).
#       Let e1,e2,e3 be the roots of F1. Then the Q-rational divisor
#       class [ [e1,0], [e2,0], [e3,0] ] is birationally equivalent
#       (via the function y - t*F1) to:
#       [ [x1,y1], [x2,y2], [x3,y3] ], where x1,x2,x3 are the roots
#       of t*F1 - t'*F1', and yi = Q(xi), where Q(x) is:
#       4*d*N(t)^2*(t2*g3+t1*h3)*((g3*h2-g2*h3)*x^2+(g3*h1-g1*h3)*x+g3*h0-g0*h3).
#
# Ex B. k = 5, d = 2,  F1 = any cubic in x over Q(sqrt(2)) [and not over Q] 
#                          with nonzero discriminant, 
#       and let C : y^2 = k*F1*F1'.  For t = t1 + t2*sqrt(d) the genus 0 
#       equation in t1,t2 is: t1^2 - d*t2^2 = k. Then in Q_5 we see
#       that there are no solutions to t1^2 - d*t2^2 = k, and furthermore
#       d is not a square in Q_5. Therefore, there is no rational
#       divisor over Q_5 in our Q-rational divisor class, and so
#       also no Q-rational divisor.
#
# Note that, at least for these examples, your instinct is correct
# about primes of good reduction never giving a contradiction.
# The only primes which might not have a Q_p soln to t1^2 - d*t2^2 = k
# are primes dividing d or k, and all of these are of bad reduction.