(* We compute in Mathematica, given a j-invariant in Q(t) and a given choice of nine nonzero field elements p1,...,p9 
satisfying the conditions of Section 3, the coefficients cijk which give the honeycomb embedding of the given elliptic curve. 
We demonstrate this code on the j-invariant 1/t^24 + t and with the choice of nine points given in Section 4. *)

(* Accuracy *)

m=25;

(* j invariant*)

jinput=Series[1/t^24+t,{t,0,m}];

(* choice of lattice length of the three bounded non-hexagon segments *)

b=1;


(* Valuation function *)

v[ser_]:=Exponent[ser,t,Min];

(* Order of 1/j*)

order=-1*v[jinput];

Clear[q,a];

(* Compute the series s3, s5, a4, a6, Delta, and j as functions of q.  
See Silverman Advanced Topics in the Arithmetic of Elliptic Curves, p.425.  *)

s3=Sum[n^3*q^n/(1-q^n),{n,1,m}];

s5=Sum[n^5*q^n/(1-q^n),{n,1,m}];

a4=-5*s3;

a6=-(5*s3+7*s5)/12;

Delta=-a6+a4^2+72 a4*a6-64 a4^3-432 a6^2;

jformula=Normal[Series[jseries=((1-48*a4)^3)/Delta,{q,0,m}]];

(* Write q as a power series in t of valuation -val(j)*)

q=Sum[a[i] t^i,{i,order,order+m}]+a[order+m+1] O[t]^(order+m+1);

(* Solve for q in terms of the input j-invariant;output stored in qout*)

qout=q/.Solve[LogicalExpand[jformula==jinput]][[1]];

(* Let mathematica pick a 6th root of q*)

r=qout^(1/6);

(* Pick a number u satisfying val(1-u)=b*)

u=1-t^b;

(* Pick the following 9 points in K^*.  As in Section 3, they produce a symmetric honeycomb embedding. 
Many other choices of p_i's are possible. *)

p0=1/(r*u);p1=1;p2=r*u;p3=r/u;

p4=r^2;p5=r^-3*u;p6=r^3/u;

p7=r^-2;p8=u/r;


(* Define the fundamental theta function *)

th0[x_]:=Expand[Product[1-qout^i*x,{i,1,m}]*Product[1-qout^i/x,{i,0,m}]];

(* Define the theta function with respect to the scalar a in K^ *)

th[a_,x_]:=th0[x/a];

(* Using theta functions,we construct meromorphic, q-periodic functions with divisors 
p3 + p4 + p5 - p0 - p1 - p2 and p6 + p7 + p8 - p0 - p1 - p2 *)

f[x_]:=th[p3,x]*th[p4,x]*th[p5,x]/(th[p0,x]*th[p1,x]*th[p2,x]);
g[x_]:=th[p6,x]*th[p7,x]*th[p8,x]/(th[p0,x]*th[p1,x]*th[p2,x]);


(* To compute an implicit equation, pick 9 points on the curve,say of the form (f(t^i),g(t^i)) and solve for the coefficients. 
Alternatively, use the explicit formulae given in Section 3. *)

tab=Table[{f[t^s]^3,f[t^s]^2*g[t^s],f[t^s]*g[t^s]^2,g[t^s]^3,f[t^s]^2,f[t^s]*g[t^s],g[t^s]^2,f[t^s],g[t^s]},{s,1,9}];

(* Assuming c003 = 1, invert the matrix above to solve for the rest of the coefficients *)

{c300,c210,c120,c030,c201,c111,c021,c102,c012,c003} =Append[Inverse[tab].{-1,-1,-1,-1,-1,-1,-1,-1,-1},1];


(* Our output is the ten coefficients of the honeycomb embedding associated to p1,...,p9. *)

{c300,c210,c120,c030,c201,c111,c021,c102,c012,c003}