Block[{s,u,q,p,a,n,k, d, x},
s=18;
u[0]=k;
u[1]=10; 
u[2]=2;
u[3]=12;
u[4]=6;
u[5]=1;
u[6]=3;
u[7]=4;
u[8]=3;
u[9]=12;
(*remaining coefficients are symmetric*)
Do[u[s-i]=u[i],{i,1,s/2-1}];
u[s]=2k;
u[s+1]=u[1];
(*define sequences p[i] and q[i] by recurrences*)
q[-1]=0;
q[0]=1;
Do[q[i]=u[i]q[i-1]+q[i-2],{i,1,s}];
p[-1]=1;
p[0]=u[0];
Do[p[i]=u[i]p[i-1]+p[i-2],{i,1,s}];
Print["The equation for D is ",Simplify[{Sqrt[d]*((k+Sqrt[d])*q[s-1]+q[s-2])==((k+Sqrt[d])p[s-1]+p[s-2])}]];
Print["The solutions of the equation are ",Reduce[{Sqrt[d]*((k+Sqrt[d])*q[s-1]+q[s-2])==((k+Sqrt[d])p[s-1]+p[s-2])},{k,d},Integers]/.C[1]->(-x),"where x is non-negative."];
k=32229463807+61360840419 x;d=1038738337292873457072+3955253970914236253087 x+3765152736925984095561 (x^2);
n[i_]:=Simplify[p[i]^2-d q[i]^2];
Print["Negative norms of convergents for i = 0, 2, 4, ..., 16: ",
Table[n[i],{i,0,s-2,2}]];
Print["Hence alpha_2 has the largest negative norm (for every x)."];
(*define norm of alpha_{i,r}*)
n[i_,r_]:=Simplify[(p[i]+r p[i+1])^2-d(q[i]+r q[i+1])^2];
Print["The list of norms N(alpha_{i,r}) of indecomposables with r=u[i+2]/2 +- 1",TableForm[Table[n[i,r],{i,-1,s-1,2},{r, Ceiling[u[i+2]/2-0.9],Floor[u[i+2]/2+0.9]}],TableHeadings->{Table[i,{i,-1,s-1,2}], None}]];
Print["Hence alpha_{7,6} is the indecomposable element with largest norm (for every x)."];
Print["Two solutions of -N[2]|D-a^2 are given by a^2=d+n[2]n[1,6]=",Simplify[d+n[2]n[1,6]]];
Print["a_0:=838263647+1595948421 x is a positive solution that is less than n[2]/2 (for every x)."];
Print["The polynomial for D is reducible and factors as ",Factor[d], "=n[8]n[7,6]"];
Print["The polynomial for D mod 4 is ",PolynomialMod[d,4],", and for x equiv 2 mod 4 we get D equiv 2 mod 4." ];
x=4y-2;
Print["After substituting x=4y-2, the polynomial for D factors as ",Factor[d]," and N=-n[2] has no fixed prime divisor, as -n[2]=",Factor[-n[2]]];
Print["The polynomials n[2], n[8], n[7,6] are pairwise coprime, as their pairwise GCDs are ",PolynomialGCD[n[2],n[8]],", ",PolynomialGCD[n[7,6],n[8]],", ",PolynomialGCD[n[2],n[7,6]]];
]