#!/bin/sh
# this archive contains:
# -rw-------   1 bruin    staff      14341 Dec 27 19:07 divcalc.mpl
# -rw-------   1 bruin    staff      60125 Dec 27 19:07 flynn.mpl
# It is produced by a very simple shar that does not protect special
# symbols, but does minimise work for manual extraction.
# If possible, treat this data as binary.
#
#--- flynn.mpl (start) ---#
cat << "#--- flynn.mpl (end) ---#" > flynn.mpl # UNIX extraction command
##############################################################################
##
## flynn.mpl
##
## version 1.0, by Nils Bruin
## date: march 26, 1997
##
## Offers an interface to various routines that can be used for calculations
## on genus-2 curves.
## special requirements are: f0 ... f6 should be set to the coefficients of
## the hyperelliptic model of the curve prior to loading this file.
## the routines defined are then valid for the described curves.
## s1 ... s15, a0 ... a15 and b0 ... b15 should remain undefined.
##
## Points on the jacobian have different representations:
## div :  divisor form. [[x,y],[u,v]]. In principle, any point can be
##        represented this way (using infinity). See divcalc.mpl for that.
##        Most routines in this file won't work, however.
## coord: projective coordinates on the P_15 model of J.
## loccoord: local coordinates around the neutral element. (a1/a0,a2/a0) if
##           [a0...a15] is the coord-representation.
## kummer: [a14,a13,a12,a5] from coord.
## theta: a13^2-4*a13*a12 (theta=0 for points on diagonal)
##
## most relevant conversion-routines (denoted by <type1>2<type2>) are available
## but might have a limited scope of validity. Converions from loccoord have
## automatically a limited precision.
##
## Apart from that, the following routines are available:
## loclog(S)=L  :  S is loccoord, L is approx. to logarithm of S
## locexp(L)=S  : approximation to inverse of loclog
## sumkummer(index,C1,C2) :  returns kummer-coordinates of two points in
## coord-representation. index chooses from 4 different versions of the
## function (some may not be defined ie. [0,0,0,0])
## The expressions in div2coord, div2loccoord and in bilinform are obtained
## from the files by E.V. Flynn, available from
## ftp://ftp.liv.ac.uk/~ftp/pub/genus2
## as are the power-series in loclog, locexp and localexpansions.
##
##############################################################################

##############################################################################
##
## div2coord
##
## converts a divisor-representation of a point on the jacobian to coordinates
## in P_15. Routine is for generic points. Special points might need some
## limits. The image returned of the generic point [[x,y],[u,v]] is correct.
## and can be used to calculate all kinds of limits.
##############################################################################

div2coord:=proc(div)
  local f0xu,f1xu,gxu,gux,hxu,hux,x,y,u,v,
       a0,a1,a2,a3,a4,a5,a6,a7,a8,a9,a10,a11,a12,a13,a14,a15;

  x:=div[1][1];y:=div[1][2];u:=div[2][1];v:=div[2][2];

  f0xu := 2*f0+f1*(x+u)+2*f2*(x*u)+f3*(x+u)*(x*u)
              +2*f4*(x*u)**2+f5*(x+u)*(x*u)**2+2*f6*(x*u)**3;
  f1xu := f0*(x+u)+2*f1*(x*u)+f2*(x+u)*(x*u)+2*f3*(x*u)**2
              +f4*(x+u)*(x*u)**2+2*f5*(x*u)**3+f6*(x+u)*(x*u)**3;
  gxu := f0*4+f1*(x+3*u)+f2*(2*x*u+2*u**2)+f3*(3*x*u**2+u**3)
          +f4*(4*x*u**3)+f5*x*(x*u**3+3*u**4)+f6*2*x*(x*u**4+u**5);
  gux := f0*4+f1*(u+3*x)+f2*(2*u*x+2*x**2)+f3*(3*u*x**2+x**3)
          +f4*(4*u*x**3)+f5*u*(u*x**3+3*x**4)+f6*2*u*(u*x**4+x**5);
  hxu := f0*2*(x+u)+f1*u*(3*x+u)+f2*4*x*u**2+f3*x*u**2*(x+3*u)
	 +f4*2*x*u**3*(x+u)+f5*x*u**4*(3*x+u)+f6*4*x**2*u**5;
  hux := f0*2*(u+x)+f1*x*(3*u+x)+f2*4*u*x**2+f3*u*x**2*(u+3*x)
          +f4*2*u*x**3*(u+x)+f5*u*x**4*(3*u+x)+f6*4*u**2*x**5;
  a15 := (x-u)**2; a14 := 1; a13 := x + u; a12 := x*u; 
  a11 := x*u*(x+u); a10 := (x*u)**2; a9 := (y-v)/(x-u);
  a8 := (u*y-x*v)/(x-u); a7 := (u**2*y-x**2*v)/(x-u);
  a6 := (u**3*y-x**3*v)/(x-u); a5 := (f0xu-2*y*v)/((x-u)**2);
  a4 := (f1xu-(x+u)*y*v)/((x-u)**2); a3 := (x*u)*a5;
  a2:=(gxu*y-gux*v)/((x-u)**3); a1:=(hxu*y-hux*v)/((x-u)**3);
  a0 :=a5**2;

  RETURN([a0,a1,a2,a3,a4,a5,a6,a7,a8,a9,a10,a11,a12,a13,a14,a15]);
end;

##############################################################################
##
## projection of coord to kummer surface
##
##############################################################################

coord2kummer:=(a)->[a[14+1],a[13+1],a[12+1],a[5+1]];

##############################################################################
##
## local coordinates from coord
## (might be not defined)
##
##############################################################################

coord2loccoord:=(a)->[a[1+1]/a[0+1],a[2+1]/a[0+1]];

##############################################################################
##
## directly from div to loccoord.
## Should work in most cases.
##
##############################################################################

div2loccoord:=proc(div)
  local gxu,gux,hxu,hux,f0xu,a012defs,x,u,y,v,S1,S2;

  x:=div[1][1];y:=div[1][2];u:=div[2][1];v:=div[2][2];

  gxu := f0*4+f1*(x+3*u)+f2*(2*x*u+2*u**2)+f3*(3*x*u**2+u**3)
	+f4*(4*x*u**3)+f5*x*(x*u**3+3*u**4)+f6*2*x*(x*u**4+u**5);
  gux := f0*4+f1*(u+3*x)+f2*(2*u*x+2*x**2)+f3*(3*u*x**2+x**3)
	+f4*(4*u*x**3)+f5*u*(u*x**3+3*x**4)+f6*2*u*(u*x**4+x**5);
  hxu := f0*2*(x+u)+f1*u*(3*x+u)+f2*4*x*u**2+f3*x*u**2*(x+3*u)
	+f4*2*x*u**3*(x+u)+f5*x*u**4*(3*x+u)+f6*4*x**2*u**5;
  hux := f0*2*(u+x)+f1*x*(3*u+x)+f2*4*u*x**2+f3*u*x**2*(u+3*x)
	+f4*2*u*x**3*(u+x)+f5*u*x**4*(3*u+x)+f6*4*u**2*x**5;
  f0xu := 2*f0+f1*(x+u)+2*f2*(x*u)+f3*(x+u)*(x*u)
	 +2*f4*(x*u)**2+f5*(x+u)*(x*u)**2+2*f6*(x*u)**3;
  a012defs:={a2=(gxu*y-gux*v)/((x-u)**3),
            a1=(hxu*y-hux*v)/((x-u)**3),
            a0=(f0xu-2*y*v)**2/((x-u)**4)};
  S1:=simplify(subs(a012defs,a1/a0));
  S2:=simplify(subs(a012defs,a2/a0));
  RETURN(subs({x=div[1][1],y=div[1][2],u=div[2][1],v=div[2][2]},[S1,S2]));
end;

##############################################################################
##
## local expansions of s1..s15 in s1,s2.
##
##############################################################################

localexpansions:=

