/* 
   This program (http://sergey.ipmu.jp/degmir.gp)
   realizes (up to details) the algorithm from
   http://sergey.ipmu.jp/degmir.pdf
*/

\\ Symplectic pairing between two vectors on Z^2
\\ normalized in a way that pair([1,0],[0,1]) = 1
pair(u,v) = v[2]*u[1] - u[2]*v[1];

\\ Let u and v be two vectors in Z^2.
\\ refl(u,v) is symplectic reflection of v in u. irefl is inverse.
\\ mupl(u,v) is a Piecewise-Linear mutation of v in u. imupl is inverse.
refl(u,v) = v + pair(u,v) * u;
irefl(u,v) = v - pair(u,v) * u;
mupl(u,v) = v + max(0,pair(u,v)) * u;
imupl(u,v) = irefl(u,mupl(-u,v));
\\ Should use honest mutation, but too lazy to apply sl2 to laurent polynomial.
tr(u,v) = imupl(-u,v);

\\ Mutate a collection of vectors V in its n'th vector.
muplv(V,n) = vector(matsize(V)[2],ii,if(ii==n,-V[ii], tr(V[n],V[ii])));

\\ Let u be a vector in Z^2 and f a rational function of the variables (x,y)
\\ Next function mutates f in u = (m,n)
\\ x^a y ^b -> x^a y^b (1 + x^n/y^m)^(an+bm)
\\ In particular, transformation (x,y) -> (x,y/(1+x)) is  f -> mu([0,1],f)
\\ I was too lazy to find the inverse, so made ad-hoc inverse in tr.
mu(u,f) = subst(subst(subst(subst(f,x, X * (1+X^u[2]/Y^u[1])^(-u[1])),y, Y * (1+X^u[2]/Y^u[1])^(-u[2])),X,x),Y,y);
\\ mu([0,1],) explicitly
mu01(f) = subst(f,y,y/(1+x));

\\ Apply SL(2,Z) transformation m to the seed.
sl2(seed,m) = [vector(matsize(seed[1])[2],ii,1/m*seed[1][ii]), \
subst(subst(subst(subst(seed[2],x,X^m[1,1]*Y^m[1,2]),y,X^m[2,1]*Y^m[2,2]),X,x),Y,y) ];

\\ Mutate the seed in its n'th vector.
museed(seed,n) = 
{
nseed=[muplv(seed[1],n), mu(seed[1][n],seed[2])];
printseed(nseed);
return(nseed);
}

/* Now the functions to compute some invariants of the seed. */

\\Given Laurent polynomial f, guess its cluster collection V (by exhaustion).
findcc(f) = {vv=[];forvec(v=[[-4,4],[-4,4]],if(v!=0,g=f;while(islaurent(g=mu(v,g))==1,vv=matconcat(vv,[v]))));return(vv)}

\\ Vector of three Markov's numbers (up to sign)
markovnumbers(V) = -1/3 * [pair(V[1],V[2]),pair(V[2],V[3]),pair(V[3],V[1])];
\\ Skew-symetric matrix of generalized Euler-Markov's numbers
markovmatrix(v) = matrix(matsize(v)[2],matsize(v)[2],ii,jj,pair(v[ii],v[jj]));

\\ Check that u is monomial.
ismonomial(u) = (u==x^poldegree(u,x)*y^poldegree(u,y));
\\ Check that w is Laurent polynomial.
islaurent(w) = ismonomial(denominator(w));

\\ The value of Laurent polynomial f at the point (x,y)=(1,1) i.e. the sum of all coefficients of the Laurent polynomial.
bigvalue(f) = subst(subst(f,x,1),y,1);

\\ Some invariants of non-zero coefficients of Laurent polynomial w: [number, max, sum, sum of inverse]
inv(w) = 
{
u = w*denominator(w);
a=[0,0,0,0]; 
for(ii=0,poldegree(u,x),
for(jj=0,poldegree(u,y),
c=polcoeff(polcoeff(u,ii,x),jj,y);
if(c!=0,a[1]++);
a[2] = max(a[2],c);
a[3] += c;
if(c>0,a[4] += 1./c);
)); 
return(a) 
}


