title: Optimizing an algorithm with 3 inner loops for searching 6-tuples of squares using Python
body: <p>The following algorithm generates 6-tuples of squares <code>[s,t,u,t+u,t+u−s,t−s]</code> up to <code>s=13572</code> (see <a href="https://github.com/Sultanow/pythagorean/blob/main/pythagorean_stu.txt" rel="nofollow noreferrer">this Textfile</a> on GitHub). Then the algorithm runs for a week without no result (I stopped it):</p>
  <pre><code>import numpy as np
  from math import sqrt
  from numba import jit
  
  @jit('void(uint64)', forceobj=True)
  def generateData(limit: np.uint64):
      for t in np.arange(17835, limit+1, dtype=np.uint64):
          for s in np.arange(13572, t, dtype=np.uint64):
              for u in np.arange(121220, limit+1, dtype=np.uint64):
                  if t != u:
                      tt = np.uint64(t*t)
                      uu = np.uint64(u*u)
                      t_u = np.uint64(tt+uu)
                      sqr = np.uint64(sqrt(t_u))
                      if sqr*sqr == t_u:
                          ss = np.uint64(s*s)
                          t_u_s = np.uint64(t_u-ss)
                          sqr = np.uint64(sqrt(t_u_s))
                          if sqr*sqr == t_u_s:
                              t_s = np.uint64(tt-ss)
                              sqr = np.uint64(sqrt(t_s))
                              if sqr*sqr == t_s:
                                  print([s, t, u, ss, tt, uu, t_u, t_u_s, t_s])
  
  limit = np.uint64(500000)
  generateData(limit)
  </code></pre>
  <p><em><strong>The motivation behind</strong></em> is <a href="https://math.stackexchange.com/questions/4342283/for-integers-xyz-why-are-these-cases-impossible-for-mengolis-six-square-pr">this question</a>. These 6-tuples form a basis for searching 4 perfect squares whose difference in each is again a square, which is a really interesting problem. The according searching routine is:</p>
  <pre class="lang-py prettyprint-override"><code>import pandas as pd
  from z3 import Ints, solve
  
  df = pd.read_csv('pythagorean_stu.txt', header=None)
  tuples = df.to_numpy()
  x, y, z, w = Ints('x y z w')
  for row in tuples:
      s=int(row[3])
      t=int(row[4])
      u=int(row[5])
      solve(x*x-w*w==s, y*y-w*w==t, z*z-y*y==u, w!=0)
  </code></pre>
  <p><em><strong>My question</strong></em>: Can we optimize this 6-tuple generation using Parallelism, CUDA or JIT? Currently it is possible to resume the algorithm. It would be great to keep this resumability feature.</p>