{s0 = 1,

s1 = s1,

s2 = s2,

s3 = s1^2-f0*s2^4-f4*s1^4-5*f4^3*s1^8+2*f4^2*s1^6-f2*f5^2*s1^8-5*f2^2*f0*s2^8-
f6^2*f0*s1^8-f6*f3^2*s1^8-f2*f5*f1*s1^4*s2^4-14*f2*f6*f0*s1^4*s2^4-16*f2*f6*f4*
s1^8-4*f2*f5*f0*s1^3*s2^5-4*f2*f0*f4*s1^2*s2^6+3*f6*f5*f1*s1^8+3*f6*f1^2*s1^4*
s2^4+f5*f1*s1^4*s2^2-f5*f3*s1^6+12*f6*f5*f0*s1^7*s2-8*f6*f1*f0*s1^3*s2^5-24*f6*
f1*f4*s1^7*s2+6*f6*f1*f3*s1^6*s2^2-52*f6*f0^2*s1^2*s2^6-40*f6*f0*f4*s1^6*s2^2+
24*f6*f0*f3*s1^5*s2^3-2*f5^2*f1*s1^7*s2-2*f5^2*f0*s1^6*s2^2-2*f5*f1*f0*s1^2*s2^
6-4*f5*f1*f4*s1^6*s2^2-20*f5*f0^2*s1*s2^7-12*f5*f0*f4*s1^5*s2^3+5*f5*f0*f3*s1^4
*s2^4+5*f5*f4*f3*s1^8+2*f1*f0*f3*s2^8-9*f0^2*f4*s2^8-6*f0*f4^2*s1^4*s2^4+4*f2*
f6*s1^6+2*f2*f0*s2^6+8*f6*f1*s1^5*s2+18*f6*f0*s1^4*s2^2+4*f5*f0*s1^3*s2^3+2*f0*
f4*s1^2*s2^4,

s4 = s1*s2-f2^2*f5*s1^4*s2^4+f2*f6*f5*s1^8-f6^2*f1*s1^8+f2*f6*f1*s1^4*s2^4+f5^
2*f1*s1^6*s2^2-f5*f0^2*s2^8+f5*f1^2*s1^2*s2^6+f6*f3*s1^6+f1*f4*s1^2*s2^4+f0*f3*
s2^6+f1*f0*f4*s2^8+f5*f0*f4*s1^4*s2^4-f1*f4^2*s1^4*s2^4+f2*f5*s1^4*s2^2-4*f2*f6
*f0*s1^3*s2^5-8*f2*f6*f4*s1^7*s2+2*f2*f6*f3*s1^6*s2^2-3*f2*f5*f1*s1^3*s2^5-10*
f2*f5*f0*s1^2*s2^6-2*f2*f5*f4*s1^6*s2^2-2*f2*f1*f4*s1^2*s2^6-8*f2*f0*f4*s1*s2^7
-3*f2*f0*f3*s2^8-4*f6^2*f0*s1^7*s2+5*f6*f5*f1*s1^7*s2+16*f6*f5*f0*s1^6*s2^2+6*
f6*f1^2*s1^3*s2^5+16*f6*f1*f0*s1^2*s2^6-10*f6*f1*f4*s1^6*s2^2+10*f6*f1*f3*s1^5*
s2^3-4*f6*f0^2*s1*s2^7-4*f6*f0*f4*s1^5*s2^3+30*f6*f0*f3*s1^4*s2^4-3*f6*f4*f3*s1
^8+6*f5^2*f0*s1^5*s2^3+5*f5*f1*f0*s1*s2^7-3*f5*f1*f4*s1^5*s2^3+2*f5*f1*f3*s1^4*
s2^4+10*f5*f0*f3*s1^3*s2^5+2*f0*f4*f3*s1^2*s2^6+4*f2*f6*s1^5*s2+9*f6*f1*s1^4*s2
^2+20*f6*f0*s1^3*s2^3+3*f5*f1*s1^3*s2^3+9*f5*f0*s1^2*s2^4+4*f0*f4*s1*s2^5,

s5 = s2^2-f2*s2^4-f6*s1^4-5*f2^3*s2^8+2*f2^2*s2^6-f6*f0^2*s2^8-f0*f3^2*s2^8-f1
^2*f4*s2^8-f5*f1*f4*s1^4*s2^4-6*f2^2*f6*s1^4*s2^4-9*f2*f6^2*s1^8-12*f2*f6*f1*s1
^3*s2^5+f5*f1*s1^2*s2^4-f1*f3*s2^6-40*f2*f6*f0*s1^2*s2^6-4*f2*f6*f4*s1^6*s2^2-4
*f2*f5*f1*s1^2*s2^6-24*f2*f5*f0*s1*s2^7+5*f2*f1*f3*s2^8-16*f2*f0*f4*s2^8-20*f6^
2*f1*s1^7*s2-52*f6^2*f0*s1^6*s2^2-2*f6*f5*f1*s1^6*s2^2-8*f6*f5*f0*s1^5*s2^3+2*
f6*f5*f3*s1^8-2*f6*f1^2*s1^2*s2^6+12*f6*f1*f0*s1*s2^7-4*f6*f1*f4*s1^5*s2^3+5*f6
*f1*f3*s1^4*s2^4-14*f6*f0*f4*s1^4*s2^4+24*f6*f0*f3*s1^3*s2^5-5*f6*f4^2*s1^8+3*
f5^2*f0*s1^4*s2^4-2*f5*f1^2*s1*s2^7+3*f5*f1*f0*s2^8+6*f5*f0*f3*s1^2*s2^6+2*f2*
f6*s1^4*s2^2+4*f6*f1*s1^3*s2^3+18*f6*f0*s1^2*s2^4+2*f6*f4*s1^6+8*f5*f0*s1*s2^5+
4*f0*f4*s2^6,

s6 = s1^3-f4*s1^5-5*f4^3*s1^9+2*f4^2*s1^7-f2*f5^2*s1^9-f6*f3^2*s1^9-f6^2*f0*s1
^9-f2*f1*f4*s1^4*s2^5+f0^2*f3*s2^9-f5*f3*s1^7-f1*f0*f3*s1*s2^8+f3*s1^4*s2+4*f2^
3*s1^3*s2^6+20*f2^2*f6*s1^7*s2^2+6*f2^2*f1*s1^2*s2^7+9*f2^2*f0*s1*s2^8+2*f2^2*
f4*s1^5*s2^4+70*f2*f6*f1*s1^6*s2^3+120*f2*f6*f0*s1^5*s2^4-22*f2*f6*f4*s1^9+12*
f2*f6*f3*s1^8*s2+7*f2*f5*f1*s1^5*s2^4+f1^2*s1*s2^6+f1*f0*s2^7+f1*f4*s1^4*s2^3+
24*f2*f5*f0*s1^4*s2^5-6*f2*f5*f4*s1^8*s2-2*f2*f5*f3*s1^7*s2^2-3*f2*f1^2*s1*s2^8
-3*f2*f1*f0*s2^9-8*f2*f1*f3*s1^3*s2^6-12*f2*f0*f3*s1^2*s2^7+4*f2*f4^2*s1^7*s2^2
+5*f6*f5*f1*s1^9+22*f6*f5*f0*s1^8*s2+66*f6*f1^2*s1^5*s2^4+231*f6*f1*f0*s1^4*s2^
5-35*f6*f1*f4*s1^8*s2+34*f6*f1*f3*s1^7*s2^2+148*f6*f0^2*s1^3*s2^6-56*f6*f0*f4*
s1^7*s2^2+96*f6*f0*f3*s1^6*s2^3+8*f5^2*f0*s1^7*s2^2+19*f5*f1^2*s1^4*s2^5+90*f5*
f1*f0*s1^3*s2^6-16*f5*f1*f4*s1^7*s2^2+6*f5*f1*f3*s1^6*s2^3+64*f5*f0^2*s1^2*s2^7
-24*f5*f0*f4*s1^6*s2^3+39*f5*f0*f3*s1^5*s2^4+5*f5*f4*f3*s1^9-2*f5*f3^2*s1^8*s2+
4*f1^2*f4*s1^3*s2^6-3*f1^2*f3*s1^2*s2^7+38*f1*f0*f4*s1^2*s2^7-2*f1*f4^2*s1^6*s2
^3-3*f1*f4*f3*s1^5*s2^4+31*f0^2*f4*s1*s2^8-14*f0*f4^2*s1^5*s2^4+10*f0*f4*f3*s1^
4*s2^5+4*f0*f3^2*s1^3*s2^6+5*f4^2*f3*s1^8*s2-2*f2^2*s1^3*s2^4+6*f2*f6*s1^7+2*f2
*f5*s1^6*s2-3*f2*f1*s1^2*s2^5-4*f2*f0*s1*s2^6-2*f2*f4*s1^5*s2^2+13*f6*f1*s1^6*
s2+30*f6*f0*s1^5*s2^2+7*f5*f1*s1^5*s2^2+18*f5*f0*s1^4*s2^3+3*f1*f3*s1^3*s2^4+10
*f0*f4*s1^3*s2^4+6*f0*f3*s1^2*s2^5-2*f4*f3*s1^6*s2+2*f2*s1^3*s2^2+3*f1*s1^2*s2^
3+3*f0*s1*s2^4,

s7 = s1^2*s2+f0*s2^5-f6^2*f1*s1^9-f2*f5^2*s1^8*s2-f6*f3^2*s1^8*s2-f5*f3*s1^6*
s2+5*f2^2*f0*s2^9+f1*s1*s2^4-f4*s1^4*s2+5*f2^2*f1*s1*s2^8+8*f2*f6*f1*s1^5*s2^4+
6*f2*f6*f0*s1^4*s2^5-16*f2*f6*f4*s1^8*s2-2*f2*f5*f1*s1^4*s2^5-8*f2*f5*f0*s1^3*
s2^6-12*f2*f0*f4*s1^2*s2^7-6*f2*f0*f3*s1*s2^8-7*f6^2*f0*s1^8*s2+2*f6*f5*f1*s1^8
*s2+8*f6*f5*f0*s1^7*s2^2+17*f6*f1^2*s1^4*s2^5+76*f6*f1*f0*s1^3*s2^6-24*f6*f1*f4
*s1^7*s2^2+8*f6*f1*f3*s1^6*s2^3+44*f6*f0^2*s1^2*s2^7-48*f6*f0*f4*s1^6*s2^3+34*
f6*f0*f3*s1^5*s2^4-2*f5^2*f1*s1^7*s2^2+3*f5*f1^2*s1^3*s2^6+30*f5*f1*f0*s1^2*s2^
7-6*f5*f1*f4*s1^6*s2^3+22*f5*f0^2*s1*s2^8-18*f5*f0*f4*s1^5*s2^4+11*f5*f0*f3*s1^
4*s2^5+5*f5*f4*f3*s1^8*s2-2*f1^2*f3*s1*s2^8+12*f1*f0*f4*s1*s2^8-2*f1*f0*f3*s2^9
+11*f0^2*f4*s2^9-10*f0*f4^2*s1^4*s2^5-5*f4^3*s1^8*s2+4*f2*f6*s1^6*s2-2*f2*f1*s1
*s2^6-2*f2*f0*s2^7+8*f6*f1*s1^5*s2^2+22*f6*f0*s1^4*s2^3+2*f5*f1*s1^4*s2^3+10*f5
*f0*s1^3*s2^4+6*f0*f4*s1^2*s2^5+2*f0*f3*s1*s2^6+2*f4^2*s1^6*s2,

s8 = s1*s2^2+f6*s1^5-f1^2*f4*s1*s2^8-f5*f0^2*s2^9-f1*f3*s1*s2^6-5*f2^3*s1*s2^8
-f0*f3^2*s1*s2^8-f2*s1*s2^4+f5*s1^4*s2-10*f2^2*f6*s1^5*s2^4+11*f2*f6^2*s1^9+12*
f2*f6*f5*s1^8*s2-18*f2*f6*f1*s1^4*s2^5-48*f2*f6*f0*s1^3*s2^6-12*f2*f6*f4*s1^7*
s2^2-6*f2*f5*f1*s1^3*s2^6-24*f2*f5*f0*s1^2*s2^7+5*f2*f1*f3*s1*s2^8-16*f2*f0*f4*
s1*s2^8+22*f6^2*f1*s1^8*s2+44*f6^2*f0*s1^7*s2^2+30*f6*f5*f1*s1^7*s2^2+76*f6*f5*
f0*s1^6*s2^3-2*f6*f5*f3*s1^9+8*f6*f1*f0*s1^2*s2^7-8*f6*f1*f4*s1^6*s2^3+11*f6*f1
*f3*s1^5*s2^4-7*f6*f0^2*s1*s2^8+6*f6*f0*f4*s1^5*s2^4+34*f6*f0*f3*s1^4*s2^5+5*f6
*f4^2*s1^9-6*f6*f4*f3*s1^8*s2+3*f5^2*f1*s1^6*s2^3+17*f5^2*f0*s1^5*s2^4-2*f5^2*
f3*s1^8*s2-2*f5*f1^2*s1^2*s2^7+2*f5*f1*f0*s1*s2^8-2*f5*f1*f4*s1^5*s2^4+8*f5*f0*
f4*s1^4*s2^5+8*f5*f0*f3*s1^3*s2^6+5*f5*f4^2*s1^8*s2+2*f2^2*s1*s2^6+6*f2*f6*s1^5
*s2^2+10*f6*f1*s1^4*s2^3+22*f6*f0*s1^3*s2^4-2*f6*f4*s1^7+2*f6*f3*s1^6*s2+2*f5*
f1*s1^3*s2^4+8*f5*f0*s1^2*s2^5-2*f5*f4*s1^6*s2+4*f0*f4*s1*s2^6,

s9 = -5*f2^3*s2^9+s2^3-f2*s2^5+2*f2^2*s2^7+f6^2*f3*s1^9-f6*f0^2*s2^9-f2*f5*f4*
s1^5*s2^4-f1^2*f4*s2^9-f6*f5*f3*s1^8*s2+f2*f5*s1^3*s2^4+f6*f5*s1^7+f5^2*s1^6*s2
-f0*f3^2*s2^9+f3*s1*s2^4-14*f2^2*f6*s1^4*s2^5-2*f2^2*f5*s1^3*s2^6+4*f2^2*f4*s1^
2*s2^7+5*f2^2*f3*s1*s2^8+31*f2*f6^2*s1^8*s2+38*f2*f6*f5*s1^7*s2^2-24*f2*f6*f1*
s1^3*s2^6-56*f2*f6*f0*s1^2*s2^7+10*f2*f6*f3*s1^5*s2^4-f1*f3*s2^7+4*f2*f5^2*s1^6
*s2^3-16*f2*f5*f1*s1^2*s2^7-35*f2*f5*f0*s1*s2^8-3*f2*f5*f3*s1^4*s2^5-6*f2*f1*f4
*s1*s2^8+5*f2*f1*f3*s2^9-22*f2*f0*f4*s2^9+2*f2*f4^2*s1^4*s2^5+64*f6^2*f1*s1^7*
s2^2+148*f6^2*f0*s1^6*s2^3+90*f6*f5*f1*s1^6*s2^3+231*f6*f5*f0*s1^5*s2^4-3*f6*f5
*f4*s1^9+8*f6*f1^2*s1^2*s2^7+22*f6*f1*f0*s1*s2^8+24*f6*f1*f4*s1^5*s2^4+39*f6*f1
*f3*s1^4*s2^5+120*f6*f0*f4*s1^4*s2^5+96*f6*f0*f3*s1^3*s2^6+9*f6*f4^2*s1^8*s2-12
*f6*f4*f3*s1^7*s2^2+4*f6*f3^2*s1^6*s2^3+19*f5^2*f1*s1^5*s2^4+66*f5^2*f0*s1^4*s2
^5-3*f5^2*f4*s1^8*s2-3*f5^2*f3*s1^7*s2^2+5*f5*f1*f0*s2^9+7*f5*f1*f4*s1^4*s2^5+6
*f5*f1*f3*s1^3*s2^6+70*f5*f0*f4*s1^3*s2^6+34*f5*f0*f3*s1^2*s2^7+6*f5*f4^2*s1^7*
s2^2-8*f5*f4*f3*s1^6*s2^3-2*f1*f4*f3*s1^2*s2^7-2*f1*f3^2*s1*s2^8+20*f0*f4^2*s1^
2*s2^7+12*f0*f4*f3*s1*s2^8+4*f4^3*s1^6*s2^3+10*f2*f6*s1^4*s2^3-2*f2*f4*s1^2*s2^
5-2*f2*f3*s1*s2^6+18*f6*f1*s1^3*s2^4+30*f6*f0*s1^2*s2^5-4*f6*f4*s1^6*s2+6*f6*f3
*s1^5*s2^2+7*f5*f1*s1^2*s2^5+13*f5*f0*s1*s2^6-3*f5*f4*s1^5*s2^2+3*f5*f3*s1^4*s2
^3+2*f1*f4*s1*s2^6+6*f0*f4*s2^7-2*f4^2*s1^4*s2^3+3*f6*s1^4*s2+3*f5*s1^3*s2^2+2*
f4*s1^2*s2^3,

s10 = f0^2*s2^8-14*f4^3*s1^10+5*f4^2*s1^8-2*f4*s1^6+s1^4-2*f2*f5^2*s1^10-10*f2
^2*f0*s1^2*s2^8-40*f2*f6*f4*s1^10-36*f2*f6*f0*s1^6*s2^4-2*f2*f5*f1*s1^6*s2^4-8*
f2*f5*f0*s1^5*s2^5-4*f2*f0^2*s2^10-12*f2*f0*f4*s1^4*s2^6-2*f6^2*f0*s1^10+6*f6*
f5*f1*s1^10+24*f6*f5*f0*s1^9*s2+6*f6*f1^2*s1^6*s2^4-32*f6*f1*f0*s1^5*s2^5-64*f6
*f1*f4*s1^9*s2+12*f6*f1*f3*s1^8*s2^2-140*f6*f0^2*s1^4*s2^6-116*f6*f0*f4*s1^8*s2
^2+48*f6*f0*f3*s1^7*s2^3-2*f6*f3^2*s1^10-4*f5^2*f1*s1^9*s2-4*f5^2*f0*s1^8*s2^2-
6*f5*f1*f0*s1^4*s2^6-10*f5*f1*f4*s1^8*s2^2-48*f5*f0^2*s1^3*s2^7-32*f5*f0*f4*s1^
7*s2^3+12*f5*f0*f3*s1^6*s2^4+12*f5*f4*f3*s1^10+4*f1*f0*f3*s1^2*s2^8-22*f0^2*f4*
s1^2*s2^8-20*f0*f4^2*s1^6*s2^4+8*f2*f6*s1^8+4*f2*f0*s1^2*s2^6+16*f6*f1*s1^7*s2+
36*f6*f0*s1^6*s2^2+2*f5*f1*s1^6*s2^2+8*f5*f0*s1^5*s2^3-2*f5*f3*s1^8+6*f0*f4*s1^
4*s2^4-2*f0*s1^2*s2^4,

s11 = 2*s1^3*s2+f5*s1^6-f6^2*f1*s1^10+f5^2*f1*s1^8*s2^2-f1*f3^2*s1^4*s2^6+f6*
f3*s1^8-f0^2*f3*s2^10-f2*f3*s1^4*s2^4-2*f2^2*f5*s1^6*s2^4+5*f2^2*f1*s1^2*s2^8-
10*f2^2*f0*s1*s2^9+2*f2^2*f3*s1^4*s2^6+14*f2*f6*f5*s1^10+12*f2*f6*f1*s1^6*s2^4-
44*f2*f6*f0*s1^5*s2^5-56*f2*f6*f4*s1^9*s2+14*f2*f6*f3*s1^8*s2^2-f1*f0*s2^8+f1*
f4*s1^4*s2^4+f1*s1^2*s2^4+f3*s1^4*s2^2-2*f2*f5^2*s1^9*s2-8*f2*f5*f1*s1^5*s2^5-
24*f2*f5*f0*s1^4*s2^6-6*f2*f5*f4*s1^8*s2^2+4*f2*f1*f0*s2^10-2*f2*f1*f4*s1^4*s2^
6-24*f2*f0*f4*s1^3*s2^7+2*f2*f4*f3*s1^6*s2^4-10*f6^2*f0*s1^9*s2+40*f6*f5*f1*s1^
9*s2+110*f6*f5*f0*s1^8*s2^2+34*f6*f1^2*s1^5*s2^5+54*f6*f1*f0*s1^4*s2^6-80*f6*f1
*f4*s1^8*s2^2+52*f6*f1*f3*s1^7*s2^3-152*f6*f0^2*s1^3*s2^7-128*f6*f0*f4*s1^7*s2^
3+162*f6*f0*f3*s1^6*s2^4-4*f6*f4*f3*s1^10-2*f6*f3^2*s1^9*s2+20*f5^2*f0*s1^7*s2^
3-3*f5^2*f3*s1^10+5*f5*f1^2*s1^4*s2^6+20*f5*f1*f0*s1^3*s2^7-20*f5*f1*f4*s1^7*s2
^3+6*f5*f1*f3*s1^6*s2^4-57*f5*f0^2*s1^2*s2^8-28*f5*f0*f4*s1^6*s2^4+46*f5*f0*f3*
s1^5*s2^5+9*f5*f4^2*s1^10+10*f5*f4*f3*s1^9*s2-2*f5*f3^2*s1^8*s2^2-2*f1^2*f3*s1^
2*s2^8+10*f1*f0*f4*s1^2*s2^8+4*f1*f0*f3*s1*s2^9-2*f1*f4^2*s1^6*s2^4-26*f0^2*f4*
s1*s2^9-20*f0*f4^2*s1^5*s2^5+12*f0*f4*f3*s1^4*s2^6-10*f4^3*s1^9*s2+5*f4^2*f3*s1
^8*s2^2+16*f2*f6*s1^7*s2+2*f2*f5*s1^6*s2^2-2*f2*f1*s1^2*s2^6+4*f2*f0*s1*s2^7+32
*f6*f1*s1^6*s2^2+76*f6*f0*s1^5*s2^3+8*f5*f1*s1^5*s2^3+23*f5*f0*s1^4*s2^4-3*f5*
f4*s1^8-2*f5*f3*s1^7*s2+12*f0*f4*s1^3*s2^5+4*f4^2*s1^7*s2-2*f4*f3*s1^6*s2^2-2*
f0*s1*s2^5-2*f4*s1^5*s2,

s12 = -f6*s1^6+s1^2*s2^2-f2*s1^2*s2^4-f0*s2^6-5*f2^3*s1^2*s2^8-10*f2^2*f6*s1^6
*s2^4-9*f2^2*f0*s2^10-2*f2^2*f4*s1^4*s2^6-13*f2*f6^2*s1^10-20*f2*f6*f1*s1^5*s2^
5-76*f2*f6*f0*s1^4*s2^6-22*f2*f6*f4*s1^8*s2^2-f2*f5^2*s1^8*s2^2-6*f2*f5*f1*s1^4
*s2^6+5*f2*f1*f3*s1^2*s2^8-22*f2*f0*f4*s1^2*s2^8-2*f2*f4^2*s1^6*s2^4-28*f6^2*f1
*s1^9*s2-71*f6^2*f0*s1^8*s2^2+3*f6*f5*f3*s1^10+f6*f1^2*s1^4*s2^6+11*f6*f1*f3*s1
^6*s2^4-71*f6*f0^2*s1^2*s2^8-76*f6*f0*f4*s1^6*s2^4+48*f6*f0*f3*s1^5*s2^5-9*f6*
f4^2*s1^10-f6*f3^2*s1^8*s2^2-2*f5^2*f1*s1^7*s2^3+f5^2*f0*s1^6*s2^4-2*f5*f1^2*s1
^3*s2^7-6*f5*f1*f4*s1^6*s2^4-28*f5*f0^2*s1*s2^9-20*f5*f0*f4*s1^5*s2^5+11*f5*f0*
f3*s1^4*s2^6+5*f5*f4*f3*s1^8*s2^2-f1^2*f4*s1^2*s2^8+3*f1*f0*f3*s2^10-13*f0^2*f4
*s2^10-10*f0*f4^2*s1^4*s2^6-f0*f3^2*s1^2*s2^8-5*f4^3*s1^8*s2^2+2*f2^2*s1^2*s2^6
+6*f2*f6*s1^6*s2^2+3*f2*f0*s2^8+f2*f4*s1^4*s2^4+12*f6*f1*s1^5*s2^3+37*f6*f0*s1^
4*s2^4+3*f6*f4*s1^8+2*f5*f1*s1^4*s2^4+12*f5*f0*s1^3*s2^5-f5*f3*s1^6*s2^2-f1*f3*
s1^2*s2^6+6*f0*f4*s1^2*s2^6+2*f4^2*s1^6*s2^2-f4*s1^4*s2^2+f2*f5*f3*s1^6*s2^4+f1
*f4*f3*s1^4*s2^6-32*f2*f5*f0*s1^3*s2^7-32*f6*f1*f4*s1^7*s2^3,

s13 = 2*s1*s2^3+f1*s2^6-f6^2*f3*s1^10+f5*f1^2*s1^2*s2^8-f5*f0^2*s2^10-f5*f3^2*
s1^6*s2^4-10*f2^3*s1*s2^9+f2*f5*s1^4*s2^4-f6*f5*s1^8+f0*f3*s2^8-f4*f3*s1^4*s2^4
+f5*s1^4*s2^2+f3*s1^2*s2^4-20*f2^2*f6*s1^5*s2^5-2*f2^2*f5*s1^4*s2^6+9*f2^2*f1*
s2^10+5*f2^2*f3*s1^2*s2^8-26*f2*f6^2*s1^9*s2+10*f2*f6*f5*s1^8*s2^2-28*f2*f6*f1*
s1^4*s2^6-128*f2*f6*f0*s1^3*s2^7-24*f2*f6*f4*s1^7*s2^3+12*f2*f6*f3*s1^6*s2^4-20
*f2*f5*f1*s1^3*s2^7-80*f2*f5*f0*s1^2*s2^8-2*f2*f5*f4*s1^6*s2^4-6*f2*f1*f4*s1^2*
s2^8+10*f2*f1*f3*s1*s2^9-56*f2*f0*f4*s1*s2^9-4*f2*f0*f3*s2^10+2*f2*f4*f3*s1^4*
s2^6-57*f6^2*f1*s1^8*s2^2-152*f6^2*f0*s1^7*s2^3+20*f6*f5*f1*s1^7*s2^3+54*f6*f5*
f0*s1^6*s2^4+4*f6*f5*f4*s1^10+4*f6*f5*f3*s1^9*s2+20*f6*f1^2*s1^3*s2^7+110*f6*f1
*f0*s1^2*s2^8-24*f6*f1*f4*s1^6*s2^4+46*f6*f1*f3*s1^5*s2^5-10*f6*f0^2*s1*s2^9-44
*f6*f0*f4*s1^5*s2^5+162*f6*f0*f3*s1^4*s2^6-10*f6*f4^2*s1^9*s2+5*f5^2*f1*s1^6*s2
^4+34*f5^2*f0*s1^5*s2^5-2*f5^2*f3*s1^8*s2^2+40*f5*f1*f0*s1*s2^9-8*f5*f1*f4*s1^5
*s2^5+6*f5*f1*f3*s1^4*s2^6+12*f5*f0*f4*s1^4*s2^6+52*f5*f0*f3*s1^3*s2^7+5*f5*f4^
2*s1^8*s2^2-2*f1^2*f4*s1*s2^9-3*f1^2*f3*s2^10+14*f1*f0*f4*s2^10-2*f1*f4^2*s1^4*
s2^6-2*f1*f3^2*s1^2*s2^8+14*f0*f4*f3*s1^2*s2^8-2*f0*f3^2*s1*s2^9+2*f4^2*f3*s1^6
*s2^4+4*f2^2*s1*s2^7+12*f2*f6*s1^5*s2^3-3*f2*f1*s2^8-2*f2*f3*s1^2*s2^6+23*f6*f1
*s1^4*s2^4+76*f6*f0*s1^3*s2^5+4*f6*f4*s1^7*s2+8*f5*f1*s1^3*s2^5+32*f5*f0*s1^2*
s2^6-2*f5*f4*s1^6*s2^2+2*f1*f4*s1^2*s2^6-2*f1*f3*s1*s2^7+16*f0*f4*s1*s2^7-2*f2*
s1*s2^5-2*f6*s1^5*s2,

s14 = s2^4-2*f2*s2^6-14*f2^3*s2^10+f6^2*s1^8+5*f2^2*s2^8-4*f6^2*f4*s1^10-20*f2
^2*f6*s1^4*s2^6-22*f2*f6^2*s1^8*s2^2-32*f2*f6*f1*s1^3*s2^7-116*f2*f6*f0*s1^2*s2
^8-12*f2*f6*f4*s1^6*s2^4-10*f2*f5*f1*s1^2*s2^8-64*f2*f5*f0*s1*s2^9+12*f2*f1*f3*
s2^10-40*f2*f0*f4*s2^10-48*f6^2*f1*s1^7*s2^3-140*f6^2*f0*s1^6*s2^4-6*f6*f5*f1*
s1^6*s2^4-32*f6*f5*f0*s1^5*s2^5+4*f6*f5*f3*s1^8*s2^2-4*f6*f1^2*s1^2*s2^8+24*f6*
f1*f0*s1*s2^9-8*f6*f1*f4*s1^5*s2^5+12*f6*f1*f3*s1^4*s2^6-2*f6*f0^2*s2^10-36*f6*
f0*f4*s1^4*s2^6+48*f6*f0*f3*s1^3*s2^7-10*f6*f4^2*s1^8*s2^2+6*f5^2*f0*s1^4*s2^6-
4*f5*f1^2*s1*s2^9+6*f5*f1*f0*s2^10-2*f5*f1*f4*s1^4*s2^6+12*f5*f0*f3*s1^2*s2^8-2
*f1^2*f4*s2^10-2*f0*f3^2*s2^10+6*f2*f6*s1^4*s2^4+8*f6*f1*s1^3*s2^5+36*f6*f0*s1^
2*s2^6+4*f6*f4*s1^6*s2^2+2*f5*f1*s1^2*s2^6+16*f5*f0*s1*s2^7-2*f1*f3*s2^8+8*f0*
f4*s2^8-2*f6*s1^4*s2^2,

s15 = 4*f6*s1^6+4*f2*s1^2*s2^4+f5^2*s1^8+f1^2*s2^8+4*f0*s2^6+f3^2*s1^4*s2^4+20
*f2^3*s1^2*s2^8+40*f2^2*f6*s1^6*s2^4-8*f2^2*f5*s1^5*s2^5+20*f2^2*f1*s1*s2^9+36*
f2^2*f0*s2^10+8*f2^2*f4*s1^4*s2^6+8*f2^2*f3*s1^3*s2^7+52*f2*f6^2*s1^10+56*f2*f6
*f5*s1^9*s2+128*f2*f6*f1*s1^5*s2^5+272*f2*f6*f0*s1^4*s2^6+24*f2*f6*f4*s1^8*s2^2
+56*f2*f6*f3*s1^7*s2^3+8*f2*f5^2*s1^8*s2^2-16*f2*f5*f4*s1^7*s2^3-4*f2*f1^2*s2^
10-16*f2*f1*f4*s1^3*s2^7-26*f2*f1*f3*s1^2*s2^8+24*f2*f0*f4*s1^2*s2^8-12*f2*f0*
f3*s1*s2^9+8*f2*f4^2*s1^6*s2^4+4*f2*f4*f3*s1^5*s2^5-2*f2*f3^2*s1^4*s2^6+108*f6^
2*f1*s1^9*s2+252*f6^2*f0*s1^8*s2^2+136*f6*f5*f1*s1^8*s2^2+352*f6*f5*f0*s1^7*s2^
3-10*f6*f5*f3*s1^10+108*f6*f1^2*s1^4*s2^6+352*f6*f1*f0*s1^3*s2^7+118*f6*f1*f3*
s1^6*s2^4+252*f6*f0^2*s1^2*s2^8+272*f6*f0*f4*s1^6*s2^4+276*f6*f0*f3*s1^5*s2^5+
36*f6*f4^2*s1^10-12*f6*f4*f3*s1^9*s2+6*f6*f3^2*s1^8*s2^2+36*f5^2*f1*s1^7*s2^3+
108*f5^2*f0*s1^6*s2^4-4*f5^2*f4*s1^10-8*f5^2*f3*s1^9*s2+36*f5*f1^2*s1^3*s2^7+
136*f5*f1*f0*s1^2*s2^8+36*f5*f1*f3*s1^5*s2^5+108*f5*f0^2*s1*s2^9+128*f5*f0*f4*
s1^5*s2^5+118*f5*f0*f3*s1^4*s2^6+20*f5*f4^2*s1^9*s2-26*f5*f4*f3*s1^8*s2^2-4*f5*
f3^2*s1^7*s2^3+8*f1^2*f4*s1^2*s2^8-8*f1^2*f3*s1*s2^9+56*f1*f0*f4*s1*s2^9-10*f1*
f0*f3*s2^10-8*f1*f4^2*s1^5*s2^5-4*f1*f3^2*s1^3*s2^7+52*f0^2*f4*s2^10+40*f0*f4^2
*s1^4*s2^6+56*f0*f4*f3*s1^3*s2^7+6*f0*f3^2*s1^2*s2^8+20*f4^3*s1^8*s2^2+8*f4^2*
f3*s1^7*s2^3-2*f4*f3^2*s1^6*s2^4-8*f2^2*s1^2*s2^6+8*f2*f6*s1^6*s2^2+8*f2*f5*s1^
5*s2^3-8*f2*f1*s1*s2^7-12*f2*f0*s2^8-4*f2*f4*s1^4*s2^4-4*f2*f3*s1^3*s2^5+16*f6*
f1*s1^5*s2^3+12*f6*f0*s1^4*s2^4-12*f6*f4*s1^8+4*f6*f3*s1^7*s2+18*f5*f1*s1^4*s2^
4+16*f5*f0*s1^3*s2^5-8*f5*f4*s1^7*s2+6*f5*f3*s1^6*s2^2+8*f1*f4*s1^3*s2^5+6*f1*
f3*s1^2*s2^6+8*f0*f4*s1^2*s2^6+4*f0*f3*s1*s2^7-8*f4^2*s1^6*s2^2-4*f4*f3*s1^5*s2
^3+4*f5*s1^5*s2+4*f1*s1*s2^5+4*f4*s1^4*s2^2+4*f3*s1^3*s2^3-2*f2*f5*f3*s1^6*s2^4
-2*f1*f4*f3*s1^4*s2^6+64*f2*f5*f0*s1^3*s2^7+64*f6*f1*f4*s1^7*s2^3}:

##############################################################################
##
## loccoord2coord
##
## approximation of full coord-vector from local coordinates (using local
## expansions)
##
##############################################################################

loccoord2coord:=proc(s)
  local vals;
  vals:=subs(localexpansions,
             [1,s1,s2,s3,s4,s5,s6,s7,s8,s9,s10,s11,s12,s13,s14,s15]);
  RETURN(subs({s1=s[1],s2=s[2]},vals));
end;  

##############################################################################
##
## loccoord2kummer 
##
## approximation of kummer-vector from local coordinates (using local
## expansions)
##
##############################################################################

loccoord2kummer:=proc(s)
  local vals;
  vals:=subs(localexpansions,[s14,s13,s12,s5]);
  RETURN(subs({s1=s[1],s2=s[2]},vals));
end;

##############################################################################
##
## loclog
##
## approximation to logarithm from local coords
##
##############################################################################

loclog:=proc(s)
  local s1,s2,LOG1,LOG2;

  s1:=s[1];s2:=s[2];

LOG1 := s1-20/7*s1^7*f4^3+6/5*s1^5*f4^2+1/3*f1*s2^3-2/3*s1^3*f4-4/7*f2*s1^7*f5^
2+10/7*f2^2*f1*s2^7-f2*f1*s2^4*s1^3*f5+2*f2*f1*s2^3*f6*s1^4-8/7*f2*s2^7*f0*f3-4
*f2*s2^6*s1*f0*f4-4*f2*s2^5*s1^2*f0*f5-4*f2*s2^4*f6*s1^3*f0-64/7*f2*f6*s1^7*f4-
4/7*f1^2*s2^7*f3+f1^2*s2^5*s1^2*f5+6*f1^2*s2^4*f6*s1^3+20/7*f1*s2^7*f0*f4+7*f1*
s2^6*s1*f0*f5+18*f1*s2^5*f6*s1^2*f0+4*f1*s2^2*f6*s1^5*f3-3*f1*s2^2*s1^5*f4*f5-
14*f1*s2*f6*s1^6*f4-f1*s2*s1^6*f5^2+13/7*f1*f6*s1^7*f5-4/7*s2^7*f0^2*f5-4*s2^6*
f6*s1*f0^2-4*s2^4*s1^3*f0*f4^2+4*s2^4*s1^3*f0*f5*f3+16*s2^3*f6*s1^4*f0*f3-8*s2^
3*s1^4*f0*f4*f5-24*s2^2*f6*s1^5*f0*f4+8*s2*f6*s1^6*f0*f5-4/7*f6^2*s1^7*f0-4/7*
f6*s1^7*f3^2+20/7*s1^7*f4*f5*f3-3/5*f2*f1*s2^5+12/5*f2*f6*s1^5+f1*s2^2*s1^3*f5+
5*f1*s2*f6*s1^4+2/5*s2^5*f0*f3+2*s2^4*s1*f0*f4+4*s2^3*s1^2*f0*f5+12*s2^2*f6*s1^
3*f0-3/5*s1^5*f5*f3:

LOG2 := s2-20/7*f2^3*s2^7+6/5*f2^2*s2^5-2/3*f2*s2^3+1/3*s1^3*f5-4*f2^2*s2^3*f6*
s1^4+20/7*f2*f1*s2^7*f3-3*f2*f1*s2^5*s1^2*f5-8*f2*f1*s2^4*f6*s1^3-64/7*f2*s2^7*
f0*f4-14*f2*s2^6*s1*f0*f5-24*f2*s2^5*f6*s1^2*f0-4*f2*s2*f6*s1^6*f4+20/7*f2*f6*
s1^7*f5-4/7*f1^2*s2^7*f4-f1^2*s2^6*s1*f5+13/7*f1*s2^7*f0*f5+8*f1*s2^6*f6*s1*f0+
4*f1*s2^3*f6*s1^4*f3-f1*s2^3*s1^4*f4*f5-4*f1*s2^2*f6*s1^5*f4+f1*s2^2*s1^5*f5^2+
7*f1*s2*f6*s1^6*f5-4/7*f1*f6^2*s1^7-4/7*s2^7*f6*f0^2-4/7*s2^7*f0*f3^2+4*s2^5*s1
^2*f0*f5*f3+16*s2^4*f6*s1^3*f0*f3+2*s2^4*s1^3*f0*f4*f5-4*s2^3*f6*s1^4*f0*f4+6*
s2^3*s1^4*f0*f5^2+18*s2^2*f6*s1^5*f0*f5-4*s2*f6^2*s1^6*f0-8/7*f6*s1^7*f4*f3+10/
7*s1^7*f4^2*f5-4/7*s1^7*f5^2*f3+2*f2*s2*f6*s1^4-3/5*f1*s2^5*f3+f1*s2^3*s1^2*f5+
4*f1*s2^2*f6*s1^3+12/5*s2^5*f0*f4+5*s2^4*s1*f0*f5+12*s2^3*f6*s1^2*f0+2/5*f6*s1^
5*f3-3/5*s1^5*f4*f5:

  RETURN([LOG1,LOG2]);
end;

##############################################################################
##
## locexp
##
## approximation to inverse log.
##
##############################################################################

locexp:=proc(s)
  local EXP1,EXP2,s1,s2;
  s1:=s[1];s2:=s[2];
  EXP1 := s1+4/315*s1^7*f4^3+2/15*s1^5*f4^2-1/3*f1*s2^3+2/3*s1^3*f4+4/7*f2*s1^7*
f5^2-2/315*f2^2*f1*s2^7-2/9*f2*f1*s2^4*s1^3*f5+2/3*f2*f1*s2^3*f6*s1^4-4/21*f2*
s2^7*f0*f3-4/3*f2*s2^6*s1*f0*f4-4*f2*s2^5*s1^2*f0*f5-12*f2*s2^4*f6*s1^3*f0-128/
35*f2*f6*s1^7*f4-1/35*f1^2*s2^7*f3+2/3*f1^2*s2^5*s1^2*f5+14/3*f1^2*s2^4*f6*s1^3
+22/105*f1*s2^7*f0*f4+2/3*f1*s2^6*s1*f0*f5+6*f1*s2^5*f6*s1^2*f0-18/5*f1*s2^2*f6
*s1^5*f3-4/15*f1*s2^2*s1^5*f4*f5-28/3*f1*s2*f6*s1^6*f4+14/9*f1*s2*s1^6*f5^2-4/
21*f1*f6*s1^7*f5+4/7*s2^7*f0^2*f5+4*s2^6*f6*s1*f0^2-4/3*s2^4*s1^3*f0*f4^2-10/3*
s2^4*s1^3*f0*f5*f3-16*s2^3*f6*s1^4*f0*f3-8/3*s2^3*s1^4*f0*f4*f5-24*s2^2*f6*s1^5
*f0*f4+4/7*f6^2*s1^7*f0+4/7*f6*s1^7*f3^2+12/35*s1^7*f4*f5*f3-1/15*f2*f1*s2^5-12
/5*f2*f6*s1^5-2/3*f1*s2^2*s1^3*f5-5*f1*s2*f6*s1^4-2/5*s2^5*f0*f3-2*s2^4*s1*f0*
f4-4*s2^3*s1^2*f0*f5-12*s2^2*f6*s1^3*f0+3/5*s1^5*f5*f3-2/15*f2*f1*f4*s1^2*s2^5+
2/9*f1^2*f4*s1*s2^6-2/9*f1*f4^2*s1^4*s2^3-f1*f5*f3*s1^4*s2^3-4/5*f0*f4*f3*s1^2*
s2^5+4*f0*f5^2*s1^5*s2^2-2/3*f1*f4*s1^2*s2^3:
EXP2 := s2+4/315*f2^3*s2^7+2/15*f2^2*s2^5+2/3*f2*s2^3-1/3*s1^3*f5-4/3*f2^2*s2^3
*f6*s1^4+12/35*f2*f1*s2^7*f3-4/15*f2*f1*s2^5*s1^2*f5-8/3*f2*f1*s2^4*f6*s1^3-128
/35*f2*s2^7*f0*f4-28/3*f2*s2^6*s1*f0*f5-24*f2*s2^5*f6*s1^2*f0-4/3*f2*s2*f6*s1^6
*f4+22/105*f2*f6*s1^7*f5+4/7*f1^2*s2^7*f4+14/9*f1^2*s2^6*s1*f5-4/21*f1*s2^7*f0*
f5-10/3*f1*s2^3*f6*s1^4*f3-2/9*f1*s2^3*s1^4*f4*f5-4*f1*s2^2*f6*s1^5*f4+2/3*f1*
s2^2*s1^5*f5^2+2/3*f1*s2*f6*s1^6*f5+4/7*f1*f6^2*s1^7+4/7*s2^7*f6*f0^2+4/7*s2^7*
f0*f3^2-18/5*s2^5*s1^2*f0*f5*f3-16*s2^4*f6*s1^3*f0*f3+2/3*s2^4*s1^3*f0*f4*f5-12
*s2^3*f6*s1^4*f0*f4+14/3*s2^3*s1^4*f0*f5^2+6*s2^2*f6*s1^5*f0*f5+4*s2*f6^2*s1^6*
f0-4/21*f6*s1^7*f4*f3-2/315*s1^7*f4^2*f5-1/35*s1^7*f5^2*f3-2*f2*s2*f6*s1^4+3/5*
f1*s2^5*f3-2/3*f1*s2^3*s1^2*f5-4*f1*s2^2*f6*s1^3-12/5*s2^5*f0*f4-5*s2^4*s1*f0*
f5-12*s2^3*f6*s1^2*f0-2/5*f6*s1^5*f3-1/15*s1^5*f4*f5-2/9*f2^2*f5*s1^3*s2^4-4/5*
f2*f6*f3*s1^5*s2^2-2/15*f2*f4*f5*s1^5*s2^2+2/9*f2*f5^2*s1^6*s2+4*f1^2*f6*s1^2*
s2^5-f1*f5*f3*s1^3*s2^4-2/3*f2*f5*s1^3*s2^2:

  RETURN([EXP1,EXP2]);
end;

##############################################################################
##
## bilinear forms, projectively equivalent to kummer of sum of coord reps
## a and b
##
##############################################################################

bilinform:=array(1..4,1..4);
bilinform[1,1] :=
a0*b14-2*a5*b5+a14*b0+

4*f5*a12*b4+8*f6*a11*b4-4*f6*a5*b10-2*f5*a3*b13+4*f5*a4*b12-2*f5*a13*b3-12*f6*
a3*b12-4*f6*a3*b15+8*f6*a4*b11-4*f6*a10*b5-12*f6*a12*b3-4*f6*b3*a15+

(4*f5*f3+8*f6*f2)*a12*b12+(f5**2-4*f4*f6)*b15*a10+(f5**2-4*f4*f6)*a15*b10+(2*f1
*f5-16*f0*f6)*a12*b14+(2*f1*f5-16*f0*f6)*a14*b12+(-2*f5**2+8*f4*f6)*a11*b11+(6*
f5**2-16*f4*f6)*a10*b12+(6*f5**2-16*f4*f6)*a12*b10+4*f5*f6*a10*b11+4*f5*f6*a11*
b10+4*f1*f6*b13*a12+4*f1*f6*a13*b12+8*f3*f6*a12*b11-4*f3*f6*a13*b10-4*f2*f6*a14
*b10+8*f3*f6*a11*b12+8*f0*f6*a13*b13-4*f3*f6*a10*b13-4*f2*f6*a10*b14-4*f0*f6*
b14*a15-4*f0*f6*a14*b15
;

bilinform[1,2] :=
b13*a0+2*a1*b9-2*a2*b8-2*a4*b5-2*a5*b4-2*a8*b2+2*a9*b1+a13*b0+

