A general method to efficiently calculate the floor of an element of $\mathbb{Q}[\sqrt{2}]$

Question

I had to decide whether to post this question here on the Mathematics Stack Exchange or on Stack Overflow, but I decided that the question was essentially a mathematical one despite being inspired by a programming problem, so I'm putting it here.

I want to write a Python class that represents numbers in $\mathbb{Q}[\sqrt{2}]$, that is to say, members of this class are numbers of the form $a + b\sqrt{2}$ where $a$ and $b$ are rational numbers. I want to fit this class into Python's numeric hierarchy, so it should fit into the tower of Abstract Base Classes in between Real and Rational. That imposes the requirement to implement certain class methods; for example, so that you can add two numbers together, and so that a member of the class can be interrogated to find out whether it is bigger or smaller than some other real number.

One of the things that you need to do to make a class a subclass of Real is to implement the __floor__() method. This is the method that tells Python how to calculate the floor function applied to your custom object. Well, it seems there are a number of related magic methods (__ceil__() and __trunc__() and __round__() and maybe others?), but the essential problem is the same.

Now the immediate thing that comes to mind is to calculate the floor by converting the object to a floating-point number. But there are two problems with this. Whereas Python's integers are arbitrary-precision, and its Fractions are essentially an integer numerator paired with an integer denominator (and therefore also arbitrary-precision), floating-point numbers are limited in precision. So the problems are:

1. The wrong answer might be obtained for some numbers (it might be off by one). The floating-point representation of $\sqrt{2}$ used is by necessity an approximation. If you try to find the floor of a number that is very very close to an integer then you might end up with the floating point representation just on the wrong side of that integer, and then you get the wrong result.
2. Some numbers will be too big and the "convert to a float" step will not work.

>>> x = 10 ** 400
>>> x
10000000000000000000000000000000000000000000000000000000000000000000000000000000
00000000000000000000000000000000000000000000000000000000000000000000000000000000
00000000000000000000000000000000000000000000000000000000000000000000000000000000
00000000000000000000000000000000000000000000000000000000000000000000000000000000
00000000000000000000000000000000000000000000000000000000000000000000000000000000
0
>>> float(x)
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
OverflowError: long int too large to convert to float

The question then is how to efficiently find the floor of one of these numbers and definitely get it right.

Let $x = a + b\sqrt{2}$, where $a, b \in \mathbb{Q}$. We first consider that we can easily find out what the sign of $x$ is. The sign of $x$ is $0$ if and only if $a = b = 0$, because $\sqrt{2}$ is irrational. Suppose now that the sign of $x$ is not zero. If $a$ and $b$ have the same sign, then that is the sign of $x$; if they have different signs then $y = a - b\sqrt{2}$ has the same sign as $a$, and the sign of $x$ is the sign of $xy$ divided by the sign of $y$. ($xy$ is $a^2 - 2b^2$, and we know how to find the sign of this because it is a Fraction). Since we know how to find out what the sign of $a + b\sqrt{2}$ is, it is also easy to compare $a + b\sqrt{2}$ with $c + d\sqrt{2}$ and find out which is bigger, because that is just finding the sign of $(a - c) + (b - d)\sqrt{2}$.

Since we can compare numbers in this way, it is clear that we can at least identify for a given integer $n$ whether $n < x$ or not. In principle this means that we can find the floor of $x$. If $x$ is positive then we try $n = 1$, $n = 10$, $n = 100$, $\ldots$ until we find an $n = 10^k$ that is bigger than $x$ and then once we find it we can do a binary search of $\{0, 1, 2, \ldots, n\}$ to identify the answer. If $x$ is negative then we do the same thing but with $n = -1$, $n = -10$, $n = -100$, $\ldots$.

We also observe that we can essentially reduce the problem to finding the floor of $b\sqrt{2}$. It is straightforward to reduce the problem to that of finding the floor of $a' + b\sqrt{2}$ where $a' = a - \lfloor a\rfloor$ is the fractional part of $a$. Then if we just ignore $a'$ our answer (after doing "the hard part") will be off by at most one, which we can easily check for and correct.

An improvement over the first method above might be to estimate $\lfloor b\sqrt{2}\rfloor$ by using an approximation of $\sqrt{2}$ as a Fraction. It should then be possible to determine the maximum possible error by checking how many digits are in the result. (The logarithm function in Python is for floating-point numbers, but in the case of big numbers we could cheat by asking Python to represent the number as a string and then asking it how many characters are in the string.) This gives a range of values which can then be checked for the true answer, like the previous method but without the initial geometric search for an upper bound to the range.

Can you suggest a better way than this?

Answer 1

Your approach using approximating fractions seems like a good one. You can use the solutions to the Pell equation to your advantage. If $\frac ef$ is one approximation to $\sqrt 2$, the next one is $\frac {e+2f}{e+f}$ and the error is of the opposite sign. You will probably be asked most often for floor $(b \sqrt 2)$ for small $b$, so I would emphasize speed in that area, while being able to handle arbitrarily large $b$. I would pick some fractional approximation and store it, also storing the largest $b$ that this approximation is guaranteed to work for. For example, we might use $\frac {1393}{985}$, which is a little low. It first fails at $b=3363$, where $b\sqrt 2 \approx 4756.00021,$ while $\frac {1393}{985}\cdot 3363\approx 4755.998985$ My approach would then be:

Store $e,f,g$, where $\frac ef$ is an approximation for $\sqrt 2$ with the error in a known direction. I will assume $\frac ef \lt \sqrt 2$ and $g$ is the first point of failure.

If $b \lt g,$ return $\lfloor \frac {eb}f \rfloor$
else, compute $\lfloor \frac {eb}f \rfloor, \lfloor \frac {(e+1)b}f \rfloor$ If they agree return it
else, $e:=e+2f,f:=e+f$ Now $\frac ef$ is too large
Compute $\lfloor \frac {eb}f \rfloor, \lfloor \frac {(e-1)b}f \rfloor$ If they agree, return it.
else, keep increasing $e,f$ until the two floors agree. You have to change the sign on the $1$ each time.

The nice thing is that if you barely go past an integer in one direction, you have almost a full unit of space the other direction.

I think I would choose $\frac ef$ to have the largest $e$ so that $eg$ fits in $64$ bits. That should let you succeed in the first test most of the time, yet have that computation fast.