\\ Print some info about the seed - if f is Laurent polynomial and some invariants.
showmarkov = 1;
printseed(seed) = print(islaurent(seed[2]), /* " (denominator = ",denominator(seed[2]), ");", */ \
," ", if(showmarkov,markovnumbers(seed[1]),markovmatrix(seed[1]))," ", inv(seed[2]));

\\ We define the initial seeds for all 10 del Pezzo surfaces.
\\ The initial seed for P^2 - 3 vectors are defined by the Laurent polynomial x+y+1/x/y.
f0 = x + y + 1/x/y;
V0 = [[1,1],[-2,1],[1,-2]];
seedp2 = [V0,f0];
seed9 = seedp2;

\\ The initial seed for P^1 x P^1:
fq = x + y + 1/x + 1/y;
Vq = [[1,1],[1,-1],[-1,-1],[-1,1]];
seedp1p1 = [Vq,fq];

\\ Let Sd be del Pezzo surface of degree d, i.e. blowup of P^2 in (9-d) points, d=1...8

\\ The initial seed for S8.
f8 = x + y + 1/x/y + x*y;
V8 = [[1,0],[0,1],[-2,1],[1,-2]];
seed8 = [V8,f8];

\\ The initial seed for S7
f7 = x + y + 1/x + 1/y + x*y;
V7 = [[-1, -1], [-1, 1], [0, 1], [1, -1], [1, 0]];
seed7 = [V7,f7];

\\ The initial seed for S6
f6 = x+y+1/x+1/y+x*y+1/x/y;
V6 = [[-1, 0], [-1, 1], [0, -1], [0, 1], [1, -1], [1, 0]];
seed6 = [V6,f6];

\\ The initial seed for S5
f5 = x*y + 2*x + 2*y + x/y + y/x + 1/x + 1/y;
V5 = [[-1, -1], [-1, 0], [0, -1], [0, 1], [0, 1], [1, 0], [1, 0]]; 
seed5 = [V5,f5];

\\ The initial seed for S4
f4 = (1+x)^2*(1+y)^2/x/y-4;
V4 = [[-1, 0], [-1, 0], [0, -1], [0, -1], [0, 1], [0, 1], [1, 0], [1, 0]];
seed4 = [V4,f4];

\\ The initial seed for S3
f3 = (1+x+y)^3/x/y-6;
V3 = [[-1, 0], [-1, 0], [-1, 0], [0, -1], [0, -1], [0, -1], [1, 1], [1, 1], [1, 1]];
seed3 = [V3,f3];

\\ The initial seed for S2
f2 = (1+x+y)^4/x/y-12;
V2 = [[-1, 0], [-1, 0], [-1, 0], [-1, 0], [0, -1], [0, -1], [0, -1], [0, -1], [1, 1], [1, 1]];
seed2 = [V2,f2];

\\The initial seed for S1
f1 = (1+x+y)^6/x/y^2-60;
V1 = [[-1, 0], [-1, 0], [-1, 0], [-1, 0], [-1, 0], [-1, 0], [0, -1], [0, -1], [0, -1], [1, 1], [1, 1]];
seed1 = [V1,f1];


/*******************************************/

\\ Now just run some computations for P^2 and Markov's triplets.

seed00 = seedp2;
printseed(seed00);		\\ 1 1 1	-> 3
seed01 = museed(seed00,1);	\\ 1 1 2	-> 5
seed02 = museed(seed01,2);	\\ 1 2 5	-> 41
seed03 = museed(seed02,1);	\\ 1 5 13	-> 14803
seed04 = museed(seed02,3);	\\ 2 5 29	-> 580179113
seed05 = museed(seed03,2);	\\ 1 13 34	-> 79979617217
\\seed06 = museed(seed04,1);	\\ 5 29 433 	-> "131 digits"  [14 seconds to compute]