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); }