0

I am writing a program that handles numbers as large as 10 ** 100, everything looks good when dealing with smaller numbers but when values get big I get these kind of problems:

>>> N = 615839386751705599129552248800733595745824820450766179696019084949158872160074326065170966642688
>>> ((N + 63453534345) / sqrt(2)) == (N / sqrt(2))
>>> True

Clearly the above comparision is false, why is this happening?

Program code:

from math import *

def rec (n):
    r = sqrt (2)
    s = r + 2
    m = int (floor (n * r))
    j = int (floor (m / s))
    if j <= 1:
        return sum ([floor (r * i) for i in range (1, n + 1)])
    assert m >= s * j and j > 1, "Error: something went wrong"
    return m * (m + 1) / 2 - j * (j + 1) - rec (j)

print rec (1e100)

Edit:

I don't think my question is a duplicate of the linked question above because the decimal points in n, m and j are not important to me and I am looking for a solution to avoid this precision issue.

Tom Zych
  • 13,329
  • 9
  • 36
  • 53
razz
  • 9,770
  • 7
  • 50
  • 68
  • 1
    Floating point numbers have limited precision. If the standard range is insufficient you must use write your own code, perhaps using modules like decimal to help you. – President James K. Polk Jan 15 '17 at 19:52
  • 4
    Python has arbitrary-precision _integers_ but not arbitrary-precision _real numbers_. When you divide `N` and `N + 63453534345` by `sqrt(2)`, both are converted to IEEE double, which doesn't have enough precision to represent the difference. – zwol Jan 15 '17 at 19:52
  • Is there a way around this, I don't care about the decimals as I need these numbers floored eventually? – razz Jan 15 '17 at 19:54
  • You need to compute square roots of 100 digits numbers? What other operations do you need to compute? – President James K. Polk Jan 15 '17 at 19:59
  • Possible duplicate of [Is floating point math broken?](http://stackoverflow.com/questions/588004/is-floating-point-math-broken) – John Coleman Jan 15 '17 at 20:39
  • 1
    The module `decimal` is the standard module for higher (though obviously still finite) precision. – John Coleman Jan 15 '17 at 20:40
  • `sqrt(2)` is your problem. You can recast everything you do in terms of exact conditions involving integers. For example `a = floor(n*sqrt(2))` if and only if `a` is the largest integer satisfying the inequality `a**2 <= 2*n**2`, and such an `a` can be found with binary search (among other ways). The module `decimal` is probably less work than going this route. It can handle hundreds of decimal places with no problem. – John Coleman Jan 15 '17 at 21:15

1 Answers1

5

You can’t retain the precision you want while dividing by standard floating point numbers, so you should instead divide by a Fraction. The Fraction class in the fractions module lets you do exact rational arithmetic.

Of course, the square root of 2 is not rational. But if the error is less than one part in 10**100, you’ll get the right result.

So, how to compute an approximation to sqrt(2) as a Fraction? There are several ways to do it, but one simple way is to compute the integer square root of 2 * 10**200, which will be close to sqrt(2) * 10**100, then just make that the numerator and make 10**100 the denominator.

Here’s a little routine in Python 3 for integer square root.

def isqrt(n):
    lg = -1
    g = (1 >> n.bit_length() // 2) + 1
    while abs(lg - g) > 1:
        lg = g
        g = (g + n//g) // 2
    while g * g > n:
        g -= 1
    return g

You should be able to take it from there.

Tom Zych
  • 13,329
  • 9
  • 36
  • 53