Prime counting function
One way to find primes is to find a polynomial time algorithm to compute [math]\displaystyle{ \pi(x) }[/math], the number of primes less than x, to reasonable accuracy (e.g. if we can do better than [math]\displaystyle{ 10^{k/2} }[/math] accuracy for k-digit x, we can break the square root barrier.) We don't necessarily have to do this for all x; just having a targeted x for which we can show [math]\displaystyle{ \pi(x+y)-\pi(x) \gt 0 }[/math] for some [math]\displaystyle{ x \sim 10^k }[/math] and [math]\displaystyle{ y=o(10^{k/2}) }[/math] would suffice.
Now, perhaps instead of trying to prove that intervals like [x, x + (\log x)^A] contain primes unconditionally, as Emmanuel suggested, we should first try to be much less ambitious and aim to show that *some* interval [y, y + \sqrt{x}] with y \in [x,2x], contains a prime number that we can discover computationally.
How? Well, let’s start by assuming that we can computationally locate all the O( \sqrt{x} \log x) zeros in the critical strip up to height \sqrt{x}. Then what we can do is some kind of “binary search” to locate an interval [y, y + \sqrt{x}] containing loads of primes: say at the ith step in the iteration we have that some interval [u,v] has loads of primes. Then, using the explicit formula for \psi(z) = \sum_{n \leq z} \Lambda_0(n) = z - \sum_{\rho : \zeta(\rho)=0} z^\rho/\rho (actually, the usual quantitative form of the identity), we can decide which of the two intervals [u, (u+v)/2] or [(u+v)/2, v] contains loads of primes (maybe both do — if so, pick either one for the next step in the iteration). We keep iterating like this until we have an interval of width around x^{1/2 + \epsilon} or so that we know contains primes.
Ok, but how do you locate the zeros of the zeta function up to height \sqrt{x}? Well, I’m not sure, but maybe one can try something like the following: if we can understand the value of \zeta(s) for enough well-spaced, but close together, points up the half-line, then by taking local interpolations with these points, we can locate the zeros to good precision. And then to evaluate \zeta(s) at these well-spaced points, maybe we can use a Dirichlet polynomial approximations, and then perhaps apply some variant of Fast Fourier Transforms (if this is even possible with Dirichlet polynomials, which are not really polynomials) to evaluate them at lots of values s = 1/2 + \delta quickly — perhaps FFTs can speed things up enough so that the whole process doesn’t take more than, say, 10^{k/2} polylog(k) bit operations. Keep in mind also that our Dirichlet polynomial approximation only needs to hold “on average” once we are sufficiently high up the half-line, so it seems quite plausible that this could work (and for s near to 1/2 we would need to be more careful, and get the sharpest approximation we can, because those terms contribute more in the explicit formula).
I realize that there are far better methods than what I describe for locating zeros of the zeta function, and in the end it might be better to use one of these other methods, assuming what I write above can work.
