4
$\begingroup$

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?

$\endgroup$
  • $\begingroup$ Have you read the Wikipedia page on square root algorothms? $\endgroup$ – Bill Dubuque Mar 26 '14 at 14:19
  • $\begingroup$ Ideas: 1) uses an arbitrary precision package like mpmath 2) Assume $a \ge 0, b > 0$, generate a sequence of fractions $p/q$ approximating $\sqrt{2}$ using continued fraction When $q$ is big enough (about the size of square root of numerator of $b$), then $\lfloor a + b \sqrt{2}\rfloor = \lfloor a + b \frac{p}{q}\rfloor$, you can use the criterion that successive terms in continued fraction approximation give same floor as a terminating condition. $\endgroup$ – achille hui Mar 26 '14 at 14:39
  • $\begingroup$ Bill, I had not. I see that one might potentially approach the problem by applying one of the square root algorithms on that page (maybe the "Babylonian method") to calculate a series of approximate square roots of $2b^2$, which when it has converged on the square root to within an integer gives a candidate for $\lfloor b\sqrt{2}\rfloor$ which should not be far off the correct answer. $\endgroup$ – Hammerite Mar 26 '14 at 14:41
4
$\begingroup$

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.

$\endgroup$
  • $\begingroup$ This is some great advice, thank you. Having the method related to the Pell numbers is interesting. I will leave off accepting for a little longer but will probably accept this answer. $\endgroup$ – Hammerite Mar 26 '14 at 16:24
  • $\begingroup$ On further thought, there is no gain in finding the smallest $e,f$ for a given $b$. The above only increases $e,f$ by a factor $3+\sqrt 2$ each step. You could take three steps at a time by $e=7e+10f,f=5e+7f$ and get a larger factor. You could also store a number of $e,f$ pairs and not compute them each time. This will speed you up for larger $b$ $\endgroup$ – Ross Millikan Mar 26 '14 at 19:05
  • $\begingroup$ Note the the update equations I give for $e,f$ are intended to be applied simultaneously-you use the old value on both right sides. Also this only works for positive $b$. For negative $b$, it is probably best to invert it, find the floor, negate it, and subtract $1$. You could sort out the changes, but then would have two copies of the code. $\endgroup$ – Ross Millikan Mar 28 '14 at 15:13
  • $\begingroup$ Yes, I understood. Are you assuming that $b$ is an integer in your first step (comparison with $g$)? I suppose (writing $b = c/d$) it is possible to find $\lfloor c\sqrt{2}\rfloor$ and then divide by $d$ and correct for any error introduced. $\endgroup$ – Hammerite Mar 28 '14 at 16:29
  • $\begingroup$ You are correct, this is written for integer $b$. I missed that it can be rational. Yes, you can write $\lfloor \frac {c\sqrt 2}d \rfloor=\left\lfloor \frac {\lfloor c\sqrt 2\rfloor }d \right\rfloor$ and use the above. $\endgroup$ – Ross Millikan Mar 28 '14 at 16:36

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.