cf(n) = {
	my(f=s=sqrtint(n), k=1, P=[0], Q=[1], R=[f]);
	while(1,
		P=concat(P, f*Q[k]-P[k]);
		Q=concat(Q, (n-P[k+1]^2)/Q[k]);
		k++;
		for(i=1, k-1, 
			if(P[i] == P[k]&&Q[i] == Q[k],
				return(R[1..k-i]);
			);
		);
		R=concat(R, f=(P[k]+s)\Q[k]);
	);
}


funit(n) = {
	my(v=cf(n)); my(q=v[#v]);
	forstep(k=#v-1, 1, -1, q=1/q+v[k]);
	v=[numerator(q), denominator(q)];
	if(v[1]^2-n*v[2]^2<0, v[2]*=2*v[1]; v[1]=2*v[1]^2+1);
	v;
}

nbs(m) = {
	my(b, c=0, r, t, u, w, x);
	if(issquarefree(m),
		fordiv(m, d,
			if(d == m,
				w=funit(d); c+=issquare(w[1])+issquare(w[1]^2+d*w[2]^2),
				if(d > 1,
					b=m/d; w=funit(d); t=w[1]; u=w[2]; r=[t%b, u%b];
					until(r == [t%b, u%b],
						c+=issquare(t/b);
						x=u*w[1]+t*w[2]; t=t*w[1]+d*u*w[2]; u=x;
					);
				);
			);
		);
	);
	c;
}

/* for A335785 */
lista1(nn) = {
	my(nb=0);
	for(m=2, nn,
		if(nbs(m), nb++; print1(m, ", "));
	);
	print; print("nb=", nb);
}


/* for A335715 */
lista2(nn) = {
	my(nb=0);
	forcomposite(m=2, nn,
		if(nbs(m) > 1, nb++; print1(m, ", "));
	);
	print; print("nb=", nb);
}

/* for A336145 */
lista3(nn) = {
	my(nb=0, rev);
	for(m=2, nn,
		if(nbs(m), if(rev == m-1, nb++; print1(m-1, ", ")); rev=m);
	);
	print; print("nb=", nb);
}