2*f4*b13*a3-4*f4*a8*b7+2*f4*a13*b3+f5*b3*a15-4*f5*a7*b7-2*f5*a11*b4+2*f5*b3*a12
+2*f5*a3*b12+f5*a5*b10+f5*a10*b5-2*f5*a4*b11+4*f6*a10*b4-4*f4*a7*b8-2*f5*a8*b6-
4*f3*a8*b8-2*f5*a6*b8+4*f6*a4*b10-2*f6*a3*b11-4*f6*a6*b7-4*f6*a7*b6-2*f6*a11*b3
-f1*a14*b5-f1*b14*a5+f3*a3*b14+f3*a14*b3-4*f4*a12*b4-2*f3*a12*b5-2*f3*a5*b12-4*
f4*a4*b12+f5*a3*b15+

(4*f1*f6+2*f2*f5)*a14*b10+(4*f6*f2+f5*f3)*a13*b10+(4*f6*f2+f5*f3)*b13*a10-f5**2
*a11*b10-f5**2*a10*b11+4*f5*f6*a10*b10-4*f3*f6*a11*b11+(4*f0*f6+f1*f5)*b14*a11+
(4*f0*f6+f1*f5)*b11*a14+(4*f1*f6+2*f2*f5)*b14*a10+2*f1*f6*a13*b11+2*f1*f6*a11*
b13-4*f0*f5*a13*b13+2*b10*f6*f3*a15-2*f1*f6*b12*a15+2*b14*f5*f0*a15+2*a10*f6*f3
*b15-2*f1*f6*a12*b15+2*a14*f5*f0*b15+(-4*f2*f5-20*f1*f6-4*f3*f4)*a12*b12+(8*f5*
f0-2*f1*f4)*a12*b14+(8*f5*f0-2*f1*f4)*a14*b12+(10*f3*f6-2*f4*f5)*a10*b12+(10*f3
*f6-2*f4*f5)*a12*b10+(-4*f0*f6-2*f1*f5)*a12*b13+(-4*f0*f6-2*f1*f5)*b12*a13+(-4*
f6*f2-2*f5*f3)*b11*a12+(-4*f6*f2-2*f5*f3)*b12*a11
;

bilinform[1,3] :=
a0*b12+b0*a12-2*a7*b2-4*a4*b4+a5*b3+a3*b5-2*b7*a2+2*b8*a1+2*a8*b1+

f5*b11*a3-2*f5*b10*a4+f5*a11*b3-2*f5*a10*b4-2*f5*a7*b6-2*f5*a6*b7+2*f6*b10*a3+2
*f6*a10*b3-4*f6*a6*b6-2*f1*a14*b4+f1*a13*b5+f1*b13*a5-2*f1*a4*b14+2*f1*a8*b9+2*
f1*a9*b8+2*f2*a12*b5+2*f2*b12*a5+4*f2*a8*b8+2*f4*a12*b3-4*f4*a7*b7+2*f4*b12*a3+
2*f0*a14*b5+2*f0*a5*b14+4*f0*a9*b9+

(4*f0*f2-2*f1**2)*a14*b14+(-4*f0*f6+f1*f5)*a11*b13+(-2*f5**2+4*f4*f6)*a10*b10+4
*f2*f6*a11*b11+2*f3*f6*a10*b11+2*f3*f6*a11*b10+(-4*f0*f6+f1*f5)*a13*b11-2*f5*f0
*a11*b14-2*f5*f0*a14*b11+2*f3*f0*a13*b14+2*f3*f0*a14*b13+f3*f5*a12*b10-2*f1*f6*
a13*b10-2*f1*f5*a14*b10+f3*f5*a10*b12+f1*f3*a14*b12-2*f1*f6*a10*b13+4*f0*f4*a13
*b13-2*f1*f5*a10*b14+f1*f3*a12*b14+2*f1*f6*b11*a15+2*f0*f5*b13*a15+2*f1*f6*a11*
b15+2*f0*f5*a13*b15+4*f0*f6*a15*b15+(68*f0*f6+8*f1*f5+4*f4*f2)*a12*b12+(2*f1*f4
+8*f5*f0)*a12*b13+(2*f1*f4+8*f5*f0)*b12*a13+(8*f1*f6+2*f2*f5)*b11*a12+(8*f1*f6+
2*f2*f5)*b12*a11+(16*f0*f6+f1*f5)*a15*b12+(16*f0*f6+f1*f5)*b15*a12
;

bilinform[1,4] :=
-2*a2*b2+b0*a5+a0*b5+

4*a6*b1*f6+4*b6*a1*f6-6*a3*b3*f6+2*a5*b5*f2+2*f5*a7*b1-2*a3*b4*f5-2*b3*a4*f5+2*
f5*b7*a1+

(-4*f4*f0+f1*f3)*a14*b5+(-2*f1*f5-8*f0*f6)*a13*b4+(-2*f1*f5-8*f0*f6)*b13*a4+(4*
f2*f5+8*f1*f6)*b8*a7+(4*f2*f5+8*f1*f6)*a8*b7+(4*f5*f3+8*f6*f2)*b7*a7+(8*f0*f6+2
*f1*f5)*a7*b9+(-8*f4*f6+f5**2)*a10*b3+(-8*f4*f6+f5**2)*b10*a3+(-2*f5**2+8*f4*f6
)*b6*a6+(8*f0*f6+2*f1*f5)*b7*a9+(4*f1*f5+8*f0*f6)*b8*a8-4*b3*a11*f6*f3-4*a4*b10
*f6*f3-8*b4*a11*f6*f2-8*a4*b11*f6*f2+8*f2*f6*a6*b8+8*f3*f6*b6*a7-4*b4*a10*f6*f3
-4*a3*b11*f6*f3+8*f3*f6*a6*b7+4*f5*f0*a8*b9+8*f2*f6*b6*a8+4*f1*f6*b9*a6+4*f1*f6
*a9*b6-4*b5*a11*f6*f1-4*a5*b11*f6*f1+4*f5*f0*b8*a9+(-f1*f5-4*f0*f6)*a3*b14+(-4*
f6*f2-2*f5*f3)*a12*b3+(-f1*f5-4*f0*f6)*b3*a14+(-4*f6*f2-2*f5*f3)*b12*a3+(-4*f4*
f0+f1*f3)*b14*a5-2*b3*a13*f6*f1-6*b5*a13*f5*f0-4*b4*a14*f5*f0-2*a3*b13*f6*f1-6*
a5*b13*f5*f0-4*a4*b14*f5*f0-8*b5*f6*f0*a15-4*b4*f6*f1*a15-8*a5*f6*f0*b15-4*a4*
f6*f1*b15+(-16*f1*f6-4*f2*f5)*b12*a4+(-16*f1*f6-4*f2*f5)*b4*a12+(-28*f0*f6-2*f1
*f5)*b5*a12+(-28*f0*f6-2*f1*f5)*b12*a5+

(-8*f0*f5*f6-f1*f5**2)*a10*b13+(-4*f1*f3*f6-2*f0*f5**2)*b13*a11+(-4*f0*f3*f6-4*
f6*f2*f1)*a11*b14+(-4*f1*f5*f6-8*f2*f4*f6-8*f6**2*f0+2*f2*f5**2)*b11*a11+(-4*f2
*f5*f6+f5**2*f3-4*f6*f4*f3-8*f1*f6**2)*a10*b11+(-4*f2*f5*f6+f5**2*f3-4*f6*f4*f3
-8*f1*f6**2)*b10*a11+(-8*f6*f4**2+2*f5**2*f4-8*f6**2*f2)*b10*a10+(-8*f1*f0*f6-8
*f0*f2*f5+f1**2*f5)*a14*b13+(-4*f0*f3*f6-4*f6*f2*f1)*b11*a14+(2*f6*f1**2-4*f0*
f3*f5-16*f0*f2*f6)*b13*a13+(-4*f1*f3*f6-2*f0*f5**2)*a13*b11+(-8*f0*f5*f6-f1*f5
**2)*b10*a13+(-8*f1*f0*f6-8*f0*f2*f5+f1**2*f5)*b14*a13+(2*f1**2*f4-4*f0*f1*f5-8
*f6*f0**2-8*f0*f2*f4+2*f0*f3**2)*b14*a14-8*b10*a14*f6*f4*f0-8*a10*b14*f6*f4*f0-
8*b13*f6*f3*f0*a15-8*b14*f6*f2*f0*a15-8*a13*f6*f3*f0*b15-8*a14*f6*f2*f0*b15+(-8
*f6*f2**2-4*f5*f1*f4+20*f0*f5**2-136*f0*f4*f6-24*f1*f3*f6-4*f3*f5*f2)*a12*b12+(
-8*f6**2*f0-2*f1*f5*f6)*b15*a10+(-8*f6**2*f0-2*f1*f5*f6)*a15*b10+(-6*f0*f3*f5-4
*f6*f1**2-2*f1*f2*f5-24*f0*f2*f6)*a12*b14+(-6*f0*f3*f5-4*f6*f1**2-2*f1*f2*f5-24
*f0*f2*f6)*a14*b12+(-12*f1*f5*f6-24*f6**2*f0-4*f3**2*f6-2*f2*f5**2)*a10*b12+(-
12*f1*f5*f6-24*f6**2*f0-4*f3**2*f6-2*f2*f5**2)*a12*b10+(-4*f6*f2*f1-2*f1*f3*f5-
4*f5*f4*f0-32*f0*f3*f6)*a12*b13+(-4*f6*f2*f1-2*f1*f3*f5-4*f5*f4*f0-32*f0*f3*f6)
*b12*a13+(-16*f0*f5*f6-8*f2*f3*f6+2*f1*f5**2-16*f1*f4*f6)*b11*a12+(-16*f0*f5*f6
-8*f2*f3*f6+2*f1*f5**2-16*f1*f4*f6)*b12*a11+(-4*f1*f3*f6+6*f0*f5**2-32*f0*f4*f6
)*a15*b12+(-4*f1*f3*f6+6*f0*f5**2-32*f0*f4*f6)*b15*a12+(-8*f0*f4*f6+2*f0*f5**2)
*a15*b15+(-4*f1*f4*f6+f1*f5**2-4*f0*f5*f6)*a15*b11+(-4*f1*f4*f6+f1*f5**2-4*f0*
f5*f6)*b15*a11
;

bilinform[2,1] :=
b13*a0-2*a1*b9+2*a2*b8-2*a4*b5-2*a5*b4+2*a8*b2-2*a9*b1+a13*b0+

2*f4*b13*a3+4*f4*a8*b7+2*f4*a13*b3+f5*b3*a15+4*f5*a7*b7-2*f5*a11*b4+2*f5*b3*a12
+2*f5*a3*b12+f5*a5*b10+f5*a10*b5-2*f5*a4*b11+4*f6*a10*b4+4*f4*a7*b8+2*f5*a8*b6+
4*f3*a8*b8+2*f5*a6*b8+4*f6*a4*b10-2*f6*a3*b11+4*f6*a6*b7+4*f6*a7*b6-2*f6*a11*b3
-f1*a14*b5-f1*b14*a5+f3*a3*b14+f3*a14*b3-4*f4*a12*b4-2*f3*a12*b5-2*f3*a5*b12-4*
f4*a4*b12+f5*a3*b15+

(4*f1*f6+2*f2*f5)*a14*b10+(4*f6*f2+f5*f3)*a13*b10+(4*f6*f2+f5*f3)*b13*a10-f5**2
*a11*b10-f5**2*a10*b11+4*f5*f6*a10*b10-4*f3*f6*a11*b11+(4*f0*f6+f1*f5)*b14*a11+
(4*f0*f6+f1*f5)*b11*a14+(4*f1*f6+2*f2*f5)*b14*a10+2*f1*f6*a13*b11+2*f1*f6*a11*
b13-4*f0*f5*a13*b13+2*b10*f6*f3*a15-2*f1*f6*b12*a15+2*b14*f5*f0*a15+2*a10*f6*f3
*b15-2*f1*f6*a12*b15+2*a14*f5*f0*b15+(-4*f2*f5-20*f1*f6-4*f3*f4)*a12*b12+(8*f5*
f0-2*f1*f4)*a12*b14+(8*f5*f0-2*f1*f4)*a14*b12+(10*f3*f6-2*f4*f5)*a10*b12+(10*f3
*f6-2*f4*f5)*a12*b10+(-4*f0*f6-2*f1*f5)*a12*b13+(-4*f0*f6-2*f1*f5)*b12*a13+(-4*
f6*f2-2*f5*f3)*b11*a12+(-4*f6*f2-2*f5*f3)*b12*a11
;

bilinform[2,2] :=
4*a0*b12+a15*b0+4*b0*a12+a0*b15-4*a5*b3-4*a3*b5+

4*f3*a12*b4-4*f5*a10*b4-4*f6*b10*a3-4*f6*a10*b3-4*f0*a5*b14-4*f1*a14*b4-4*f1*a4
*b14-4*f2*a14*b3-4*f2*a3*b14-2*f3*a5*b11+4*f3*a4*b12-2*f3*a3*b13-4*f5*b10*a4-2*
f3*a13*b3-2*f3*a11*b5-4*f4*a10*b5-4*f4*a5*b10-4*f0*a14*b5

-4*f2*f6*a10*b15-24*f0*f6*a12*b15-24*f0*f6*a15*b12-4*f2*f6*a15*b10-8*f1*f6*a12*
b11-8*f1*f6*a11*b12-4*f3*f6*a10*b11-4*f3*f6*a11*b10-8*f5*f0*a12*b13-8*f5*f0*a13
*b12-4*f3*f0*a13*b14-4*f3*f0*a14*b13-4*b13*a11*f1*f5-4*b11*a13*f1*f5-4*f1*f6*
b11*a15-4*f0*f5*b13*a15-4*f1*f6*a11*b15-4*f0*f5*a13*b15-8*f0*f6*a15*b15-4*b14*
f4*f0*a15-4*f4*f0*a14*b15+(4*f3**2+8*f1*f5-64*f0*f6)*a12*b12+(-16*f4*f0-2*f1*f3
)*a12*b14+(-16*f4*f0-2*f1*f3)*a14*b12+(-16*f6*f2-2*f5*f3)*a10*b12+(-16*f6*f2-2*
f5*f3)*a12*b10+(-4*f1*f6-4*f2*f5)*b10*a13+(f3**2-4*f0*f6-4*f4*f2-6*f1*f5)*b10*
a14+(-4*f1*f6-4*f2*f5)*b13*a10+(f3**2-4*f0*f6-4*f4*f2-6*f1*f5)*b14*a10+(-4*f5*
f0-4*f1*f4)*b14*a11+(-4*f5*f0-4*f1*f4)*b11*a14+(-8*f4*f6-2*f5**2)*b10*a10+(-8*
f0*f2-2*f1**2)*a14*b14
;
bilinform[2,3] :=
-2*a3*b4-2*a4*b3+2*a1*b7-2*a2*b6+a0*b11-2*a6*b2+2*a7*b1+a11*b0+

4*f1*a8*b8-2*f1*a4*b13+2*f2*a5*b11+2*f1*a7*b9-2*f0*a13*b5+4*f3*a7*b7+f1*a5*b15+
f1*a15*b5+2*f2*a11*b5-2*f1*a13*b4+f3*a10*b5+2*f1*a12*b5-f5*a3*b10-f5*a10*b3+f3*
a5*b10+2*f1*a5*b12+2*f1*a9*b7+4*f0*a9*b8+4*f0*a8*b9+4*f2*a8*b7+4*f2*a7*b8+4*f0*
a4*b14-2*f0*a5*b13+4*f0*a14*b4-4*f2*a12*b4-2*f3*a12*b3+f1*a14*b3-2*f3*a3*b12-4*
f2*a4*b12+f1*a3*b14+

(-2*f1*f3-4*f4*f0)*a12*b13+(2*f1*f4+4*f5*f0)*b14*a10-2*b15*f5*f0*a12-2*a15*f5*
f0*b12+2*f0*f3*a14*b15+2*f0*f3*a15*b14-4*f6*f1*a11*b11+(-4*f2*f3-20*f5*f0-4*f1*
f4)*a12*b12+(-4*f0*f6-2*f1*f5)*a11*b12+(-4*f0*f6-2*f1*f5)*a12*b11+2*f5*f0*a11*
b13+2*f5*f0*a13*b11+(4*f0*f6+f1*f5)*b10*a13+(2*f1*f4+4*f5*f0)*a14*b10+(f1*f3+4*
f4*f0)*a14*b11+(-2*f1*f3-4*f4*f0)*a13*b12+(4*f0*f6+f1*f5)*a10*b13-f1**2*a14*b13
+(f1*f3+4*f4*f0)*a11*b14-f1**2*a13*b14-4*f3*f0*a13*b13+4*f1*f0*a14*b14+2*f6*f1*
a10*b15+2*f6*f1*b10*a15+(10*f3*f0-2*f2*f1)*a12*b14+(10*f3*f0-2*f2*f1)*a14*b12+(
8*f1*f6-2*f2*f5)*a10*b12+(8*f1*f6-2*f2*f5)*a12*b10
;
bilinform[2,4] :=
2*b0*a4+2*a0*b4-2*a1*b2-2*b1*a2

