# Solving $x^2 + y^2 = 1 + z^4$ with (x,y,z) = 1 and z < x < y

I have a computer programming problem where I need to find n many sets of integers that meet the condition $x^2 + y^2 = 1 + z^4$ with (x,y,z) = 1 and z < x < y

I can do this relatively easily with a brute force algorithm that increments z, then finds an x which has a GCD of 1 with z, then finds a y that has a GCD of 1 with z & x. Which then checks to see if the co-prime tripple of x, y & z satisfy the above equation.

However the program is unbearably slow making millions of function calls once I've discovered the first 5 - 10 sets.

How can I restate the problem to quickly find a candidate x & y given some z? Is there some other approach I could consider over this slow brute force method?

• Possible hint: rewrite the equation to use difference of squares. – user2469 Jan 17 '16 at 3:08
• You can first solve the equation. $x^2+y^2=z^2+1$ Then make choices when there are $z -$ square. The formula there. math.stackexchange.com/questions/351491/… math.stackexchange.com/questions/74931/… – individ Jan 17 '16 at 4:34
• @barrycarter So: $x^2 = 1 + (z^2)^2 - y^2$ ... or $x^2 = 1 + (z^2 + y)(z^2 - y)$ ... ? I appreciate the hint but I'm still lost here. – Peanutter Jan 17 '16 at 4:40
• @individ Thanks for the links you provided. I read through both answers on those links but don't understand them. Unfortunately you may be overestimating my maths knowledge and experience. Reading through the wiki page on Pell's equation was no help either. – Peanutter Jan 17 '16 at 4:49
• x^2 - 1 = z^4 - y^2 is a difference of square on both sides. Not sure if this helps, just a random thought. – user2469 Jan 17 '16 at 5:40

You don't have to worry about $(x,y,z)=1$ because it's automatic: if $x\equiv y\equiv z\equiv 0 \pmod{p}$ then your condition $x^2+y^2=z^4+1$ becomes $0\equiv1\pmod{p}$.

Edit: Modified code to cater for missed constraint $z<x$ as per @alex.jordan's comment.

This brute force python fragment produces about 380 solutions in a few seconds on my laptop:

for z in range(1,300+1):
target = z*z*z*z+1
for x in range(z+1,int((target/2)**.5)+1):
y = int((target-x*x)**.5)
if x*x + y*y == target:
print("Solution: z=",z,"x=",x,"y=",y)


(Here (blah)**.5 is python for $\sqrt{\text{blah}}$.) And int(blah) is floor(blah).

Basically, for a given $z$ and viable $x$ values, you test the single plausible $y$ value of $\left\lfloor\sqrt{z^4-1-x^2}\right\rfloor$.

If you need to find even more solutions efficiently, I think you'll have to factor $z^4+1$ into primes of the form $4k+1$ (plus an optional factor of 2), decompose the primes into sums of squares, then recombine those solutions. It'd be faster, but take a lot more code and you'd still be limited by the size of the numbers you can factor which probably won't go much beyond the $300^4+1$ in the code fragment above.

Incidentally, there will always be at least one solution to $x^2+y^2=z^4+1$ for each $z$ since $z^4+1$ is a product of primes of the form $4k+1$, plus an optional factor of 2. For a prime $p$ of the form $4k+3$, $z^4+1\equiv0 \pmod{p}$ has no solutions since $x^2\equiv-1\pmod{p}$ has no solutions. And $z^4+1\equiv0\pmod{4}$ has no solutions since $x^2\equiv3\pmod{4}$ has no solutions, so there is at most 1 factor of 2 in $z^4+1$. Any number can be written as the sum of two squares iff is it the product of powers of 2 and $4k+1$ type primes and even powers of $4k+3$ primes. The constraint $z<x$ doesn't always hold in the solutions, though (thanks to alex for remark in comment below).

• Regarding the last paragraph, OP is requiring $z<x<y$, so there is not always such a solution for each $z$. For instance with $z=2$, one of $x,y$ must equal $\pm1$. Is the condition $z<x<y$ built in to your search algorithm? – alex.jordan Jan 17 '16 at 7:03
• @alex.jordan: You're absolutely correct. I'd missed that completely. Will update the code and comments! – Frentos Jan 17 '16 at 7:08
• $x$ can start at $z+1$ and need only go up to $\left\lfloor\sqrt{\frac{1+z^4}{2}}\right\rfloor$. – alex.jordan Jan 17 '16 at 7:10
• For the most part there will be quite a bit of flexibility in the presentation as a sum of two squares, at least if there enough prime factors. So I am not overly concerned about finding $x,y$ pairs that are both larger than $z$. I guess that prime factorization will prove to be a bottleneck when you go for larger choicec of $z$ though. +1 for the observations that 1) checking gcd's is a waste of effort,2) using prime factorizations – Jyrki Lahtonen Jan 17 '16 at 7:11
• Great answer, thank you. – Peanutter Jan 17 '16 at 9:25

Write for the equation.

$$x^2+y^2=z^4+1$$

One simple formula.

$$x=(\frac{a^4+1}{2}t+2a)ta^2+1$$

$$y=\frac{a^8-1}{4}t^2+a(a^4-1)t+a^2$$

$$z=\frac{a^4+1}{2}t+a$$

If you use the solutions of the equation Pell. Where $t -$ ask yourself.

$$p^2-2(t^2+1)s^2=1$$

Make the change. $$a=2t(p+ts)s$$

Then decisions can be recorded.

$$x=(2a+1)(t^2+1)-1$$

