Pythagorean quadruples
One form of Pythagorean quadruples is (for positive integers a, b, c, and d):
- a2 + b2 + c2 = d2
- Task
For positive integers up 2,200 (inclusive), for all values of a,
b, and c,
find (and show here) those values of d that can't be represented.
Show the values of d on one line of output (optionally with a title).
REXX
brute force
This version is a brute force algorithm, with some optimization (to save compute time)
which pre-computes the squares of the positive integers used in the search.
<lang rexx>/*REXX program computes/displays integers that aren't possible for: a² + b² + c² = d² */
parse arg hi . /*obtain optional argument from the CL.*/
if hi== | hi=="," then hi=2200 /*Not specified? Then use the default.*/
@.=. /*array of integers to be squared. */
!.=. /* " " " squared. */
do j=1 for hi**2 + (hi-2)**2 + (hi-3)**2 /*precompute possible squares (to max).*/
jj=j*j; !.jj=j; if j<=hi then @.j=jj /*define a square; D value; squared # */
end /*j*/
d.=. /*array of possible solutions (D) */
do a=1 for hi-2 /*go hunting for solutions to equation.*/
do b=a to hi-1; ab = @.a + @.b /*calculate sum of 2 (A,B) squares.*/
do c=b to hi; abc= ab + @.c /* " " " 3 (A,B,C) " */
if !.abc==. then iterate /*Not a (perfect) square? Then skip it.*/
s=!.abc; d.s= /*define this D solution as being found*/
end /*c*/
end /*b*/
end /*a*/
say say 'Not possible positive integers for d ≤' hi " using equation: a² + b² + c² = d²" say $= /* [↓] find all the "not possibles". */
do p=1 for hi; if d.p==. then $=$ p /*Not possible? Then add it to the list*/
end /*p*/ /* [↓] display list of not-possibles. */
say substr($, 2) /*stick a fork in it, we're all done. */</lang>
- output when using the default input:
Not possible positive integers for d ≤ 2200 using equation: a² + b² + c² = d² 1 2 4 5 8 10 16 20 32 40 64 80 128 160 256 320 512 640 1024 1280 2048
optimized
This REXX version is an optimized version, it solves the formula:
- a2 + b2 = d2 - c2
This REXX version is over thirty times faster then the previous version. <lang rexx>/*REXX program computes/displays integers that aren't possible for: a² + b² + c² = d² */ parse arg hi . /*obtain optional argument from the CL.*/ if hi== | hi=="," then hi=2200 /*Not specified? Then use the default.*/ high=hi * 3 /*D can be three times the HI (max).*/ @.=. /*array of integers to be squared. */ !.=. /* " " " squared. */ dbs.=. /* differences between squares.*/
do s=1 for high /*process up to three times the maximum*/
ss=s*s; r.ss=s /*precompute squares and square roots. */
if s<=hi then @.s=ss /*only define @.s if within range.*/
end /*s*/
do c=1 for high; cc=c*c /*precompute possible differences. */
do d=c+1 to high; dif=d*d - cc /*process D squared; calc differences*/
if dbs.dif==. then do; dbs.dif=cc; iterate; end /*No list yet? Add it.*/
if wordpos(cc,dbs.dif)==0 then dbs.dif=dbs.dif cc /*Not in list? Add it.*/
end /*d*/
end /*c*/
d.=. /*array of possible solutions (D) */
do a=1 for hi-2 /*go hunting for solutions to equation.*/
do b=a to hi-1; ab=@.a + @.b /*calculate sum of 2 (A,B) squares.*/
if dbs.ab==. then iterate /*Not a difference? Then ignore it. */
do n=1 for words(dbs.ab) /*handle all ints that satisfy equation*/
x=word(dbs.ab, n) /*get a C integer from the DBS.AB list.*/
abc=ab + x /*add the C² integer to A² + B² */
_=r.abc /*retrieve the square root of C² */
d._= /*mark the D integer as being found. */
end /*n*/
end /*b*/
end /*a*/
say say 'Not possible positive integers for d ≤' hi " using equation: a² + b² + c² = d²" say $= /* [↓] find all the "not possibles". */
do p=1 for hi; if d.p==. then $=$ p /*Not possible? Then add it to the list*/
end /*p*/ /* [↓] display list of not-possibles. */
say substr($, 2) /*stick a fork in it, we're all done. */</lang>
- output is the same as the 1st REXX version.