-4*f2*a8*b2-4*a1*f4*b7+f3*a12*b0+b0*f1*a14-2*f3*a8*b1+b3*f3*a5+a0*f5*b10+4*a3*
f4*b4-2*a1*f5*b6-4*b1*f4*a7+2*a3*b3*f5-2*f3*a7*b2+a0*f3*b12+f5*a10*b0+f3*a3*b5+
4*f3*a4*b4-4*f2*b8*a2+2*a5*b5*f1-2*f1*b9*a2-2*f3*b8*a1+4*b3*f4*a4-2*a9*b2*f1-2*
f3*b7*a2+4*a5*b4*f2+4*b5*a4*f2-2*b1*f5*a6+a0*f1*b14
+
(4*f3*f0+2*f2*f1)*a14*b5+(-4*f0*f6+2*f1*f5)*a11*b5+(-2*f1*f3-8*f4*f0)*a8*b9+(4*
f6*f2+f5*f3)*a11*b3+(-4*f0*f6+2*f1*f5)*b11*a5+(2*f4*f5+4*f3*f6)*b10*a3+(-2*f1*
f5-8*f0*f6)*a9*b6+(-2*f1*f3-8*f4*f0)*b8*a9+(4*f6*f2+f5*f3)*b11*a3+(-2*f1*f5-8*
f0*f6)*b9*a6+(2*f4*f5+4*f3*f6)*a10*b3+(-8*f1*f5-8*f4*f2-8*f0*f6)*b8*a7+(-8*f1*
f5-8*f4*f2-8*f0*f6)*a8*b7+(-8*f5*f0-4*f1*f4)*a7*b9+(-8*f5*f0-4*f1*f4)*b7*a9+(-8
*f6*f2-2*f5*f3)*b6*a7+(-8*f1*f6-8*f2*f5-4*f3*f4)*b7*a7+(-8*f6*f2-2*f5*f3)*a6*b7
+(-8*f1*f4-8*f5*f0-4*f2*f3)*b8*a8+(-8*f1*f6-4*f2*f5)*b6*a8+(-8*f1*f6-4*f2*f5)*
a6*b8-4*f3*f0*a9*b9+4*a4*b11*f5*f2-4*f3*f6*a6*b6+4*b10*f3*f5*a4+4*f3*f5*a10*b4+
2*b5*a10*f5*f2+4*b4*a11*f5*f2+2*a5*b10*f5*f2+(4*f1*f6+4*f2*f5+2*f3*f4)*b3*a12+(
4*f1*f4+4*f5*f0+2*f2*f3)*b12*a5+(-4*f0*f6+2*f1*f5)*a13*b3+(f1*f3+4*f4*f0)*b5*
a13+(-4*f0*f6+2*f1*f5)*b13*a3+(f1*f3+4*f4*f0)*a5*b13+(4*f3*f0+2*f2*f1)*b14*a5+4
*b4*a13*f4*f1+4*b4*a14*f1*f3+2*b3*f4*f1*a14+4*a4*b13*f4*f1+4*a4*b14*f1*f3+2*a3*
f4*f1*b14+2*b5*f0*f5*a15+2*b4*f1*f5*a15+2*b3*f1*f6*a15+2*a3*f1*f6*b15+2*a4*f1*
f5*b15+2*a5*f0*f5*b15+(4*f1*f6+4*f2*f5+2*f3*f4)*b12*a3+(8*f1*f5+8*f4*f2-8*f0*f6
)*b12*a4+(8*f1*f5+8*f4*f2-8*f0*f6)*b4*a12+(4*f1*f4+4*f5*f0+2*f2*f3)*b5*a12+

(4*f5*f4*f0-12*f0*f3*f6+4*f6*f2*f1+f1*f3*f5)*a13*b11+(4*f0*f5**2+2*f5*f1*f4-8*
f0*f4*f6)*a10*b13+(4*f5*f4*f0-12*f0*f3*f6+4*f6*f2*f1+f1*f3*f5)*a11*b13+(4*f0*f5
*f6+f1*f5**2)*b15*a10+(4*f1*f5*f6+2*f2*f5**2+8*f6**2*f0+2*f3**2*f6)*a10*b11+(8*
f1*f6**2+2*f5**2*f3+4*f6*f4*f3)*b10*a10+(4*f2*f3*f6+8*f0*f5*f6+2*f1*f5**2)*b11*
a11+(4*f1*f5*f6+2*f2*f5**2+8*f6**2*f0+2*f3**2*f6)*b10*a11+(8*f4**2*f0-8*f0*f2*
f6+4*f1*f2*f5+12*f6*f1**2+2*f4*f3*f1)*b12*a13+(32*f6*f2*f1+8*f2**2*f5+8*f1*f4**
2+4*f3*f2*f4-12*f0*f3*f6+32*f5*f4*f0)*a12*b12+2*b10*f3*f5*f1*a14+(4*f0*f5**2+2*
f5*f1*f4-8*f0*f4*f6)*b10*a13+(4*f6*f1**2-8*f0*f2*f6+2*f1*f2*f5)*b11*a14+(2*f1**
2*f5+4*f4*f3*f0+8*f1*f0*f6)*b13*a13+(2*f1**2*f4+8*f6*f0**2+2*f0*f3**2+4*f0*f1*
f5)*a14*b13+(4*f6*f1**2-8*f0*f2*f6+2*f1*f2*f5)*a11*b14+(2*f1**2*f4+8*f6*f0**2+2
*f0*f3**2+4*f0*f1*f5)*b14*a13+(4*f0*f2*f3+2*f1**2*f3+8*f0**2*f5)*a14*b14+2*a10*
f5*f1*f3*b14+4*f0*f3*f6*a15*b15+(4*f0*f5*f6+f1*f5**2)*a15*b10+(12*f1*f0*f6-8*f0
*f2*f5+4*f4*f3*f0+8*f1**2*f5+4*f1*f4*f2+f1*f3**2)*a12*b14+(12*f1*f0*f6-8*f0*f2*
f5+4*f4*f3*f0+8*f1**2*f5+4*f1*f4*f2+f1*f3**2)*a14*b12+(12*f0*f5*f6+8*f1*f5**2+4
*f2*f3*f6+4*f5*f2*f4+f3**2*f5-8*f1*f4*f6)*a10*b12+(12*f0*f5*f6+8*f1*f5**2+4*f2*
f3*f6+4*f5*f2*f4+f3**2*f5-8*f1*f4*f6)*a12*b10+(8*f4**2*f0-8*f0*f2*f6+4*f1*f2*f5
+12*f6*f1**2+2*f4*f3*f1)*a12*b13+(8*f6*f2**2-8*f0*f4*f6+4*f5*f1*f4+12*f0*f5**2+
2*f3*f5*f2)*b11*a12+(8*f6*f2**2-8*f0*f4*f6+4*f5*f1*f4+12*f0*f5**2+2*f3*f5*f2)*
b12*a11+(4*f1*f0*f6+f1**2*f5)*a15*b14+(4*f1*f0*f6+f1**2*f5)*b15*a14+(4*f6*f2*f1
+f1*f3*f5+4*f5*f4*f0+8*f0*f3*f6)*a15*b12+(4*f6*f2*f1+f1*f3*f5+4*f5*f4*f0+8*f0*
f3*f6)*b15*a12+(2*f1*f3*f6+2*f0*f5**2)*a15*b11+(2*f6*f1**2+2*f0*f3*f5)*a15*b13+
(2*f1*f3*f6+2*f0*f5**2)*b15*a11+(2*f6*f1**2+2*f0*f3*f5)*b15*a13
;

bilinform[3,1] :=
b0*a12-2*a8*b1+2*a7*b2+a5*b3-4*a4*b4+a3*b5+2*b7*a2-2*b8*a1+a0*b12+

2*f0*a5*b14+2*f6*a10*b3+4*f6*a6*b6+2*f6*b10*a3+f5*b11*a3-2*f5*b10*a4+f5*a11*b3-
2*f5*a10*b4+2*f5*a7*b6+2*f5*a6*b7+2*f4*b12*a3+2*f4*a12*b3+4*f4*a7*b7-4*f2*a8*b8
+2*f2*b12*a5+2*f2*a12*b5-2*f1*a4*b14-2*f1*a8*b9-2*f1*a14*b4-2*f1*a9*b8+f1*b13*
a5+2*f0*a14*b5-4*f0*a9*b9+f1*a13*b5+

(4*f0*f2-2*f1**2)*a14*b14+(-4*f0*f6+f1*f5)*a11*b13+(-2*f5**2+4*f4*f6)*a10*b10+4
*f2*f6*a11*b11+2*f3*f6*a10*b11+2*f3*f6*a11*b10+(-4*f0*f6+f1*f5)*a13*b11-2*f5*f0
*a11*b14-2*f5*f0*a14*b11+2*f3*f0*a13*b14+2*f3*f0*a14*b13+f3*f5*a12*b10-2*f1*f6*
a13*b10-2*f1*f5*a14*b10+f3*f5*a10*b12+f1*f3*a14*b12-2*f1*f6*a10*b13+4*f0*f4*a13
*b13-2*f1*f5*a10*b14+f1*f3*a12*b14+2*f1*f6*b11*a15+2*f0*f5*b13*a15+2*f1*f6*a11*
b15+2*f0*f5*a13*b15+4*f0*f6*a15*b15+(68*f0*f6+8*f1*f5+4*f4*f2)*a12*b12+(2*f1*f4
+8*f5*f0)*a12*b13+(2*f1*f4+8*f5*f0)*b12*a13+(8*f1*f6+2*f2*f5)*b11*a12+(8*f1*f6+
2*f2*f5)*b12*a11+(16*f0*f6+f1*f5)*a15*b12+(16*f0*f6+f1*f5)*b15*a12
;

bilinform[3,2] :=
-2*a3*b4-2*a4*b3-2*a1*b7+2*a2*b6+a0*b11+2*a6*b2-2*a7*b1+a11*b0

-4*f1*a8*b8-2*f1*a4*b13+2*f2*a5*b11-2*f1*a7*b9-2*f0*a13*b5-4*f3*a7*b7+f1*a5*b15
+f1*a15*b5+2*f2*a11*b5-2*f1*a13*b4+f3*a10*b5+2*f1*a12*b5-f5*a3*b10-f5*a10*b3+f3
*a5*b10+2*f1*a5*b12-2*f1*a9*b7-4*f0*a9*b8-4*f0*a8*b9-4*f2*a8*b7-4*f2*a7*b8+4*f0
*a4*b14-2*f0*a5*b13+4*f0*a14*b4-4*f2*a12*b4-2*f3*a12*b3+f1*a14*b3-2*f3*a3*b12-4
*f2*a4*b12+f1*a3*b14+

(-2*f1*f3-4*f4*f0)*a12*b13+(2*f1*f4+4*f5*f0)*b14*a10-2*b15*f5*f0*a12-2*a15*f5*
f0*b12+2*f0*f3*a14*b15+2*f0*f3*a15*b14-4*f6*f1*a11*b11+(-4*f2*f3-20*f5*f0-4*f1*
f4)*a12*b12+(-4*f0*f6-2*f1*f5)*a11*b12+(-4*f0*f6-2*f1*f5)*a12*b11+2*f5*f0*a11*
b13+2*f5*f0*a13*b11+(4*f0*f6+f1*f5)*b10*a13+(2*f1*f4+4*f5*f0)*a14*b10+(f1*f3+4*
f4*f0)*a14*b11+(-2*f1*f3-4*f4*f0)*a13*b12+(4*f0*f6+f1*f5)*a10*b13-f1**2*a14*b13
+(f1*f3+4*f4*f0)*a11*b14-f1**2*a13*b14-4*f3*f0*a13*b13+4*f1*f0*a14*b14+2*f6*f1*
a10*b15+2*f6*f1*b10*a15+(10*f3*f0-2*f2*f1)*a12*b14+(10*f3*f0-2*f2*f1)*a14*b12+(
8*f1*f6-2*f2*f5)*a10*b12+(8*f1*f6-2*f2*f5)*a12*b10
;

bilinform[3,3] :=
a0*b10-2*a3*b3+a10*b0+

4*f1*a4*b12-2*f1*a5*b11-2*f1*a11*b5+4*f1*a12*b4-4*f0*a3*b14+8*f0*a4*b13-12*f0*
a5*b12-4*f0*b15*a5-4*f0*a15*b5-12*f0*a12*b5+8*f0*a13*b4-4*f0*a14*b3+

(4*f1*f3+8*f4*f0)*a12*b12-4*f6*f0*a10*b15-4*f6*f0*b10*a15+(-16*f0*f2+6*f1**2)*
a12*b14+(-16*f0*f2+6*f1**2)*a14*b12+8*f6*f0*a11*b11+(2*f1*f5-16*f0*f6)*a10*b12+
(2*f1*f5-16*f0*f6)*a12*b10+8*f3*f0*a12*b13+8*f3*f0*a13*b12+4*f5*f0*a12*b11-4*f4
*f0*a14*b10+4*f5*f0*a11*b12+(-2*f1**2+8*f0*f2)*a13*b13-4*f4*f0*a10*b14+(f1**2-4
*f0*f2)*a15*b14+(f1**2-4*f0*f2)*b15*a14-4*f3*f0*a11*b14-4*f3*f0*a14*b11+4*f1*f0
*a13*b14+4*f1*f0*a14*b13
;

bilinform[3,4] :=
a0*b3-2*a1*b1+b0*a3+

2*a3*b3*f4+2*f1*b8*a2-2*b5*a4*f1-2*a5*b4*f1+2*f1*a8*b2+4*b9*a2*f0-6*a5*b5*f0+4*
a9*b2*f0

-8*a4*b13*f0*f4-4*a5*b13*f0*f3-4*b5*a13*f0*f3-4*b3*a13*f0*f5-8*b4*a13*f0*f4+(f1
**2-8*f0*f2)*a14*b5+(f1**2-8*f0*f2)*b14*a5+(-2*f1*f3-4*f4*f0)*a5*b12+(-2*f1*f3-
4*f4*f0)*b5*a12-8*b3*f0*f6*a15-4*b4*a15*f5*f0-4*a4*b15*f5*f0-8*a3*f0*f6*b15-4*
a3*b13*f0*f5-4*b4*a14*f0*f3+(-4*f6*f2+f5*f3)*a10*b3+(8*f5*f0+4*f1*f4)*b8*a7+(8*
f0*f6+2*f1*f5)*b6*a8+(-f1*f5-4*f0*f6)*b10*a5+(-2*f1*f5-8*f0*f6)*b11*a4+(-f1*f5-
4*f0*f6)*a10*b5+(4*f1*f3+8*f4*f0)*b8*a8+(-2*f1**2+8*f0*f2)*b9*a9+(4*f1*f5+8*f0*
f6)*b7*a7+(8*f0*f6+2*f1*f5)*a6*b8+(-2*f1*f5-8*f0*f6)*a11*b4+(8*f5*f0+4*f1*f4)*
a8*b7+(-4*f6*f2+f5*f3)*b10*a3+8*f4*f0*a9*b7+8*f4*f0*b9*a7+8*f3*f0*a9*b8+8*f3*f0
*b9*a8+4*f5*f0*a6*b9-2*a5*b11*f0*f5+4*f5*f0*b6*a9-2*b5*a11*f0*f5+4*f1*f6*b7*a6-
6*a3*b11*f1*f6-4*a4*b10*f1*f6-4*b4*a10*f1*f6-6*b3*a11*f1*f6+4*f1*f6*a7*b6-4*a4*
b14*f0*f3+(-16*f5*f0-4*f1*f4)*b4*a12+(-28*f0*f6-2*f1*f5)*a3*b12+(-16*f5*f0-4*f1
*f4)*b12*a4+(-28*f0*f6-2*f1*f5)*b3*a12+

(-8*f1*f0*f6-f1**2*f5)*a11*b14-8*b10*a15*f6*f4*f0-8*a10*b15*f6*f4*f0+(-8*f0*f2*
f6+2*f6*f1**2)*b15*a15+(-4*f1*f0*f6-4*f0*f2*f5+f1**2*f5)*a15*b13+(-4*f1*f0*f6-4
*f0*f2*f5+f1**2*f5)*b15*a13+(-8*f6*f0**2-2*f0*f1*f5)*b15*a14+(-8*f0*f5*f6+f1*f5
**2-8*f1*f4*f6)*b10*a11+(-4*f1*f5*f6-8*f6**2*f0-8*f2*f4*f6+2*f2*f5**2+2*f3**2*
f6)*b10*a10+(-8*f0*f5*f6+f1*f5**2-8*f1*f4*f6)*a10*b11+(-16*f0*f4*f6+2*f0*f5**2-
4*f1*f3*f6)*b11*a11+(-8*f6*f0**2-2*f0*f1*f5)*a15*b14+(-4*f6*f2*f1-2*f1*f3*f5-4*
f5*f4*f0-32*f0*f3*f6)*b11*a12+(-16*f1*f0*f6+2*f1**2*f5-16*f0*f2*f5-8*f4*f3*f0)*
b12*a13+(-4*f6*f2*f1-2*f1*f3*f5-4*f5*f4*f0-32*f0*f3*f6)*b12*a11+(-24*f0*f4*f6-2
*f5*f1*f4-6*f1*f3*f6-4*f0*f5**2)*a12*b10+(-16*f1*f0*f6+2*f1**2*f5-16*f0*f2*f5-8
*f4*f3*f0)*a12*b13+(-4*f0*f3**2-24*f6*f0**2-12*f0*f1*f5-2*f1**2*f4)*a14*b12+(-
24*f0*f4*f6-2*f5*f1*f4-6*f1*f3*f6-4*f0*f5**2)*a10*b12+(-4*f0*f3**2-24*f6*f0**2-
12*f0*f1*f5-2*f1**2*f4)*a12*b14+(-24*f0*f3*f5-8*f4**2*f0+20*f6*f1**2-4*f4*f3*f1
-136*f0*f2*f6-4*f1*f2*f5)*a12*b12-8*b14*a10*f0*f2*f6+(6*f6*f1**2-4*f0*f3*f5-32*
f0*f2*f6)*a15*b12+(6*f6*f1**2-4*f0*f3*f5-32*f0*f2*f6)*b15*a12-8*b11*f0*f3*f6*
a15+(-4*f0*f3*f6-4*f5*f4*f0)*b10*a13+(-8*f0**2*f4-8*f0*f2**2+2*f1**2*f2)*a14*
b14+(-8*f0**2*f5-4*f0*f2*f3-4*f4*f1*f0+f1**2*f3)*a14*b13+(-8*f0**2*f5-4*f0*f2*
f3-4*f4*f1*f0+f1**2*f3)*b14*a13+(-4*f0*f1*f5+2*f1**2*f4-8*f0*f2*f4-8*f6*f0**2)*
b13*a13+(-4*f0*f3*f6-4*f5*f4*f0)*a10*b13+(-2*f6*f1**2-4*f0*f3*f5)*a11*b13-8*a14
*b10*f0*f2*f6+(-8*f1*f0*f6-f1**2*f5)*b11*a14+(-2*f6*f1**2-4*f0*f3*f5)*a13*b11-8
*a11*f0*f3*f6*b15
;

bilinform[4,1] :=
b0*a5+2*a2*b2+a0*b5+

2*a5*b5*f2-4*a6*b1*f6-4*b6*a1*f6-2*f5*b7*a1-2*a3*b4*f5-2*b3*a4*f5-2*f5*a7*b1-6*
a3*b3*f6+