$$y=2a(a+1)(t^2+1)-2a-1$$

$$z=t(p^2+2tps+2(t^2+1)s^2)$$

• Nice (FWIW I verified this)! Where does this come from? You do need either $a$ odd or $t$ even to get integers, and $a>1$ to get $x<y$. – Jyrki Lahtonen Jan 19 '16 at 11:23

Here is a quick and dirty method for generating infinitely many solutions. The idea is very crude in the sense that it only produces solutions, where $y=z^2-3$. Restricting ourselves in this way reduces the problem to a Pell equation, and we can solve those efficiently. My write-up does not assume any knowledge of Pell equations, but those familiar with their theory will be bored.

Given all this the algorithm below only finds a tiny fraction of all the solutions of your equation, but it does find this subset very quickly.

If plug in $y=z^2-3$ we get the following equation $$x^2-6z^2=-8.\qquad(*)$$ So the question is how to find $(x,z)$ pairs satisfying this. To that end let me introduce the following functions:

• The function $N:\Bbb{Z}^2\to\Bbb{Z}$ that takes two integers as inputs and spews out a single integer given by $$N(a,b)=a^2-6b^2$$
• The function $M:(\Bbb{Z}^2)^2\to\Bbb{Z}^2$ that take two pairs of integers as inputs and outputs one pair of integers. The formula is $$M((a,b);(c,d))=(ac+6bd,ad+bc).$$

Remark: The function $N$ gives the norm of the element $a+b\sqrt6$ in the ring $\Bbb{Z}[\sqrt6]$, and the function $M$ calculates the product of the numbers $(a+b\sqrt6)(c+d\sqrt6)$ in that ring.

Key Lemma. We have for all integers $a,b,c,d$ the identity $$N(M((a,b);(c,d)))=N(a,b) N(c,d).$$

Proof. You do this. Or you can just take my word for it.

The Key Lemma allows us to find solutions to $(*)$, i.e. $N(x,z)=-8$, if we can locate one solution to it, and if we can also locate a solution to the equation $N(a,b)=1$. These are easy to find with a bit of searching among small integers. Namely $$N(4,2)=4^2-6\cdot2^2=-8$$ and $$N(5,2)=5^2-6\cdot2^2=1.$$

Here is then an algorithm for producing infinitely many solutions to $(*)$. We first produce a sequence of pairs $(a_n,b_n)$, $n=1,2,\ldots$, such that $N(a_n,b_n)=1$:

1. The starting point is $(a_1,b_1)=(5,2)$, i.e. the solution found above.
2. Given a solution $(a_k,b_k)$ we then find another solution $(a_{k+1},b_{k+1})$ with the recurrence formula $$(a_{k+1},b_{k+1})=M((a_k,b_k),(5,2))=(5a_k+12b_k,5b_k+2a_k).$$

The point is that the Key Lemma implies that if $a^2-6b^2=1$ and $c^2-6d^2=1$, then also $(ac+6bd)^2-6(ad+bc)^2=1$.

On with producing solutions to $(*)$. We saw that $(x_0,z_0)=(4,2)$ is one solution. We can then use the pairs produced above to get others by setting $$(x_k,z_k)=M((x_0,z_0);(a_k,b_k))=(4a_k+12b_k,2a_k+4b_k)$$ for $k=1,2,\ldots$. Observe that $(*)$ implies $x\approx z\sqrt6$, so the inequality $x_k>z_k$ is more or less given. Because $y_k=z_k^2-3$, the other requirement $y_k>x_k$ is also more or less clear, at least for all sufficiently large values of $k$. A quadratic polynomial grows faster than a linear one.

Here's how the list of solutions begins. I also present the calculations using $\Bbb{Z}[\sqrt6]$ arithmetic. The numbers will grow very quickly here. That is a price for it being deterministic and fast (no search is involved). You have to judge whether this is a problem.

1. The starting solution to $(*)$, namely $(x_0,z_0)=(4,2)$ must be discarded because here $y_0=z_0^2-3=1$ is too small.
2. Because $$(4+2\sqrt6)(5+2\sqrt6)=44+18\sqrt{6}$$ we get as the next solution $(x_1,z_1)=(44,18)$. This gives $y_1=z_1^2-3=321$ and a valid solution $(x,y,z)=(44,321,18)$.
3. Because $$(44+18\sqrt6)(5+2\sqrt6)=436+178\sqrt6$$ we get $(x_2,z_2)=(436,178)$. Here $y_2=178^2-3=31681$, and we have the solution $(x,y,z)=(436,31681,178)$.

Essentially we keep multiplying by $5+2\sqrt6$ to get to the next solution of this form.

• We see that for all the solutions produced here both $x$ and $z$ will be even. As is evident from Frentos' answer this is not a problem, because $y$ will always be odd, and thus $\gcd(x,y,z)$ can be $1$ (and, in fact, is $=1$ as per Frentos' argument). – Jyrki Lahtonen Jan 17 '16 at 7:49
• It would not surprise me, if there were another method for producing a lot of smaller solutions as well. If you can factor the numbers $z^4+1$ produced here, and there are several factors, you can probably find more $(x,y)$ pairs for the given $z$ using the idea described in the last paragraph of Frentos' answer. – Jyrki Lahtonen Jan 17 '16 at 8:12
• Each iteration here multiplies $x$ and $z$ by roughly ten, and $y$ by a factor of roughly one hundred. Exponential growth cannot be avoided, when we reduce the problem to a specific Pell equation with a choice like $y=z^2-3$. – Jyrki Lahtonen Jan 17 '16 at 8:16