CODEWITHSUNDEEP
pythonnumpyfloating-pointprecisionfloating-accuracy
spore234
spore234

Reputation: 3640

numpy precision with large numbers

I want to factorize a large number using Fermat's factorization method. This is how I implemented it:

import numpy as np

def fac(n):
    x = np.ceil(np.sqrt(n))
    y = x*x - n
    while not np.sqrt(y).is_integer():
        x += 1
        y = x*x - n
    return(x + np.sqrt(y), x - np.sqrt(y))

Using this method I want to factor N into its components. Note that N=p*q, where p and q are prime.

I chose the following values to compute N:

p = 34058934059834598495823984675767545695711020949846845989934523432842834738974239847294083409583495898523872347284789757987987387543533846141.0
q = 34058934059834598495823984675767545695711020949846845989934523432842834738974239847294083409583495898523872347284789757987987387543533845933.0

and defined N

N = p*q

Now I factor N:

r = fac(n)

However, the factorization seems to not be correct:

int(r[0])*int(r[1]) == N

It does work for smaller ints:

fac(65537)
Out[1]: (65537.0, 1.0)

I'm quite sure the reason is numerical precision at some point.

I tried calculating N in numpy using object types:

N = np.dot(np.array(p).astype(object), np.array(q).astype(object))

but it doesn't help. Still, the numpy requires a float for the sqrt function.

I also tried using the math library instead of numpy, this library seems to not require a float for its sqrt function, but ultimately running into precision issues as well.

Upvotes: 0

Views: 1231

Answers (1)

Serge Ballesta
Serge Ballesta

Reputation: 148890

Python int are multiple precision numbers. But numpy is a wrapper around C low level libraries to speed up operations. The downside is that it cannot handle those multi-precision numbers. Worse, if you try to use np.sqrt on them, they will be converted to floating point numbers (C double or numpy float64) what have a precision of about 15 decimal digits.

But as Python int type is already a multiprecision type, you could use math.sqrt to get an approximative value of the true square root, and then use Newton to find a closer value:

def isqrt(n):
    x = int(math.sqrt(n))
    old = None
    while True:
        d = (n - x * x) // (2 * x)
        if d == 0: break
        if d == 1:                 # infinite loop prevention
            if old is None:
                old = 1
            else: break
        x += d
    return x

Using it, your fac function could become:

def fac(n):
    x = isqrt(n)
    if x*x < n: x += 1
    y = x * x - n
    while True:
        z = isqrt(y)
        if z*z == y: break
        x += 1
        y = x*x -n
    return x+z, x-z

Demo:

p = 34058934059834598495823984675767545695711020949846845989934523432842834738974239847294083409583495898523872347284789757987987387543533846141
q = 34058934059834598495823984675767545695711020949846845989934523432842834738974239847294083409583495898523872347284789757987987387543533845933
N = p*q
print(fac(N) == (p,q))

prints as expected True

Upvotes: 3

Related Questions