(-4*f4*f0+f1*f3)*a14*b5+(-2*f1*f5-8*f0*f6)*a13*b4+(-2*f1*f5-8*f0*f6)*b13*a4+(-8
*f4*f6+f5**2)*a10*b3+(-8*f4*f6+f5**2)*b10*a3-4*b3*a11*f6*f3-4*a4*b10*f6*f3-8*b4
*a11*f6*f2-8*a4*b11*f6*f2-8*f2*f6*a6*b8-8*f3*f6*b6*a7-4*b4*a10*f6*f3-4*a3*b11*
f6*f3-8*f3*f6*a6*b7-4*f5*f0*a8*b9-8*f2*f6*b6*a8-4*f1*f6*b9*a6-4*f1*f6*a9*b6-4*
b5*a11*f6*f1-4*a5*b11*f6*f1-4*f5*f0*b8*a9+(-f1*f5-4*f0*f6)*a3*b14+(-8*f0*f6-4*
f1*f5)*b8*a8+(-8*f1*f6-4*f2*f5)*b8*a7+(-4*f5*f3-8*f6*f2)*b7*a7+(-2*f1*f5-8*f0*
f6)*a7*b9+(-8*f1*f6-4*f2*f5)*b7*a8+(-2*f1*f5-8*f0*f6)*b7*a9+(2*f5**2-8*f4*f6)*
b6*a6+(-4*f6*f2-2*f5*f3)*a12*b3+(-f1*f5-4*f0*f6)*b3*a14+(-4*f6*f2-2*f5*f3)*b12*
a3+(-4*f4*f0+f1*f3)*b14*a5-2*b3*a13*f6*f1-6*b5*a13*f5*f0-4*b4*a14*f5*f0-2*a3*
b13*f6*f1-6*a5*b13*f5*f0-4*a4*b14*f5*f0-8*b5*f6*f0*a15-4*b4*f6*f1*a15-8*a5*f6*
f0*b15-4*a4*f6*f1*b15+(-16*f1*f6-4*f2*f5)*b12*a4+(-16*f1*f6-4*f2*f5)*b4*a12+(-
28*f0*f6-2*f1*f5)*b5*a12+(-28*f0*f6-2*f1*f5)*b12*a5+

(-8*f0*f5*f6-f1*f5**2)*a10*b13+(-4*f1*f3*f6-2*f0*f5**2)*b13*a11+(-4*f0*f3*f6-4*
f6*f2*f1)*a11*b14+(-4*f1*f5*f6-8*f2*f4*f6-8*f6**2*f0+2*f2*f5**2)*b11*a11+(-4*f2
*f5*f6+f5**2*f3-4*f6*f4*f3-8*f1*f6**2)*a10*b11+(-4*f2*f5*f6+f5**2*f3-4*f6*f4*f3
-8*f1*f6**2)*b10*a11+(-8*f6*f4**2+2*f5**2*f4-8*f6**2*f2)*b10*a10+(-8*f1*f0*f6-8
*f0*f2*f5+f1**2*f5)*a14*b13+(-4*f0*f3*f6-4*f6*f2*f1)*b11*a14+(2*f6*f1**2-4*f0*
f3*f5-16*f0*f2*f6)*b13*a13+(-4*f1*f3*f6-2*f0*f5**2)*a13*b11+(-8*f0*f5*f6-f1*f5
**2)*b10*a13+(-8*f1*f0*f6-8*f0*f2*f5+f1**2*f5)*b14*a13+(2*f1**2*f4-4*f0*f1*f5-8
*f6*f0**2-8*f0*f2*f4+2*f0*f3**2)*b14*a14-8*b10*a14*f6*f4*f0-8*a10*b14*f6*f4*f0-
8*b13*f6*f3*f0*a15-8*b14*f6*f2*f0*a15-8*a13*f6*f3*f0*b15-8*a14*f6*f2*f0*b15+(-8
*f6*f2**2-4*f5*f1*f4+20*f0*f5**2-136*f0*f4*f6-24*f1*f3*f6-4*f3*f5*f2)*a12*b12+(
-8*f6**2*f0-2*f1*f5*f6)*b15*a10+(-8*f6**2*f0-2*f1*f5*f6)*a15*b10+(-6*f0*f3*f5-4
*f6*f1**2-2*f1*f2*f5-24*f0*f2*f6)*a12*b14+(-6*f0*f3*f5-4*f6*f1**2-2*f1*f2*f5-24
*f0*f2*f6)*a14*b12+(-12*f1*f5*f6-24*f6**2*f0-4*f3**2*f6-2*f2*f5**2)*a10*b12+(-
12*f1*f5*f6-24*f6**2*f0-4*f3**2*f6-2*f2*f5**2)*a12*b10+(-4*f6*f2*f1-2*f1*f3*f5-
4*f5*f4*f0-32*f0*f3*f6)*a12*b13+(-4*f6*f2*f1-2*f1*f3*f5-4*f5*f4*f0-32*f0*f3*f6)
*b12*a13+(-16*f0*f5*f6-8*f2*f3*f6+2*f1*f5**2-16*f1*f4*f6)*b11*a12+(-16*f0*f5*f6
-8*f2*f3*f6+2*f1*f5**2-16*f1*f4*f6)*b12*a11+(-4*f1*f3*f6+6*f0*f5**2-32*f0*f4*f6
)*a15*b12+(-4*f1*f3*f6+6*f0*f5**2-32*f0*f4*f6)*b15*a12+(-8*f0*f4*f6+2*f0*f5**2)
*a15*b15+(-4*f1*f4*f6+f1*f5**2-4*f0*f5*f6)*a15*b11+(-4*f1*f4*f6+f1*f5**2-4*f0*
f5*f6)*b15*a11
;

bilinform[4,2] :=
2*a0*b4+2*b1*a2+2*a1*b2+2*b0*a4+

4*f2*a8*b2+4*a1*f4*b7+f3*a12*b0+b0*f1*a14+2*f3*a8*b1+b3*f3*a5+a0*f5*b10+4*a3*f4
*b4+2*a1*f5*b6+4*b1*f4*a7+2*a3*b3*f5+2*f3*a7*b2+a0*f3*b12+f5*a10*b0+f3*a3*b5+4*
f3*a4*b4+4*f2*b8*a2+2*a5*b5*f1+2*f1*b9*a2+2*f3*b8*a1+4*b3*f4*a4+2*a9*b2*f1+2*f3
*b7*a2+4*a5*b4*f2+4*b5*a4*f2+2*b1*f5*a6+a0*f1*b14+

(4*f3*f0+2*f2*f1)*a14*b5+(-4*f0*f6+2*f1*f5)*a11*b5+(4*f6*f2+f5*f3)*a11*b3+(-4*
f0*f6+2*f1*f5)*b11*a5+(2*f4*f5+4*f3*f6)*b10*a3+(4*f6*f2+f5*f3)*b11*a3+(2*f4*f5+
4*f3*f6)*a10*b3+4*f3*f0*a9*b9+4*a4*b11*f5*f2+4*f3*f6*a6*b6+4*b10*f3*f5*a4+4*f3*
f5*a10*b4+2*b5*a10*f5*f2+4*b4*a11*f5*f2+2*a5*b10*f5*f2+(4*f1*f6+4*f2*f5+2*f3*f4
)*b3*a12+(8*f5*f0+4*f1*f4)*b9*a7+(8*f2*f5+4*f3*f4+8*f1*f6)*b7*a7+(8*f4*f0+2*f1*
f3)*b8*a9+(2*f5*f3+8*f6*f2)*b7*a6+(8*f4*f0+2*f1*f3)*b9*a8+(8*f5*f0+8*f1*f4+4*f2
*f3)*b8*a8+(8*f4*f2+8*f1*f5+8*f0*f6)*b8*a7+(8*f0*f6+2*f1*f5)*a9*b6+(4*f2*f5+8*
f1*f6)*b6*a8+(8*f5*f0+4*f1*f4)*b7*a9+(4*f2*f5+8*f1*f6)*a6*b8+(2*f5*f3+8*f6*f2)*
b6*a7+(8*f0*f6+2*f1*f5)*b9*a6+(8*f4*f2+8*f1*f5+8*f0*f6)*a8*b7+(4*f1*f4+4*f5*f0+
2*f2*f3)*b12*a5+(-4*f0*f6+2*f1*f5)*a13*b3+(f1*f3+4*f4*f0)*b5*a13+(-4*f0*f6+2*f1
*f5)*b13*a3+(f1*f3+4*f4*f0)*a5*b13+(4*f3*f0+2*f2*f1)*b14*a5+4*b4*a13*f4*f1+4*b4
*a14*f1*f3+2*b3*f4*f1*a14+4*a4*b13*f4*f1+4*a4*b14*f1*f3+2*a3*f4*f1*b14+2*b5*f0*
f5*a15+2*b4*f1*f5*a15+2*b3*f1*f6*a15+2*a3*f1*f6*b15+2*a4*f1*f5*b15+2*a5*f0*f5*
b15+(4*f1*f6+4*f2*f5+2*f3*f4)*b12*a3+(8*f1*f5+8*f4*f2-8*f0*f6)*b12*a4+(8*f1*f5+
8*f4*f2-8*f0*f6)*b4*a12+(4*f1*f4+4*f5*f0+2*f2*f3)*b5*a12+

(4*f5*f4*f0-12*f0*f3*f6+4*f6*f2*f1+f1*f3*f5)*a13*b11+(4*f0*f5**2+2*f5*f1*f4-8*
f0*f4*f6)*a10*b13+(4*f5*f4*f0-12*f0*f3*f6+4*f6*f2*f1+f1*f3*f5)*a11*b13+(4*f0*f5
*f6+f1*f5**2)*b15*a10+(4*f1*f5*f6+2*f2*f5**2+8*f6**2*f0+2*f3**2*f6)*a10*b11+(8*
f1*f6**2+2*f5**2*f3+4*f6*f4*f3)*b10*a10+(4*f2*f3*f6+8*f0*f5*f6+2*f1*f5**2)*b11*
a11+(4*f1*f5*f6+2*f2*f5**2+8*f6**2*f0+2*f3**2*f6)*b10*a11+(8*f4**2*f0-8*f0*f2*
f6+4*f1*f2*f5+12*f6*f1**2+2*f4*f3*f1)*b12*a13+(32*f6*f2*f1+8*f2**2*f5+8*f1*f4**
2+4*f3*f2*f4-12*f0*f3*f6+32*f5*f4*f0)*a12*b12+2*b10*f3*f5*f1*a14+(4*f0*f5**2+2*
f5*f1*f4-8*f0*f4*f6)*b10*a13+(4*f6*f1**2-8*f0*f2*f6+2*f1*f2*f5)*b11*a14+(2*f1**
2*f5+4*f4*f3*f0+8*f1*f0*f6)*b13*a13+(2*f1**2*f4+8*f6*f0**2+2*f0*f3**2+4*f0*f1*
f5)*a14*b13+(4*f6*f1**2-8*f0*f2*f6+2*f1*f2*f5)*a11*b14+(2*f1**2*f4+8*f6*f0**2+2
*f0*f3**2+4*f0*f1*f5)*b14*a13+(4*f0*f2*f3+2*f1**2*f3+8*f0**2*f5)*a14*b14+2*a10*
f5*f1*f3*b14+4*f0*f3*f6*a15*b15+(4*f0*f5*f6+f1*f5**2)*a15*b10+(12*f1*f0*f6-8*f0
*f2*f5+4*f4*f3*f0+8*f1**2*f5+4*f1*f4*f2+f1*f3**2)*a12*b14+(12*f1*f0*f6-8*f0*f2*
f5+4*f4*f3*f0+8*f1**2*f5+4*f1*f4*f2+f1*f3**2)*a14*b12+(12*f0*f5*f6+8*f1*f5**2+4
*f2*f3*f6+4*f5*f2*f4+f3**2*f5-8*f1*f4*f6)*a10*b12+(12*f0*f5*f6+8*f1*f5**2+4*f2*
f3*f6+4*f5*f2*f4+f3**2*f5-8*f1*f4*f6)*a12*b10+(8*f4**2*f0-8*f0*f2*f6+4*f1*f2*f5
+12*f6*f1**2+2*f4*f3*f1)*a12*b13+(8*f6*f2**2-8*f0*f4*f6+4*f5*f1*f4+12*f0*f5**2+
2*f3*f5*f2)*b11*a12+(8*f6*f2**2-8*f0*f4*f6+4*f5*f1*f4+12*f0*f5**2+2*f3*f5*f2)*
b12*a11+(4*f1*f0*f6+f1**2*f5)*a15*b14+(4*f1*f0*f6+f1**2*f5)*b15*a14+(4*f6*f2*f1
+f1*f3*f5+4*f5*f4*f0+8*f0*f3*f6)*a15*b12+(4*f6*f2*f1+f1*f3*f5+4*f5*f4*f0+8*f0*
f3*f6)*b15*a12+(2*f1*f3*f6+2*f0*f5**2)*a15*b11+(2*f6*f1**2+2*f0*f3*f5)*a15*b13+
(2*f1*f3*f6+2*f0*f5**2)*b15*a11+(2*f6*f1**2+2*f0*f3*f5)*b15*a13
;

bilinform[4,3] :=
a0*b3+2*a1*b1+b0*a3+

2*a3*b3*f4-2*a5*b4*f1-2*b5*a4*f1-4*a9*b2*f0-6*a5*b5*f0-4*b9*a2*f0-2*f1*a8*b2-2*
f1*b8*a2

-8*a4*b13*f0*f4-4*a5*b13*f0*f3-4*b5*a13*f0*f3-4*b3*a13*f0*f5-8*b4*a13*f0*f4+(f1
**2-8*f0*f2)*a14*b5+(f1**2-8*f0*f2)*b14*a5+(-2*f1*f3-4*f4*f0)*a5*b12+(-2*f1*f3-
4*f4*f0)*b5*a12-8*b3*f0*f6*a15-4*b4*a15*f5*f0-4*a4*b15*f5*f0-8*a3*f0*f6*b15-4*
a3*b13*f0*f5-4*b4*a14*f0*f3+(-4*f6*f2+f5*f3)*a10*b3+(-f1*f5-4*f0*f6)*b10*a5+(-2
*f1*f5-8*f0*f6)*b11*a4+(-f1*f5-4*f0*f6)*a10*b5+(-2*f1*f5-8*f0*f6)*a11*b4+(-4*f6
*f2+f5*f3)*b10*a3-8*f4*f0*a9*b7-8*f4*f0*b9*a7-8*f3*f0*a9*b8-8*f3*f0*b9*a8-4*f5*
f0*a6*b9-2*a5*b11*f0*f5-4*f5*f0*b6*a9-2*b5*a11*f0*f5-4*f1*f6*b7*a6-6*a3*b11*f1*
f6-4*a4*b10*f1*f6-4*b4*a10*f1*f6-6*b3*a11*f1*f6-4*f1*f6*a7*b6-4*a4*b14*f0*f3+(-
8*f4*f0-4*f1*f3)*b8*a8+(-2*f1*f5-8*f0*f6)*b6*a8+(-8*f0*f6-4*f1*f5)*b7*a7+(-8*f5
*f0-4*f1*f4)*a8*b7+(-8*f5*f0-4*f1*f4)*b8*a7+(-2*f1*f5-8*f0*f6)*a6*b8+(-8*f0*f2+
2*f1**2)*b9*a9+(-16*f5*f0-4*f1*f4)*b4*a12+(-28*f0*f6-2*f1*f5)*a3*b12+(-16*f5*f0
-4*f1*f4)*b12*a4+(-28*f0*f6-2*f1*f5)*b3*a12+

(-8*f1*f0*f6-f1**2*f5)*a11*b14-8*b10*a15*f6*f4*f0-8*a10*b15*f6*f4*f0+(-8*f0*f2*
f6+2*f6*f1**2)*b15*a15+(-4*f1*f0*f6-4*f0*f2*f5+f1**2*f5)*a15*b13+(-4*f1*f0*f6-4
*f0*f2*f5+f1**2*f5)*b15*a13+(-8*f6*f0**2-2*f0*f1*f5)*b15*a14+(-8*f0*f5*f6+f1*f5
**2-8*f1*f4*f6)*b10*a11+(-4*f1*f5*f6-8*f6**2*f0-8*f2*f4*f6+2*f2*f5**2+2*f3**2*
f6)*b10*a10+(-8*f0*f5*f6+f1*f5**2-8*f1*f4*f6)*a10*b11+(-16*f0*f4*f6+2*f0*f5**2-
4*f1*f3*f6)*b11*a11+(-8*f6*f0**2-2*f0*f1*f5)*a15*b14+(-4*f6*f2*f1-2*f1*f3*f5-4*
f5*f4*f0-32*f0*f3*f6)*b11*a12+(-16*f1*f0*f6+2*f1**2*f5-16*f0*f2*f5-8*f4*f3*f0)*
b12*a13+(-4*f6*f2*f1-2*f1*f3*f5-4*f5*f4*f0-32*f0*f3*f6)*b12*a11+(-24*f0*f4*f6-2
*f5*f1*f4-6*f1*f3*f6-4*f0*f5**2)*a12*b10+(-16*f1*f0*f6+2*f1**2*f5-16*f0*f2*f5-8
*f4*f3*f0)*a12*b13+(-4*f0*f3**2-24*f6*f0**2-12*f0*f1*f5-2*f1**2*f4)*a14*b12+(-
24*f0*f4*f6-2*f5*f1*f4-6*f1*f3*f6-4*f0*f5**2)*a10*b12+(-4*f0*f3**2-24*f6*f0**2-
12*f0*f1*f5-2*f1**2*f4)*a12*b14+(-24*f0*f3*f5-8*f4**2*f0+20*f6*f1**2-4*f4*f3*f1
-136*f0*f2*f6-4*f1*f2*f5)*a12*b12-8*b14*a10*f0*f2*f6+(6*f6*f1**2-4*f0*f3*f5-32*
f0*f2*f6)*a15*b12+(6*f6*f1**2-4*f0*f3*f5-32*f0*f2*f6)*b15*a12-8*b11*f0*f3*f6*
a15+(-4*f0*f3*f6-4*f5*f4*f0)*b10*a13+(-8*f0**2*f4-8*f0*f2**2+2*f1**2*f2)*a14*
b14+(-8*f0**2*f5-4*f0*f2*f3-4*f4*f1*f0+f1**2*f3)*a14*b13+(-8*f0**2*f5-4*f0*f2*
f3-4*f4*f1*f0+f1**2*f3)*b14*a13+(-4*f0*f1*f5+2*f1**2*f4-8*f0*f2*f4-8*f6*f0**2)*
b13*a13+(-4*f0*f3*f6-4*f5*f4*f0)*a10*b13+(-2*f6*f1**2-4*f0*f3*f5)*a11*b13-8*a14
*b10*f0*f2*f6+(-8*f1*f0*f6-f1**2*f5)*b11*a14+(-2*f6*f1**2-4*f0*f3*f5)*a13*b11-8
*a11*f0*f3*f6*b15
;

