pellsolve := proc(D::posint)
  local P, Q, a, A, B, i;
  if type(sqrt(D), integer) then
    error("D must be a nonsquare integer");
  end if;
  P := 0;
  Q := 1;
  a := floor(sqrt(D));
  A := 1, a;
  B := 0, 1;
  for i from 1 do
    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 A[1], B[1];
end;