genpellsolvealt := proc(D::posint, N::integer) local P, Q, a, A, B, i, sols, f; if type(sqrt(D), integer) then error("D must be a nonsquare integer"); end if; if not(0 < N^2 and N^2 < D) then error("0 < N^2 < D not satisfied"); end if; P := 0; Q := 1; a := floor(sqrt(D)); A := 1, a; B := 0, 1; sols := {}; for i from 1 do f := sqrt(N/(A[2]^2-D*B[2]^2)); if type(f, integer) then sols := sols union {[f*A[2], f*B[2]]}; end if; P := a*Q - P; Q := (D - P^2)/Q; a := floor((P+sqrt(D))/Q); A := A[2], a*A[2]+A[1]; B := B[2], a*B[2]+B[1]; if Q = 1 and i mod 2 = 0 then break; end if; end do: return sols; end;