Faster Sage Code for Diophantine Equation? 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

 A: 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.

  
*
  
*Find an integer z$\color{blue}{^{[2]}}$ such that $0 < z < p/2$ and
  $$z^2 + 1 \equiv 0 \pmod p$$
  
*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$.
A: 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=2m+1
        y=2n
        w=(x^2)+(y^2)
        if is_prime(w)==True and w<i:
            solution=(x,y,w)
            primeswithsolutions.append(solution)
print primeswithsolutions

