On Pell's equation A post was made (Reduction from factoring to solving Pell equation) seeking clarification to solving $$x^2-Dy^2=1$$ to factoring when $D>0$. An answer was posted stating that to factor $N$, it suffices to do trials of $D=N,2N,\dots$.
My questions:
1. How many such trials are needed?
2. What is complexity for each trial?
3. Is there an efficient method to solving equations of form $x^2+Dy^2=m$ when $D>0$?
 A: Partial answer to (2).
The smallest solution in some cases might be prohibitively large.
Check this paper pp 3,4.
It uses negative Pell, but its solution is smaller that the wanted.
The paper shows the smallest solution might be exponential in $D$.

Efficient solution to (3) would give probabilistic factoring algorithm.
Computing square roots modulo composites gives probabilistic factoring,
since for composite $m$ if a root exists, in general there are more than two
by the chinese remainder theorem.
If you solve $r_1^2\equiv r_2^2 \pmod{m}$ and $r_1 \not \equiv \pm r_2 \pmod{m}$,
you know proper divisor.
Let $n$ be integer you want to factor. Pick random integer $d$
and set $D= (-d^2) \mod m$. Set $m=d^2+D$. Observe that $n$ divides $m$.
There is at at least one solution $(d,1)$ and there might be more.
Solutions modulo $m$ might give divisor of $n$, but not necessarily.
Solution means $x^2 \equiv -D y^2 \pmod{m}$. If $y$ is not invertible modulo $n$ this might divisor.
If it is invertible, you have found square root of $-D$ and already know one random.
If this fail to give divisor, of $n$ chose another $d$ and repeat.
Here is pari/gp toy implementation:
? n=11*17;d=52;D=-d^2%n;m=d^2+D;K=bnfinit(x^2+D);nor=bnfisintnorm(K,m)
%18 = [x + 52, -2*x + 49, -2*x - 49, -x + 52]
? t=49/2%n;[gcd(d-t,n),gcd(d+t,n)]
%19 = [11, 17]
? n=11*17;d=53;D=-d^2%n;m=d^2+D;K=bnfinit(x^2+D);nor=bnfisintnorm(K,m)
%26 = [x + 53, -4*x - 8, 4*x - 8, -x + 53]
? t=8/4%n;[gcd(d-t,n),gcd(d+t,n)]
%27 = [17, 11]

To solve it, you must factor both $D,m$, since similar argument applies for
composite $D$.

If you need this in practice, pari/gp's bnfisintnorm() does it,
but has additional to factoring overhead.

I believe even if factoring is easy this would be hard for negative $D$,
since it is a Pell equation for $m=1$ and it is not easy for prime $D$.
For $D$ positive likely it will be efficient assuming factoring oracle.
