Faster Sage Code for Diophantine Equation?

Question

I'm having trouble with the computation time. Does anyone have any ideas for faster code?

#prints a list of tuples containing the solutions to x^2+y^2=p where p is a prime less than a million
i=10^6
primeswithsolutions=[]
for x in range(i):
    for y in range(i):
        w=(x^2)+(y^2)
        if is_prime(w)==True and w<i:
            solution=(x,y,w)
            primeswithsolutions.append(solution)
print primeswithsolutions

Thank you for any help you can provide.

---

Answered: Realized it was just the range that was messed up.

Thank you!

#prints a list of tuples containing the solutions to x^2+y^2=p where p is a prime less than a million
i=10^3
primeswithsolutions=[(1, 1, 2)]
for x in range(i):
    if is_even(x)==True:
        for y in range(x,i):
            if is_odd(y)==True:
                w=(x^2)+(y^2)
                if is_prime(w)==True and w<(i^2):
                    solution=(x,y,w)
                    primeswithsolutions.append(solution)
print primeswithsolutions

Answer 1

Putting those comments together:

prints a list of tuples containing the solutions to x^2+y^2=p where p is a prime less than a million

i=10^6
primeswithsolutions=[]
for m in range(i/2):
    for n in range(m,i/2):
        x=2*m+1
        y=2*n
        w=(x^2)+(y^2)
        if is_prime(w)==True and w<i:
            solution=(x,y,w)
            primeswithsolutions.append(solution)
print(primeswithsolutions)

For example $i=100$, it makes

[(1, 2, 5), (1, 4, 17), (1, 6, 37), (3, 2, 13), (3, 8, 73), (5, 4, 41), (5, 6, 61), (5, 8, 89)]

Answer 2

It seems the real problem of your code lies in the range of $x$ and $y$. Notice $$x^2 + y^2 = p < i = 10^6 \quad\implies x, y < \sqrt{i} = 10^3$$

The range to use should be something like $\text{range}(\sqrt{i})$ instead of $\text{range}(i)$.
( Not sure about the correct SAGE syntax).

In any event, instead of looping over $x$ and $y$, an alternate approach is loop over the set of primes $p$ of the form $4k+1$ and uses following algorithm$\color{blue}{^{[1]}}$ to express $p$ as a sum of squares.

1. Find an integer z$\color{blue}{^{[2]}}$ such that $0 < z < p/2$ and $$z^2 + 1 \equiv 0 \pmod p$$

2. Carry out Euclidean algorithm on $p/z$ (not $z/p$), producing a sequence of remainders $R_1, R_2, \ldots,$ to the point where $R_k$ is first less than $\sqrt{p}$, then

$$p = \begin{cases} R_{k}^2 + R_{k+1}^2,&\text{if} R_1 > 1\\ z^2 + 1,&\text{if} R_1 = 1 \end{cases}$$

For more information, look at the paper cited below and answers for a similar question.

Notes

• $\color{blue}[1]$ Notes on Representing a Prime as a Sum of Two Squares, John Brillhart, Mathematics of Computation, Volume 26, Number 120, October 1972.

  An online copy can be found here.

• $\color{blue}[2]$ If $c$ is any quadratic non-residue for $p = 4k+1$, then $z = c^k$ satisfies $$z^2 + 1 \equiv 0 \pmod p$$ Since half of the numbers $1, 2, \ldots, p-1$ are quadratic non-residues, we can randomly pick a $c$ and test whether $c^k$ works or not. On average, we only need to do this twice to find a working $z$.