2

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.

Oliver Jack Baynham
  • 73
  • 5
  • 1
    This probably stems from the fact that `complex` uses two double precision numbers internally for the real and imag part. You could implement your own complex numbers using [`decimal`](https://docs.python.org/2/library/decimal.html), or `round` your result to an appropriate number of decimal places. – Jan Christoph Terasa May 08 '17 at 15:18
  • Ugh, what is `if IsGaussian(n / complex(x, y)) == n:` supposed to do? It _can_ only be True for n = 0 or n = 1. You probably want to get rid of the `== n`. – Jan Christoph Terasa May 08 '17 at 15:44

2 Answers2

3

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)
Efron Licht
  • 488
  • 3
  • 13
1

You could define an epsilon, by using round to round to the desired number of decimal places/precision (e.g. 10):

def IsGaussian(z, prec=10):
    # returns True if the complex number is a Gaussian integer
    # rounds the input number to the `prec` number of digits
    z = complex(round(z.real,prec), round(z.imag,prec))
    return complex(int(z.real), int(z.imag)) == z

Your code has another issue though:

if IsGaussian(n / complex(x, y)) == n:

This will only give results for n = 0 or n = 1. You probably want to remove the check for equality.

Jan Christoph Terasa
  • 5,781
  • 24
  • 34