Fast solutions for $x^2+y^2+z^2=d^2$ for large $d>1000$

Question

I need a fast algorithm to calculate solutions to $x^2+y^2+z^2=d^2$ for $d>1000$.

I tried Pythagorean quadruple (Alternate parametrization approach) for $d<400$ and it isn't that fast.

UPDATE

def EfficientDivisors(n):
    factors = []
    i = 2
    while (i < sqrt(n)):
        if (n%i == 0):
            factors.append(i)
            if (i != (n/i)):
                factors.append(n//i)
        i += 1
    return factors
#alternative parameterization: a and b are both even
def fncTest4(d):
    c = 0
    nSolutionNumber = 0
    for a in range(2, d, 2):
        l = a//2
        for b in range(2, d, 2):
            m = b//2
            nTemp1 = l ** 2 + m ** 2
            factors = EfficientDivisors(nTemp1)
            for n in factors:
                if n ** 2 < nTemp1:
                    nTemp2 = (nTemp1 + n ** 2) // n
                    if d == nTemp2:
                        c = d - 2 * n
                        nSolutionNumber += 1
                        print(nSolutionNumber, ":", a, b, c, d)
                        break

fncTest4(2000) # 4 minutes

In terms of time-complexity I see an $O(N^3)$ solution easily via three nested loops. How fast are you looking for? — unseen_rider, Feb 02 '18 at 16:19
It may may easier to find a pre-computed table for solutions $x, y, z$ and $d$, then use it as a lookup table for finding corresponding values of $d$. I believe there would be pre-computed tables online — unseen_rider, Feb 02 '18 at 16:21
Value of $d$ is sufficiently large, for example, $10^4$. I don't know if a pre-computed table exists — yW0K5o, Feb 02 '18 at 16:34
@yW0K5o I would call $10^4$ rather small, but if you do need to work with large numbers then there is a reference here for a conjecturally polynomial-time (meaning $O(\log^m d)$, not $O(d^m)$) algorithm: https://mathoverflow.net/a/110277 — Erick Wong, Feb 02 '18 at 16:59
try https://en.wikipedia.org/wiki/Pythagorean_quadruple#Parametrization_of_primitive_quadruples — Will Jagy, Feb 02 '18 at 19:12
@yW0K5o It doesn't look inherently slow to me. But there is enough generality in the description that you could have a very slow implementation of the concept. — Erick Wong, Feb 03 '18 at 05:24
I assume you want $(a,b,c)>0$. Are you after solutions or primitive solutions? — Old Peter, Feb 05 '18 at 16:40
@OldPeter Yes, I want the primitive unique solutions. $(a,b,c)>0$ . Then I use permutations to get all solutions (order with signs). — yW0K5o, Feb 05 '18 at 18:42
I've just ran $d=1001$ to $2000$ in under 60 seconds, producing $168505$ solutions, using VB6. If this is fast enough, I'll post my method tomorrow. — Old Peter, Feb 05 '18 at 19:53
OK. It is fast enough. Looking forward to see your solution tomorrow. Thank you !!! — yW0K5o, Feb 05 '18 at 20:05
@yW0K5o I think your simplest solution can be found in another post about matching sides of Pythagorean triples. — poetasis, May 27 '19 at 15:36

Old Peter · Accepted Answer · 2018-02-08T12:03:06.623

1

Here’s the outline.

WLOG, $a>=b>=c$.

Then $a$ has a maximum value of $d-1$ and a minimum value of $\sqrt{\frac{d^2}{3}}$, rounded down.

For each $a$, calculate $u=d^2-a^2$

Then maximum $b$ is $\sqrt{u}$ rounded up, but also $<=a$ and minimum $b$ is $\sqrt{\frac{u}{2}}$, rounded down.

For each $b$, calculate $v=u-b^2$

When $v$ is a square, $c=\sqrt{v}$ and we have a solution.

Good luck with your project. Please let me know if you need any more details.

I’ll add the VB6 code, if I ever find out how to do it.

Some test results, $(c,b,a,d)$

$$(20,28,29,45)$$ $$(15,30,30,45)$$ $$(6,30,33,45)$$ $$(20,20,35,45)$$ $$(4,28,35,45)$$ $$(16,20,37,45)$$ $$(13,16,40,45)$$ $$(8,19,40,45)$$ $$(5,20,40,45)$$ $$(6,15,42,45)$$ $$(5,8,44,45)$$ $$(22,31,34,51)$$ $$(17,34,34,51)$$ $$(24,27,36,51)$$ $$(3,36,36,51)$$ $$(14,31,38,51)$$ $$(1,34,38,51)$$ $$(14,17,46,51)$$ $$(1,22,46,51)$$ $$(14,14,47,51)$$ $$(10,10,49,51)$$ $$(2,14,49,51)$$ $$(1,10,50,51)$$ $$(480,600,640,1000)$$ $$(192,640,744,1000)$$ $$(424,480,768,1000)$$ $$(280,576,768,1000)$$ $$(224,600,768,1000)$$ $$(24,640,768,1000)$$ $$(360,480,800,1000)$$ $$(168,576,800,1000)$$ $$(192,480,856,1000)$$ $$(352,360,864,1000)$$ $$(152,480,864,1000)$$ $$(96,480,872,1000)$$ $$(96,360,928,1000)$$ $$(168,224,960,1000)$$ $$(960,1200,1280,2000)$$ $$(384,1280,1488,2000)$$ $$(848,960,1536,2000)$$ $$(560,1152,1536,2000)$$ $$(448,1200,1536,2000)$$ $$(48,1280,1536,2000)$$ $$(720,960,1600,2000)$$ $$(336,1152,1600,2000)$$ $$(384,960,1712,2000)$$ $$(704,720,1728,2000)$$ $$(304,960,1728,2000)$$ $$(192,960,1744,2000)$$ $$(192,720,1856,2000)$$ $$(336,448,1920,2000)$$

edited Feb 08 '18 at 12:03

answered Feb 06 '18 at 14:31

Old Peter

1,357

Peter, how could you verify if your method gives the same solutions which described in Wikipedia[Pythagorian quadruple]? – yW0K5o Feb 06 '18 at 18:46
I've tested small values against those independently produced on a spreadsheet, and I've used this general method for years. If you give me a value for $d<46340$ I'll post the values for you to check against your results. Please will you tell me how to add my code, as per the grey box in your question? – Old Peter Feb 06 '18 at 19:14
With $d<30$ my code produces the 31 solutions shown in your link, plus 18 others with $gcd(a,b,c)>1$ – Old Peter Feb 06 '18 at 19:48
Peter, I use Chrome browser. When you edit your answer put/click your VB6 code in curly braces. – yW0K5o Feb 06 '18 at 20:10
To verify your solutions use d=50 it gives 2 unique solutions, and canbe verified in this link https://www.wolframalpha.com/input/?i=PowersRepresentations%5B50%5E2,3,2%5D [discard 0-solutions] – yW0K5o Feb 06 '18 at 20:16
Here's the image how to add your code: https://image.ibb.co/h6AOox/code.jpg – yW0K5o Feb 06 '18 at 20:21
Program produces $(c,b,a)=(24,30,32),(18,24,40) for $d=50$. It's easy to change if you want the results with zero given by your link. – Old Peter Feb 07 '18 at 11:56
Try other value for test 45, 51, 1000, 2000. – yW0K5o Feb 07 '18 at 14:56
I implemented your method and it gives less solutions than Wikepedia's one. Could you also provide proof why your method can find all solutions? – yW0K5o Feb 07 '18 at 15:01
Test results added to answer. Please give examples of missing results. I'm still unable to add my code despite many attempts; the code could help. I'll look at a proof. – Old Peter Feb 07 '18 at 15:25
Test for 45: My code: [ 4 28 35 ] [ 6 30 33 ] [ 6 15 42 ] [ 8 19 40 ] [ 5 8 44 ] [ 16 20 37 ] [ 13 16 40 ] [ 20 20 35 ] [ 20 28 29 ] [ 5 20 40 ] [ 15 30 30 ] – yW0K5o Feb 08 '18 at 16:13
Test for 45: Peter's code: def fncTestPeterMethod(d): c = 0 dSquared = d ** 2 nAUpperRange = d - 1 nALowerRange = int(sqrt(dSquared/3)) for a in range(nALowerRange, nAUpperRange, 1): u = dSquared - a**2 nBUpperRange = int(sqrt(u)) nBLowerRange = int(sqrt(u/2)) for b in range(nBLowerRange, nBUpperRange, 1): v = u - b**2 nTemp1 = CheckThePefectSquare(v) if nTemp1 != -1: c = nTemp1 print(a, b, c) – yW0K5o Feb 08 '18 at 16:14
Test for 45:Peter's code results: [ 20 28 29 ][ 15 30 30 ][ 20 20 35 ][ 16 20 37 ][ 13 16 40 ][ 8 19 40 ][ 6 15 42 ] – yW0K5o Feb 08 '18 at 16:14
I fixed my errors it seems like your method works ! – yW0K5o Feb 08 '18 at 16:32
Yes, I've used this method for more than ten years, and would be devastated to find it flawed. Well done for correcting your code. Perhaps you would care to up-vote my answer? – Old Peter Feb 08 '18 at 16:43

Fast solutions for $x^2+y^2+z^2=d^2$ for large $d>1000$

1 Answers1