3

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
Community
  • 1
st1221
  • 31
  • 1
    Perhaps you should give more detail about the math problem being solved. It seems that you are taking all the possible sums of two squares and checking the sum to see if it is prime. Certainly it would seem more effective to consider primes congruent to 1 mod 4, and find all the ways to express as a sum of two squares (if that is really your goal). [NB: Yes, 2=1+1 is a special case.] – hardmath Apr 11 '15 at 01:08
  • Well, you certainly don't need $y$ to go over the entire range; you're testing values like $(3, 5)$ and then $(5, 3)$. So rather than y in range(i), you can have the statement y in range(x, i). You could also do some divisibility tests; only execute the code if $\gcd(x, y) = 1$, otherwise their sum of squares will be divisible by that $\gcd$. (Note: I have no idea if divisibility tests are any quicker than the tests you're running in the loop; it might not save time, but it probably would). – pjs36 Apr 11 '15 at 01:09
  • Also, here's William Stein's Sage for Power Users guide; section 3 about using Cython might be helpful, but I've never attempted it myself. – pjs36 Apr 11 '15 at 01:14
  • @pjs36 Thank you. This save me a lot of time. – st1221 Apr 11 '15 at 01:18
  • 1
    As stated above, there is one special case where both x and y are odd ($1^2+1^2=2$). Apart from that all primes are odd, so if $x$ is odd, $y$ must be even. You can therefore reduce calculation time by restricting $x$ to odd numbers in the range and $y$ to even numbers in the range. – tomi Apr 11 '15 at 01:23
  • @will-jagy I just wanted to write a program that would give solutions. I did so, but figured the run time could be improved considering I've let it run for approximately 30 minutes and still no output. – st1221 Apr 11 '15 at 01:24
  • 1
    Also, using range(n) instead of xrange(n) will suck up a lot of memory. xrange only generates the next number as needed whereas range(n) puts the whole list into memory each time it is called. – TravisJ Apr 11 '15 at 01:31
  • presumably the is_prime() function is something reasonably optimized? That could potentially also be a bottleneck (I'm not a sage expert, but I code a lot in python) – TravisJ Apr 11 '15 at 01:32
  • @TravisJ Yes, is_prime() is optimized. I've found that running most code with range 10^4+ is pretty slow in sage in most cases. – st1221 Apr 11 '15 at 01:53
  • A better place for such questions is Ask SageMath at http://ask.sagemath.org/ – Samuel Lelièvre Apr 11 '15 at 11:13

2 Answers2

1

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$.

Community
  • 1
achille hui
  • 122,701
0

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)]
user682141
  • 836
tomi
  • 9,594