genpellsolve := proc(D::posint, N::integer)
  local t, u, L1, L2, sols, x, y;
  if type(sqrt(D), integer) then
    error("D must be a nonsquare integer");
  end if;
  t, u := pellsolve(D);
  if N > 0 then
    L1 := 0;
    L2 := floor(sqrt(N*(t-1)/(2*D)));
  elif N < 0 then
    L1 := ceil(sqrt(-N/D));
    L2 := floor(sqrt((-N)*(t+1)/(2*D)));
  else
    return {[0, 0]};
  end if;
  sols := {};
  for y from L1 to L2 do
    x := sqrt(N+D*y^2);
    if type(x, integer) then
      sols := sols union {[x, y]};
      if x^2+y^2*D mod N <> 0 or 2*x*y mod N <> 0 then
        sols := sols union {[-x, y]};
      end if;
    end if;
  end do;
  return sols;
end;