Stop Approximation in Complex Division in Python

Question

I've been writing some code to list the Gaussian integer divisors of rational integers in Python. (Relating to Project Euler problem 153)

I seem to have reached some trouble with certain numbers and I believe it's to do with Python approximating the division of complex numbers.

Here is my code for the function:

def IsGaussian(z):
    #returns True if the complex number is a Gaussian integer
    return complex(int(z.real), int(z.imag)) == z

def Divisors(n):
    divisors = []

    #Firstly, append the rational integer divisors
    for x in range(1, int(n / 2 + 1)):
        if n % x == 0:
            divisors.append(x)

    #Secondly, two for loops are used to append the complex Guassian integer divisors
    for x in range(1, int(n / 2 + 1)):
        for y in range(1, int(n / 2 + 1)):
            if IsGaussian(n / complex(x, y)) == n:
                divisors.append(complex(x, y))
                divisors.append(complex(x, -y))

    divisors.append(n)

    return divisors

When I run Divisors(29) I get [1, 29], but this is missing out four other divisors, one of which being (5 + 2j), which can clearly be seen to divide into 29.

On running 29 / complex(5, 2), Python gives (5 - 2.0000000000000004j)

This result is incorrect, as it should be (5 - 2j). Is there any way to somehow bypass Python's approximation? And why is it that this problem has not risen for many other rational integers under 100?

Thanks in advance for your help.

Answer 1

Internally, CPython uses a pair of double-precision floats for complex numbers. The behavior of numerical solutions in general is too complicated to summarize here, but some error is unavoidable in numerical calculations.

EG:

>>>print(.3/3)
0.09999999999999999

As such, it is often correct to use approximate equality rather than actual equality when testing solutions of this kind.

The isclose function in the cmath module is available for this exact reason.

>>>print(.3/3 == .1)
False
>>>print(isclose(.3/3, .1))
True

This kind of question is the domain of Numerical Analysis; this may be a useful tag for further questions on this subject.

Note that it is considered 'pythonic' for function identifiers to be in snake_case.

from cmath import isclose
def is_gaussian(z):
    #returns True if the complex number is a Gaussian integer
    rounded = complex(round(z.real), round(z.imag))
    return isclose(rounded, z)