Root to digits of π

Note: I'm not using the leaderboard format, for reasons that will become apparent later.

Let's suppose we have a rational approximation $d = \frac a b$ to some real number that we want to be, modulo 1, close to the square root of a natural number - in this case, $.314159...$ or $\frac\pi{10}$. Then we want $N := (x + d)^2$ to be as close to an integer as possible, and in fact

$$ N = (x + d)^2 = (x + \frac a b)^2 = x^2 + 2 x \frac a b + {\frac a b}^2 = x^2 + \frac{2 x a b + a^2}{b^2} $$

Since we know $x^2$ is an integer and we want $N$ to approximate one, we can ignore that and deal with $N' := \frac{2 x a b + a^2}{b^2}$. For this to be an integer, we need

$$ b^2 \text{ divides } 2xab + a^2 \iff -a^2 = 2 x a b \text{ mod } b^2 \iff x = -a^2 \cdot \left(2 a b\right)^{-1} \text{ mod } b^2 $$

except this can't quite work: $2 a b$ isn't coprime to $b^2$, so $\left(2 a b\right)^{-1}$ won't exist modulo $b^2$ - and indeed, if this did work, we would be able to find integer square roots with arbitrary finite decimal expansions, which can't happen - $(x + .5)^2$ is never an integer, for instance. What we can do, however, is cheat slightly and add a small perturbation - take $\left(2 a b - 1\right)^{-1}$, for instance (which could also share a factor with $b^2$ and not be divisible, but sufficiently nice choices of $a$ and $b$ will prevent this)

Note 2: I don't actually know why this works - upon further thought, due to how modular arithmetic works, a small change in $2 a b$ could result in a large change in the resulting value of $x$? But for some reason it in fact works very well, as will become apparent later.

This gives us a value for $x$, and from that we can compute $N = x^2 + \left\lfloor\frac{2 x a b + a^2}{b^2}\right\rfloor$, where the floor division is necessary to fix our approximations.

For an example of this, let's approximate $\pi$ as $\frac{22}7$, i.e. $d = \frac{\frac{22}7}{10} = \frac {11}{ 35}$, i.e. $a = 11, b=35$. Then we get $ x = 11^2 \cdot (2 \cdot 11 \cdot 35 - 1)^{-1} \text{ mod } 35^2 $, which comes out to $1034$, after which $N$ comes out as $1069806$ and $\sqrt N$ as $1034.314265$ - pretty good!

Note 3: My original calculation used $-a^2$, but for some reason $a^2$ works better - not entirely sure why.

Next, for a larger example, let's try plugging in the first 200 decimal places of $\pi$, and after a bit my code spits out

937523293002658776348704769842539305648144222902391041804198307023876611679007339563904576704517047026597642317245013563261350766724787513251817694301543229194444616131823525148446741635864093056342088369126394578418147892002004186052415814273659670054499001041125811242784584390954000298648838953612272499808590813639231701499889493344119553570944163927982067912689004203386142800764496807779136844467819936800975933272481669347865507385952118302992371399189910547547628198410940202325281436034655796414880027985007133712574299608836521224536654484339688199950743647254454272876064422303005862463315882536889249083787990355940641001191524146159193403485502061273254506726884093838731770801055842611988398619589116686583248837626948667051332081102169659142504795515415306915108033323241905620107674

whose square root's decimal expansion looks like

[integer part].314159265358979323846264338327950288419716939937510582097494459230781640628620899862803482534211706798214808651328230664709384460955058223172535940812848111745028410270193852110555964462294895493038189999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999848068430337174906334649751655113636032420403882399880832538599596915793497396586679196815633320326311365700306099279385136990931183996110519484150914211076274575736266077990160878090363392867830509955

which is 199 digits of $\pi$. Interestingly, we actually get $400$ digits of precision out of this - the $200$ digits following the digits of $\pi$ are all $9$s, which sorta makes sense considering that the term we introduce an error into, $2 a b$, has 400 digits of precision. In any case, this is accurate to 199 digits (the final digit is rounded down, so isn't technically correct), and also makes the online verifier hang.

So, let's take this to its logical conclusion: I downloaded 4 million digits of pi and plugged (most of) them into my script. My code had to be optimised a fair bit - in particular, modulo operations are absurdly slow so I had to use string slices (later made redundant), and computing the modular inverse of something that large is effectively impossible, so I instead simply divide $(b^2)^2$, roughly, by the value that needs to be inverted, which works almost as well. I later ended up switching to gmpy2, which was a massive speedup and allowed me to actually verify things.

The result is a 12566368-digit number that has the first 3141592 digits after its square root's decimal point equal to the corresponding digits of pi.

---

As for how far this can go, there isn't really any practical limit. With gmpy2, my laptop can process and verify the 3141592-digit case in about 3 seconds, and a 10-million-digit case in about 13. Algorithmically, I can produce a 4N-digit number that gives N decimal places of accuracy, and it's absolutely possible to improve on that constant factor (although a simple counting argument shows that you can't do any better asymptotically). It's difficult to get an idea of the algorithm's speed since most of it comes down to the underlying bigint representations, but it does seem relatively close to linear in the number of digits. Just for fun, I let my computer loose on a 100-million-digit testcase and it took 3 minutes, so with enough memory, scores of a billion and beyond definitely seem feasible.