Here is Magma code that gets you the answer in a few seconds. I made a special case for the bad primes, and did them by hand.
_<x> := PolynomialRing(Rationals());
f5 := 344 + 3106*x - 1795*x^2 - 780*x^3 - x^4 + x^5;
g24 := 14488688572801 - 2922378139308818*x^2 + 134981448876235615*x^4 -
1381768039105642956*x^6 + 4291028045077743465*x^8 -
2050038614413776542*x^10 + 287094814384960835*x^12 -
9040633522810414*x^14 + 63787035668165*x^16 - 158664037068*x^18 +
152929135*x^20 - 50726*x^22 + x^24;
K := NumberField(f5);
_,D := IsSquare(Integers()!Discriminant(f5));
prec := 30;
CHAR_TABLE := CharacterTable(GaloisGroup(g24));
chi := CHAR_TABLE[2];
BAD_FACTORS :=
[ <2,Polynomial([1,-1,1])>,
<3,Polynomial([1,-ComplexField(prec)!chi[9],1])>,
<5,Polynomial([1,0,1])>,
<7,Polynomial([1,0,1])>,
<71,Polynomial([1,0,1])>,
<137,Polynomial([1,1,1])>,
<163,Polynomial([1,-ComplexField(prec)!chi[5],1])>,
<1951,Polynomial([1])>,
<16061,Polynomial([1,-2,1])>,
<889289,Polynomial([1,-ComplexField(prec)!chi[8],1])> ];
BAD := [bf[1] : bf in BAD_FACTORS];
FACTORS := [bf[2] : bf in BAD_FACTORS];
function LOCAL(p,d : Precision:=prec)
if p in BAD then return FACTORS[Position(BAD,p)]; end if;
R := Roots(ChangeRing(f5,GF(p)));
if #R eq 1 then return Polynomial([1,0,1]); end if;
if #R eq 2 then
ord := Lcm([Degree(f[1]) : f in Factorization(Polynomial(GF(p),g24))]);
return Polynomial([1,ord eq 3 select 1 else -1,1]); end if;
if #R eq 5 then
ord := Lcm([Degree(f[1]) : f in Factorization(Polynomial(GF(p),g24))]);
return Polynomial([1,ord eq 1 select -2 else 2,1]); end if;
r := Roots(ChangeRing(f5,GF(p^5)));
x := r[1][1];
prod := GF(p)!&*[x^(p^i)-x^(p^j) : j in [(i+1)..4], i in [0..4]];
wh := prod eq GF(p)!D;
ord := Lcm([Degree(f[1]) : f in Factorization(Polynomial(GF(p),g24))]);
if ord eq 10 then class := wh select 8 else 9; // compatible with FACTORS
else class := wh select 6 else 5; end if;
return Polynomial([1,-ComplexField(prec)!chi[class],1]);
end function;
L := LSeries(1, [0,0], 1951^2, LOCAL : Precision:=prec);
// s->1-s, Gamma(s/2)^2
psi := DirichletGroup(1951, CyclotomicField(10)).1;
p1951 := Polynomial([1,-ComplexField(prec)!CyclotomicField(5).1]);
TP := TensorProduct(L, LSeries(psi : Precision:=prec), [<1951, 1, p1951>]);
CheckFunctionalEquation(TP);
Here is the special values:
ev := Evaluate(TP,0); // 2-1.453085056...
rel := PowerRelation(ev,4 : Al:="LLL");
NF := NumberField(rel);
Q5<zeta5> := CyclotomicField(5);
assert IsIsomorphic(NF,Q5);
Q5!NF.1;
So $L(\rho,0)=-4\zeta_5(1+\zeta_5)$ for Marty. I get $L(\rho_0,-1)=32(48723\sqrt{5} - 778741)$ as an algebraic. I get $L(\rho,-2)=8800\zeta_5^3 - 14444\zeta_5^2 + 35604\zeta_5 + 17412$ with more precision. I determined the TensorProduct factor at 1951 via trial and error, making the obvious guesses until one worked (the failure is at 100-110 digits). With this, I take it to 240 digits and I can even get $$L(\rho,-4)=-18475535360\zeta_5^3 - 11142861380\zeta_5^2 - 12091894020\zeta_5 - 7107607296$$ and $$L(\rho,-6)=25255057273186244\zeta_5^3 - 1015274469604000\zeta_5^2 - 15695788409197884\zeta_5 +
9459547822189412$$ The precision can go higher if you want more.
Finally, the Maass form:
function MaassEval(L,z)
x:=Real(z); y:=Imaginary(z);
printf "Using %o coefficients\n", Ceiling(11/y);
C := LGetCoefficients(L,Ceiling(11/y));
pi := Pi(RealField());
a := Sqrt(y)*&+[C[n]*KBessel(0,2*pi*n*y)*Sin(2*pi*n*x) : n in [1..#C]];
return a;
end function;
zz:=0.0001+0.0001*ComplexField().1;
MaassEval(TP,zz);
// Using 110000 coefficients
// -1.71477211817772949974178783985E-8 + 9.01673609747756708674470686948E-9*i
MaassEval(TP,zz/(1951*zz+1));
// Using 161297 coefficients
// -1.71477211817772949974179078240E-8 + 9.01673609747756708674496293450E-9*i
The basic and classical methods, apart from brute force, are
continued fraction expansions (regular, nearest integer, etc.) or, equivalently,
some form of reduction theory for indefinite binary quadratic forms;
computing many elements of small norm in a quadratic number field, which often
is a lot more effective; the technique used here is also used for factoring
integers.
For a detailed algorithmic description see Jacobson & Williams (Solving the Pell Equation)
or Buchmann & Vollmer (Binary Quadratic Forms).
In addition, you can compute a power of the fundamental unit from the class number formulas, which essentially consists in taking norms of suitable cyclotomic units.
Kronecker has shown how to solve the Pell equation using elliptic and modular functions,
and Girstmair (Kronecker's solution of the Pell equation on a computer {Kroneckers Lösung der Pellschen Gleichung auf dem Computer], Math. Semesterber. 53, 45-64 (2006)) has shown that it can be made to work in practice.
You can also imitate the theory of descent on elliptic curves; I have sketched connections with classical tricks in some preprints on higher descent on Pell conics.
Best Answer
Based on the comments, it looks like this is not a question specific to Pell's equation, and that you just want to evaluate a single binomial with big inputs as quickly as possible.
If you check the equation directly using fast multiplication algorithms (e.g., Schönhage-Strassen), you can expect the calculation to require about $(\log z )\cdot (\log \log z)$ operations, where $z = \max \{ x, y, D \}$.
If you want a way to quickly find a negative answer, you can check relative sizes and leading digits, then try reduction modulo small integers. If you start with a random negative answer, you can expect it to be eliminated after at most a few divisions (i.e., about $\log z$ operations).
To find a positive answer using modular arithmetic, you can use the Chinese remainder theorem. To prove that the identity holds, it suffices to check it modulo $n$, for $n$ ranging over a collection of positive integers whose least common denominator is larger than $x^2$ and $Dy^2$. It is common to check modulo a large collection of small primes, and this will require about $\log z$ primes and $(\log z)^2$ operations. Another natural choice with a binary computer is Fermat numbers, of the form $2^{2^n}+1$, since the division-with-remainder can be optimized - this ends up looking a lot like direct calculation.
In summary, the advantage of checking modulo small primes is that it lets you quickly eliminate negative answers, and the disadvantage is that (if I'm not mistaken) it is roughly quadratically slower than direct calculation when you have a positive answer. You can choose your method depending on exactly what sort of calculation you plan to do.