bilinform[4,4] :=
a0*b0+

0+

(8*f4*f0-2*f1*f3)*a5*b5+8*a5*f5*f0*b4+8*b5*f5*f0*a4+8*b3*f6*f1*a4+8*a3*f6*f1*b4
+8*a3*b5*f6*f0+8*a5*b3*f6*f0+32*f6*f0*a4*b4+(8*f6*f2-2*f5*f3)*a3*b3+

(16*f0*f2*f5-4*f1**2*f5)*a5*b13+(8*f0*f3*f5-8*f6*f1**2+32*f0*f2*f6)*b4*a13+16*
f6*f3*f0*a4*b15+16*f6*f3*f0*b4*a15+(4*f1*f3*f5+96*f0*f3*f6)*b12*a4+(-16*f6*f1**
2+64*f0*f2*f6+12*f0*f3*f5)*b5*a12+(-16*f6*f1**2+64*f0*f2*f6+12*f0*f3*f5)*b12*a5
+(-4*f1*f5**2+16*f0*f5*f6+16*f1*f4*f6)*b4*a10+(-4*f1*f5**2+16*f0*f5*f6+16*f1*f4
*f6)*b10*a4+(4*f1*f5*f6-4*f3**2*f6-4*f2*f5**2+16*f2*f4*f6)*b10*a3+16*a14*b3*f6*
f2*f0+16*a13*b3*f6*f3*f0+16*a3*b13*f6*f3*f0+16*a3*b14*f6*f2*f0+(16*f0*f2*f5-4*
f1**2*f5)*a13*b5+(-4*f1**2*f5+16*f1*f0*f6+16*f0*f2*f5)*b4*a14+(-4*f1**2*f4+4*f0
*f1*f5-4*f0*f3**2+16*f0*f2*f4)*b5*a14+(8*f0*f3*f5-8*f6*f1**2+32*f0*f2*f6)*a4*
b13+(-4*f1**2*f4+4*f0*f1*f5-4*f0*f3**2+16*f0*f2*f4)*b14*a5+(-4*f1**2*f5+16*f1*
f0*f6+16*f0*f2*f5)*b14*a4+(16*f1*f4*f6-4*f1*f5**2)*a11*b3+(16*f1*f4*f6-4*f1*f5
**2)*a3*b11+(8*f1*f3*f6+32*f0*f4*f6-8*f0*f5**2)*b4*a11+(4*f1*f5*f6-4*f3**2*f6-4
*f2*f5**2+16*f2*f4*f6)*b3*a10+(8*f1*f3*f6+32*f0*f4*f6-8*f0*f5**2)*a4*b11+16*a10
*b5*f6*f4*f0+16*a5*b10*f6*f4*f0+16*a5*b11*f6*f3*f0+16*a11*b5*f6*f3*f0+(-4*f6*f1
**2+16*f0*f2*f6)*a15*b5+(-4*f6*f1**2+16*f0*f2*f6)*b15*a5+(12*f1*f3*f6+64*f0*f4*
f6-16*f0*f5**2)*b3*a12+(4*f1*f3*f5+96*f0*f3*f6)*b4*a12+(-4*f0*f5**2+16*f0*f4*f6
)*a15*b3+(-4*f0*f5**2+16*f0*f4*f6)*b15*a3+(12*f1*f3*f6+64*f0*f4*f6-16*f0*f5**2)
*a3*b12+

(-4*f5*f2*f1**2+32*f6*f2*f1*f0-16*f6*f3*f0**2-8*f6*f1**3+16*f5*f2**2*f0-4*f5*f3
*f1*f0+16*f5*f4*f0**2)*b13*a14+(80*f6*f5*f0**2+8*f5*f3**2*f0+96*f6*f3*f2*f0-20*
f6*f3*f1**2-16*f6*f4*f1*f0+8*f5**2*f1*f0)*b12*a13+(-2*f5*f3*f1**2+24*f6*f3*f1*
f0+16*f5*f3*f2*f0-16*f6*f2*f1**2+64*f6*f2**2*f0+16*f5**2*f0**2+32*f6*f4*f0**2)*
a14*b12+(8*f6*f5*f1*f0+16*f6*f3**2*f0+32*f6**2*f0**2)*a13*b11+(16*f6**2*f1*f0+
16*f6*f5*f2*f0+16*f6*f4*f3*f0-4*f6*f5*f1**2)*b10*a13+(16*f6*f3*f2*f0+16*f6*f4*
f1*f0-4*f5**2*f1*f0+16*f6*f5*f0**2)*b11*a14+(8*f6*f5*f1*f0+16*f6*f3**2*f0+32*f6
**2*f0**2)*a11*b13+(16*f6**2*f1**2+32*f6*f4**2*f0-8*f6*f5*f3*f0+8*f6*f4*f3*f1-8
*f5**2*f4*f0-2*f5**2*f3*f1)*a11*b11+(-8*f5**3*f0-4*f5**2*f4*f1-4*f6*f5*f3*f1-16
*f6**2*f3*f0+32*f6*f5*f4*f0+16*f6**2*f2*f1+16*f6*f4**2*f1)*b10*a11+(-8*f5**3*f0
-4*f5**2*f4*f1-4*f6*f5*f3*f1-16*f6**2*f3*f0+32*f6*f5*f4*f0+16*f6**2*f2*f1+16*f6
*f4**2*f1)*b11*a10+(-4*f6*f4*f3**2+16*f6**2*f2**2-16*f6**2*f3*f1+4*f6*f5**2*f0-
6*f5**3*f1+24*f6*f5*f4*f1-8*f6*f5*f3*f2+16*f6**2*f4*f0+16*f6*f4**2*f2-4*f5**2*
f4*f2+f5**2*f3**2)*b10*a10+(16*f6**2*f1*f0+16*f6*f5*f2*f0+16*f6*f4*f3*f0-4*f6*
f5*f1**2)*b13*a10+(16*f5**2*f0**2-2*f5*f3*f1**2-8*f6*f2*f1**2+8*f5*f3*f2*f0-8*
f6*f3*f1*f0+32*f6*f2**2*f0)*a13*b13+(-16*f5*f3*f0**2+f3**2*f1**2+4*f6*f1**2*f0+
16*f4*f2**2*f0-4*f3**2*f2*f0+16*f4**2*f0**2-8*f4*f3*f1*f0-4*f4*f2*f1**2-6*f5*f1
**3+24*f5*f2*f1*f0+16*f6*f2*f0**2)*b14*a14+(16*f6*f3*f2*f0+16*f6*f4*f1*f0-4*f5
**2*f1*f0+16*f6*f5*f0**2)*b14*a11+(8*f6*f5*f1*f0-4*f6*f3**2*f0-2*f5**2*f1**2+16
*f6**2*f0**2+32*f6*f4*f2*f0)*b14*a10+(-4*f5*f2*f1**2+32*f6*f2*f1*f0-16*f6*f3*f0
**2-8*f6*f1**3+16*f5*f2**2*f0-4*f5*f3*f1*f0+16*f5*f4*f0**2)*b14*a13+(8*f6*f5*f1
*f0-4*f6*f3**2*f0-2*f5**2*f1**2+16*f6**2*f0**2+32*f6*f4*f2*f0)*b10*a14+(80*f6*
f5*f0**2+8*f5*f3**2*f0+96*f6*f3*f2*f0-20*f6*f3*f1**2-16*f6*f4*f1*f0+8*f5**2*f1*
f0)*a12*b13+(-4*f6*f2*f1**2+16*f6*f2**2*f0+16*f6*f4*f0**2)*b15*a14+(-16*f6*f5*
f2*f0+8*f6*f3**2*f1+96*f6*f4*f3*f0+8*f6*f5*f1**2-20*f5**2*f3*f0+80*f6**2*f1*f0)
*b11*a12+(-4*f6*f2*f1**2+16*f6*f2**2*f0+16*f6*f4*f0**2)*a15*b14+(64*f6**2*f0**2
+4*f5**2*f1**2+8*f6*f5*f1*f0+64*f6*f4*f2*f0-16*f5**2*f2*f0-16*f6*f4*f1**2+16*f6
*f3**2*f0)*a15*b12+(64*f6**2*f0**2+4*f5**2*f1**2+8*f6*f5*f1*f0+64*f6*f4*f2*f0-
16*f5**2*f2*f0-16*f6*f4*f1**2+16*f6*f3**2*f0)*b15*a12+(-2*f5*f3*f1**2+24*f6*f3*
f1*f0+16*f5*f3*f2*f0-16*f6*f2*f1**2+64*f6*f2**2*f0+16*f5**2*f0**2+32*f6*f4*f0**
2)*a12*b14+(16*f6**2*f1**2-16*f5**2*f4*f0+24*f6*f5*f3*f0-2*f5**2*f3*f1+32*f6**2
*f2*f0+16*f6*f4*f3*f1+64*f6*f4**2*f0)*a10*b12+(16*f6**2*f1**2-16*f5**2*f4*f0+24
*f6*f5*f3*f0-2*f5**2*f3*f1+32*f6**2*f2*f0+16*f6*f4*f3*f1+64*f6*f4**2*f0)*a12*
b10+(16*f6*f5*f0**2-4*f6*f3*f1**2+16*f6*f3*f2*f0)*b15*a13+(-16*f6*f5*f2*f0+8*f6
*f3**2*f1+96*f6*f4*f3*f0+8*f6*f5*f1**2-20*f5**2*f3*f0+80*f6**2*f1*f0)*b12*a11+(
-4*f5**2*f4*f0+16*f6*f4**2*f0+16*f6**2*f2*f0)*b15*a10+(16*f6*f4*f3*f0+16*f6**2*
f1*f0-4*f5**2*f3*f0)*b15*a11+(16*f6*f4*f3*f0+16*f6**2*f1*f0-4*f5**2*f3*f0)*a15*
b11+(-4*f5**2*f4*f0+16*f6*f4**2*f0+16*f6**2*f2*f0)*a15*b10+(-4*f6*f4*f1**2+f5**
2*f1**2+16*f6*f4*f2*f0-4*f5**2*f2*f0+16*f6**2*f0**2)*b15*a15+(16*f6*f5*f0**2-4*
f6*f3*f1**2+16*f6*f3*f2*f0)*a15*b13+(288*f6**2*f0**2+144*f6*f3**2*f0+64*f6*f5*
f1*f0-80*f6*f4*f1**2+8*f5*f4*f3*f0+24*f5**2*f1**2+4*f3**2*f5*f1-80*f5**2*f2*f0+
8*f6*f3*f2*f1+288*f6*f4*f2*f0)*a12*b12:

##############################################################################
##
## sumkummer
##
## given two points in coord-form, calculates kummer of sum. Point returned
## might be [0,0,0,0]. idx can be chosen from [1..4] and can be used to select
## different projective representations.
##
##############################################################################

sumkummer:=proc(idx,av,bv)
  local vals,sbs,i;
  vals:=[bilinform[idx,i]$i=1..4];
  sbs:={seq(a.i=av[i+1],i=0..15),seq(b.i=bv[i+1],i=0..15)};
  RETURN(subs(sbs,vals));
end;

##############################################################################
##
## kummer2theta
##
## sends a point on kummer to theta (=0 if point is on diagonal)
##
##############################################################################

kummer2theta:=k->k[2]^2-4*k[1]*k[3];

##############################################################################
##
## hpol
##
## returns the polynomial to check for way galois acts. offset can be used to
## get rid of square factors.
##
##############################################################################

hpol:=proc(offset)
  local sextic,i,t0,t1,t2,t3,t4,t5,t6,k0,k1,k2,k3,k4,k5,
        h0,h1,h2,h3,h4,h5,h6,h7,h8,h9;
  sextic:=collect(sum('f.i'*(x+offset)^i,i=0..6),x);

  t6 := coeff(sextic,x,6):
  t5 := coeff(sextic,x,5):
  t4 := coeff(sextic,x,4):
  t3 := coeff(sextic,x,3):
  t2 := coeff(sextic,x,2):
  t1 := coeff(sextic,x,1):
  t0 := coeff(sextic,x,0):

  k5 := t5/t6:
  k4 := t4/t6:
  k3 := t3/t6:
  k2 := t2/t6:
  k1 := t1/t6:
  k0 := t0/t6:

  h0 := k0^4*k5^6+k0^3*k2^2*k5^4-2*k0^3*k1*k3*k5^4-4*k0^4*k4*k5^4+k0^2*k1^2*k3^2*
  k5^2-2*k0^3*k2*k3^2*k5^2-2*k0^2*k1^2*k2*k4*k5^2+8*k0^3*k1*k3*k4*k5^2+2*k0^2*k1^
  3*k5^3-4*k0^3*k1*k2*k5^3+8*k0^4*k3*k5^3+k0^3*k3^4-2*k0^2*k1^2*k3^2*k4+k0*k1^4*
  k4^2-2*k0*k1^4*k3*k5+8*k0^2*k1^2*k2*k3*k5-12*k0^3*k1*k3^2*k5-4*k0^2*k1^3*k4*k5+
  k1^6-4*k0*k1^4*k2+8*k0^2*k1^3*k3:

  h1 := k1^5*k4+8*k0^2*k1^3-8*k0^3*k3^3+8*k0^4*k5^3+k0^3*k2*k5^5-3*k0^3*k3^2*k5^3
  -3*k0*k1^3*k3^2+2*k0^3*k1*k5^4+2*k0*k1^4*k5+k0^2*k1*k3^3*k5+k0^2*k1*k2*k3*k5^3+
  k0*k1^3*k3*k4*k5+8*k0^2*k1*k2*k3^2+24*k0^3*k1*k3*k5-8*k0^2*k1^2*k3*k4-3*k0*k1^3
  *k2*k5^2-16*k0^2*k1^2*k2*k5+6*k0^2*k1^2*k3*k5^2+8*k0^2*k1^2*k4^2*k5-8*k0^3*k2*
  k3*k5^2-16*k0^3*k1*k4*k5^2-8*k0^2*k1*k2*k3*k4*k5+8*k0^3*k3^2*k4*k5-3*k0^2*k1^2*
  k4*k5^3+8*k0^2*k1*k2^2*k5^2:

  h2 := -k1^4*k2+18*k0^3*k3^2+k1^4*k3*k5-k0^3*k4*k5^4-7*k0^2*k1^2*k5^2-4*k0*k1^3*
  k3-6*k0^2*k1^2*k4+8*k0^3*k4^2*k5^2+8*k0*k1^2*k2^2-6*k0^3*k2*k5^2-4*k0^3*k3*k5^3
  +36*k0^3*k1*k5-4*k0*k1^3*k5^3+k0^2*k2*k3^2*k5^2+k0*k1^2*k2*k4*k5^2+k0*k1^2*k3^2
  *k4+k0^2*k1*k3*k5^4-24*k0^3*k3*k4*k5-4*k0^2*k1*k2*k4*k5-8*k0*k1^2*k2*k3*k5+10*
  k0*k1^3*k4*k5-24*k0^2*k1*k2*k3+8*k0^2*k1*k3*k4^2+8*k0^2*k2^2*k3*k5+8*k0^2*k1*k3
  ^2*k5-2*k0*k1^2*k2*k4^2+10*k0^2*k1*k2*k5^3-2*k0^2*k2*k3^2*k4-2*k0^2*k2^2*k4*k5^
  2-8*k0^2*k1*k3*k4*k5^2:

  h3 := k0^2*k3^2*k5^3+k0*k1^2*k4*k5^3-2*k0^2*k2*k4*k5^3-k0^2*k1*k5^4+k0*k1*k2*k3
  *k4*k5-3*k0^2*k3^2*k4*k5-3*k0*k1^2*k4^2*k5+k1^3*k2*k5^2-3*k0*k1*k2^2*k5^2-6*k0*
  k1^2*k3*k5^2+11*k0^2*k2*k3*k5^2+14*k0^2*k1*k4*k5^2-9*k0^3*k5^3+k1^3*k3^2-3*k0*
  k1*k2*k3^2+3*k0^2*k3^3-2*k1^3*k2*k4+11*k0*k1^2*k3*k4-k1^4*k5+14*k0*k1^2*k2*k5-
  39*k0^2*k1*k3*k5-9*k0*k1^3:

  h4 := k0*k1*k3*k4*k5^2-k0^2*k4^2*k5^2+k1^3*k5^3-3*k0*k1*k2*k5^3+k0^2*k3*k5^3+k0
  *k2^2*k4^2-2*k0*k1*k3*k4^2-4*k0^2*k4^3+k1^2*k2*k3*k5-2*k0*k2^2*k3*k5-k0*k1*k3^2
  *k5-3*k1^3*k4*k5+16*k0^2*k3*k4*k5+6*k0*k1^2*k5^2-2*k0^2*k2*k5^2-k1^2*k2^2-4*k0*
  k2^3+k1^3*k3+16*k0*k1*k2*k3-15*k0^2*k3^2-2*k0*k1^2*k4+18*k0^2*k2*k4-30*k0^2*k1*
  k5-27*k0^3:

  h5 := k0*k2*k4^2*k5+k1^2*k3*k5^2-2*k0*k2*k3*k5^2-2*k0*k1*k4*k5^2+2*k0^2*k5^3+k1
  *k2^2*k4-2*k1^2*k3*k4-k0*k2*k3*k4-4*k0*k1*k4^2-2*k1^2*k2*k5-4*k0*k2^2*k5+16*k0*
  k1*k3*k5+3*k0^2*k4*k5+2*k1^3+3*k0*k1*k2+9*k0^2*k3:

  h6 := k0*k4^3+k1*k2*k4*k5-3*k0*k3*k4*k5-k1^2*k5^2+k0*k2*k5^2+k2^3-3*k1*k2*k3+3*
  k0*k3^2+k1^2*k4-9*k0*k2*k4+9*k0*k1*k5+27*k0^2:

  h7 := k1*k4^2+k2^2*k5-2*k1*k3*k5-k0*k4*k5-k1*k2-6*k0*k3:

  h8 := k2*k4-k1*k5-9*k0:

  h9 := k3:
 
  RETURN(expand(t6**10*(x**10 + h9*x**9 + h8*x**8 + h7*x**7 + h6*x**6 + h5*x**5 + h4*x**4
         + h3*x**3 + h2*x**2 + h1*x + h0))):
end: 
#--- flynn.mpl (end) ---#
#--- divcalc.mpl (start) ---#
cat << "#--- divcalc.mpl (end) ---#" > divcalc.mpl # UNIX extraction command
##############################################################################
##############################################################################
##
## divcalc.mpl
##
## for Maple Vr3 or Vr4
##
## version 1.0, January 28 1997
##
## Nils Bruin, nbruin@wi.leidenuniv.nl
##
## Maple-routines for doing calculations on the jacobian of a genus 2 curve.
##
## Used representation:
## []: Principal divisor class
## [P,Q]: Class represented by the divisor {P+Q-(infty^+)-(infty^-)]
## Points at infinity are represented by [infinity,Y/X^3] i.e.
##   [infinity, +/- sqrt(f6)].
## Quadratic points are represented by RootOf's.
##
## adddiv(F,D1,D2) calculates D1+D2 op y^2=F(x)
## adddivp(p,F,D1,D2) calculates D1+D2 mod p (p good)
## muldiv(F,n,D) calculates n*D
## muldivp(p,F,,n,D) calculates n*D mod p.
##
## Global variables: x,y,g0,g1,g2,g3
##   these should remain unassigned.
##
## (C) This code is regarded as a scientific publication. Work (partially)
## based on results obtained using this code should contain an appropriate
## reference. You are free to copy this file unaltered for non-commercial
## purposes.
## This code is provided as is, with no claim of fitness for any particular
## purpose. Although great care has been taken to eliminate all bugs, there is
## no guarantee that results produced are correct.
##############################################################################
##############################################################################

############################################################################################################################################################
##
## sortdiv
##
## Sorts a divisor, for easier handling. Also does some trivial simplification.
## The result still represents the same divisor class.
##
## IN: list of points (an effective divisor)
## OUT: List sorted on x, with multiplicity. [x,y]+[x,-y]=0 is applied.
## 
## Returned value is of the form
## [ [[x1,y1],e1], ..., [[xr,yr],er] ]
## The ei are positive.
##############################################################################

sortdiv:=proc(d)
  local res,xhad,ref,e,div;
  res:=NULL;
  xhad:=NULL;
  for ref in d do
    if not(member(ref[1],{xhad})) then
      xhad:=xhad,ref[1];
      e:=0;
      for div in d do
        if div[1]=ref[1] then
          if ref[2]=0 then
            e:=1-e;
          elif div[2]+ref[2]=0 then
            e:=e-1;
          else
            e:=e+1;
          fi;
        fi;
      od;
      if e>0 then
        res:=res,[ref,e];
      elif e<0 then
        res:=res,[[ref[1],-ref[2]],-e];
      fi;
    fi;
  od;
  RETURN([res]);
end:

##############################################################################
##
## genrel
##
## returns a comma-separated sequence of relations on g0..g3 (global variables)
## produceert relaties op g (=g0+g1*x+g2*x^2+g3*x^3) zodat y-g(x) een e-de orde nulpunt
## op y^2=F(x) heeft in (x0,y0).
##
## Als y0=0, dan moet e=1.
##
##############################################################################

genrel:=proc(Fin,Gin,x0in,y0in,e)
  local F1,F2,F3,G1,G2,G3,rels,F,G,x0,y0;
  F:=Fin;G:=Gin;x0:=x0in;y0:=y0in;
  if e>1 and y0=0 then
    ERROR(`2-division point with high multiplicity`);
  fi;
  if x0=infinity then
    # als x0=infty, ga over op z=1/x; w=y/x^3 ("y0" is al w in dat geval)
    F:=collect(x^6*subs(x=1/x,F),x);
    G:=collect(x^3*subs(x=1/x,G),x);
    x0:=0;
  fi;
  rels:=NULL;
  F1:=diff(F,x);
  F2:=diff(F1,x);
  F3:=diff(F2,x);
  G1:=diff(G,x);
  G2:=diff(G1,x);
  G3:=diff(G2,x);
  
  if e>=1 then
    rels:=rels,subs(x=x0,y=y0,y-G);
  fi;
  
  # als (x0,y0) geen weierstrass-punt, dan is x-x0 een uniformizer. Differentieer in
  # die gevallen y-g(x) voldoende malen naar x & stel 0.
  
  if e>=2 then
    rels:=rels,subs(x=x0,y=y0,F1-2*y*G1);
  fi;
  if e>=3 then
    rels:=rels,subs(x=x0,y=y0,2*F2*y^2-F1^2-4*G2*y^3);
  fi;
  if e=4 then
    rels:=rels,subs(x=x0,y=y0,F2*F1+2*F*F3-6*G2*F1*y-2*y*G1*F2-4*G3*F*y);
  fi;
  RETURN(rels);
end:

##############################################################################
##
## adddiv
##
## IN: 2 divisoren div1 en div2 (effectief, graad <=2)
## UIT: Een divisor div3 (efectief, graad <=2) zodat div1+div2=div3 in
## pic0-zin.
##
##############################################################################

adddiv:=proc(F,div1,div2)
  option remember;
  local div,relations,pnt,xi,yi,ei,x1,x2,y1,y2,sume,knowninfty,gx,H;
  div:=sortdiv([op(div1),op(div2)]);
  relations:=NULL;
  sume:=0;
  knowninfty:=0;
  for pnt in div do
    xi:=pnt[1][1];
    yi:=pnt[1][2];
    ei:=pnt[2];
    sume:=sume+ei;
    if xi=infinity then knowninfty:=knowninfty+ei; fi;
    relations:=relations,genrel(F,g0+g1*x+g2*x^2+g3*x^3,xi,yi,ei);
  od;
  if sume<=2 then
    RETURN([seq(pnt[1]$pnt[2],pnt=div)]);
  fi;
  
  gx:=subs(solve({relations}),g0+g1*x+g2*x^2+g3*x^3);
  H:=F-gx^2;
  for pnt in div do
    if pnt[1][1]<>infinity then
      H:=simplify(H/(x-pnt[1][1])^pnt[2]);
    fi; 
  od;   
  if degree(H,x)=0 then
    if H=0 then
      RETURN([`strange result`]);
    else
      RETURN([[infinity,-coeff(gx,x,3)]$2]);
    fi;
  elif degree(H,x)=1 then
    xi:=-coeff(H,x,0)/coeff(H,x,1);yi:=subs(x=xi,gx);
    RETURN([[xi,-yi],[infinity,-coeff(gx,x,3)]]);
  elif irreduc(H) then
    x1:=RootOf(H); y1:=evala(subs(x=x1,gx));
    x2:=evala(evala(Trace(x1))-x1);y2:=evala(subs(x=x2,gx));
    RETURN([[x1,-y1],[x2,-y2]]);
  else
    x1:=solve(H);
    if nops([x1])=2 then
      x2:=op(2,[x1]);x1:=op(1,[x1]);
    else
      x2:=x1;
    fi;
    y1:=subs(x=x1,gx);y2:=subs(x=x2,gx);
    RETURN([[x1,-y1],[x2,-y2]]);
  fi;
end:
 
##############################################################################
##
## sortdivp
##
## IN: lijst punten ('n divisor)
## UIT: de lijst gesorteerd op x, met multipliciteit. De vereenvoudiging
## [x,y]+[x,-y]=0 wordt gebruikt.
## Output-format is dus:
## [ [[x1,y1],e1], ..., [[xr,yr],er] ]
## De ei zijn alle positief.
##
##############################################################################

sortdivp:=proc(p,d)
  local res,xhad,ref,e,div;
  res:=NULL;
  xhad:=NULL;
  for ref in d do
    if ref[1]=infinity then
      ref:=[ref[1], ref[2] mod p];
    else
      ref:=[ref[1] mod p, ref[2] mod p];
    fi;
    if not(member(ref[1],{xhad})) then
      xhad:=xhad,ref[1];
      e:=0;
      for div in d do
        if (div[1]=infinity and ref[1]=infinity)or(div[1]-ref[1])mod p =0 then
          if ref[2]=0 then
            e:=1-e;
          elif (div[2]+ref[2]) mod p=0 then
            e:=e-1;
          else
            e:=e+1;
          fi;
        fi;
      od;
      if e>0 then
        res:=res,[ref,e];
      elif e<0 then
        res:=res,[[ref[1],(-ref[2])mod p],-e];
      fi;
    fi;
  od;
  RETURN([res]);
end:

##############################################################################
##
## invdiv
##
## inverts a divisor
##############################################################################

invdiv:=div->map(u->[u[1],-u[2]],div):

##############################################################################
##
## adddivp
##
## IN: 2 divisoren div1 en div2 (effectief, graad <=2)
## UIT: Een divisor div3 (efectief, graad <=2) zodat div1+div2=div3 in
## pic0-zin.
##
##############################################################################

adddivp:=proc(p,F,div1,div2)
  local div,relations,pnt,xi,yi,ei,x1,x2,y1,y2,sume,knowninfty,gx,H,known;
  div:=sortdivp(p,[op(div1),op(div2)]);
  relations:=NULL;
  sume:=0;
  knowninfty:=0;
  for pnt in div do
    xi:=pnt[1][1];
    yi:=pnt[1][2];
    ei:=pnt[2];
    sume:=sume+ei;
    if xi=infinity then knowninfty:=knowninfty+ei; fi;
    relations:=relations,genrel(F,g0+g1*x+g2*x^2+g3*x^3,xi,yi,ei);
  od;
  if sume<=2 then
    RETURN([seq(pnt[1]$pnt[2],pnt=div)]);
  fi;
  
  gx:=subs(op(1,[solve({relations})])mod p,(g0+g1*x+g2*x^2+g3*x^3));
  gx:=Normal(evala(gx))mod p;
  H:=(F-gx^2)mod p;
  known:=1;
  for pnt in div do
    if pnt[1][1]<>infinity then
      known:=Expand(known*(x-pnt[1][1])^pnt[2]) mod p;
    fi; 
  od;
  H:=Normal(H/known) mod p;   
  if degree(H,x)=0 then
    if H=0 then
      RETURN([`strange result`]);
    else
      RETURN([[infinity,(-coeff(gx,x,3))mod p]$2]);
    fi;
  elif degree(H,x)=1 then
    xi:=(-coeff(H,x,0)/coeff(H,x,1))mod p;yi:=(subs(x=xi,gx))mod p;
    RETURN([[xi,(-yi)mod p],[infinity,(-coeff(gx,x,3))mod p]]);
  elif Irreduc(H)mod p then
    x1:=RootOf(H)mod p; y1:=Normal(subs(x=x1,gx))mod p;
    x2:=Normal(evala(Trace(x1))-x1)mod p;y2:=Normal(subs(x=x2,gx))mod p;
    RETURN([[x1,(-y1)mod p],[x2,(-y2)mod p]]);
  else
    x1:=map(u->u[1]$u[2],Roots(H)mod p);
    x2:=x1[2];x1:=x1[1];
    y1:=subs(x=x1,gx)mod p;y2:=subs(x=x2,gx)mod p;
    RETURN([[x1,(-y1)mod p],[x2,(-y2)mod p]]);
  fi;
end:

##############################################################################
##
## muldiv
##
## returns the multiple of a divisor
##############################################################################

muldiv:=proc(F,n,div)
  local div2;
  option remember;
  if n<0 then
    RETURN(invdiv(muldiv(F,-n,div)));
  elif n=0 then
    RETURN([]);
  elif n=1 then
    RETURN(div);
  elif n mod 2 = 0 then
    div2:=muldiv(F,n/2,div);
    RETURN(adddiv(F,div2,div2));
  else
    div2:=muldiv(F,(n-1)/2,div);
    div2:=adddiv(F,div2,div2);
    RETURN(adddiv(F,div,div2));
  fi;
end:

##############################################################################
##
## muldivp
##
## returns the multiple of a divisor mod p
##############################################################################

muldivp:=proc(p,F,n,div)
  local div2;
  option remember;
  if n<0 then
    RETURN(invdiv(muldivp(p,F,-n,div))mod p);
  elif n=0 then
    RETURN([]);
  elif n=1 then
    RETURN(div);
  elif n mod 2 = 0 then
    div2:=muldivp(p,F,n/2,div);
    RETURN(adddivp(p,F,div2,div2));
  else
    div2:=muldivp(p,F,(n-1)/2,div);
    div2:=adddivp(p,F,div2,div2);
    RETURN(adddivp(p,F,div,div2));
  fi;
end:
  
##############################################################################
##
## initcurve
##
## initializes routines for effective chabauty-computations
## requires availability of "flynn.mpl"
##############################################################################

initcurve:=proc(F)
  global f0,f1,f2,f3,f4,f5,f6,ad,adp,ml,mlp;
  local i;
  for i from 0 to 6 do f.i:=coeff(F,x,i); od;
  read `flynn.mpl`;
  ad:=subs(Ft=F,(A,B)->adddiv(Ft,A,B));
  adp:=subs(Ft=F,(p,A,B)->adddivp(p,Ft,A,B));
  ml:=subs(Ft=F,(n,A)->muldiv(Ft,n,A));
  mlp:=subs(Ft=F,(p,n,A)->muldivp(p,Ft,n,A));
  print(`current curve is `,y^2=F);
end:

##############################################################################
##
## countcrvp
##
## counts number of points of y^2=F(x) mod p
##############################################################################

countcrvp:=proc(p,F)
  local f6,N,xt,y,f;
  option remember;
  f6:=coeff(F,x,6) mod p;
  N:=0;
  if f6=0 then
    N:=1;
  elif Irreduc(y^2-f6)mod p then
    N:=0;
  else
    N:=2;
  fi;
  for xt from 0 to p-1 do
    f:=subs(x=xt,F)mod p;
    if f=0 then
      N:=N+1;
    elif not(Irreduc(y^2-f)mod p) then
      N:=N+2;
    fi;
  od;
  RETURN(N);
end:

##############################################################################
##
## countcrvpsqr
##
## counts number of points of y^2=F(x) over F_(p^2)
##############################################################################

countcrvpsqr:=proc(p,F)
  local alpha,f6,N,xt1,xt2,f;
  option remember;
  alpha:=RootOf(Randprime(2,_Z) mod p);
  f6:=coeff(F,x,6) mod p;
  N:=0;
  if f6=0 then
    N:=1;
  elif Irreduc(y^2-f6,alpha)mod p then
    N:=0;
  else
    N:=2;
  fi;
  for xt1 from 0 to p-1 do
    for xt2 from 0 to p-1 do
      f:=Normal(subs(x=xt1*alpha+xt2,F))mod p;
      if f=0 then
        N:=N+1;
      elif not(Irreduc(y^2-f,alpha)mod p) then
        N:=N+2;
      fi;
    od;
  od;
  RETURN(N);
end:

##############################################################################
##
## njac
##
## calculates number of points of jac(y^2-F(x))
##############################################################################

njac:=(p,F)->countcrvpsqr(p,F)/2+countcrvp(p,F)^2/2-p:

##############################################################################
##
## quadseek
##
## search for quadratic points on y^2=F(x)
##############################################################################

readlib(issqr):

quadseek:=proc(F,Ab,Bb,Cb)
  local Al,Bl,Cl,Au,Bu,Cu,a,b,c,alpha,Y,stp,Result;

  Result:=NULL;

  if nargs=5 then
    stp:=args[5];
  else
    stp:=true;
  fi;

  if type(Ab,`..`) then
    Al:=op(1,Ab);Au:=op(2,Ab);
  else
    Al:=1;Au:=Ab;
  fi;
  if type(Bb,`..`) then
    Bl:=op(1,Bb);Bu:=op(2,Bb);
  else
    Bl:=-Bb;Bu:=Bb;
  fi;
  if type(Cb,`..`) then
    Cl:=op(1,Cb);Cu:=op(2,Cb);
  else
    Cl:=-Cb;Cu:=Cb;
  fi;
  
  for a from Al to Au do 
    for b from Bl to Bu do
      for c from Cl to Cu do
        if not(issqr(b^2-4*a*c)) then
          alpha:=RootOf(a*_Z^2+b*_Z+c);
          if not(evala(irreduc(Y^2-subs(x=alpha,F)),alpha)) then
            if stp then
              RETURN([alpha,evala(Roots(Y^2-subs(x=alpha,F)),alpha)[1][1]]);
            else
              Result:=Result,
                 [alpha,evala(Roots(Y^2-subs(x=alpha,F)),alpha)[1][1]];
            fi;
          fi;
        fi;
      od;
    od;
  od;
  RETURN([Result]);
end:

##############################################################################
##
## maketablep
##
## produces a table of t1+nG...tm+nG for n=1...ord_p(G)
##############################################################################

maketablep:=proc(p,F,Tlist,G)
  local Result,n,nG,tmp,T;

  Result:=[0,op(Tlist)];
  n:=1;
  nG:=G;
  while nG<>[] do
    tmp:=[n,seq(adddivp(p,F,nG,T),T=Tlist)];
    #print(tmp);
    Result:=Result,tmp;
    n:=n+1;
    nG:=adddivp(p,F,nG,G);
  od;
  RETURN(array(1..nops([Result]),1..1+nops(Tlist),[Result]));
end:

##############################################################################
##
## conj
##
## conjugates a list of quadratic points
##############################################################################

conj:=proc(L)
  local f;
  f:=proc(x)
    local t;
    t:=evala(Trace(x));
    if t=x then
      RETURN(x);
    else
      RETURN(t-x);
    fi;
  end;
  RETURN(map(f,L));
end:
#--- divcalc.mpl (end) ---#