misc, Pythagorean triples

Just for fun, I have been playing with a QBasic program that generates
"Pythagorean triples" - sets of three integers, A, B, and C, such that
A^2 + B^2 = C^2. Obviously, these could be the lengths of the sides of
a right-angled triangle in which all the lengths are integers.

There are a couple of other "rules of the game". The program must *not*
output triples which are multiples of other valid triples, such as
6,8,10, which is just twice 3,4,5. Also, the triples must be output in
order, sorted by increasing A (the shortest side of the triangle), and
then by increasing B (the mid-length side).

The following version is the latest. It will generate all 1457 triples
in which A<=1000 in just a few seconds. (The limiting value for A is
the number L, which is set at the beginning of the program.) I have run
it as far as A=10,000, but by that point it gets boringly slow. By
then, 19,271 triples have been output.

Anyway... Does anyone here have anything much faster?

dow

----------------------------------------------------

' pythagorean triples
DEFLNG A-Z
L = 1000 ' limiting length of shortest side
DIM B(1 TO SQR(L)), C(1 TO SQR(L))
P = 0
CLS
FOR A = 3 TO L
Z = 0
IF A MOD 2 THEN
FOR D = 1 TO SQR(A) STEP 2
E# = A / D
E = E#
B = (E * E - D * D) \ 2
GOSUB XX
NEXT
END IF
FOR D = 1 TO A \ 2 STEP 2
E# = SQR(A + A + D * D)
E = E#
B = D * E
GOSUB XX
NEXT
FOR S = Z - 1 TO 1 STEP -1
FOR T = 1 TO S
IF B(T) > B(T + 1) THEN
SWAP B(T), B(T + 1)
SWAP C(T), C(T + 1)
END IF
NEXT
NEXT
FOR S = 1 TO Z
P = P + 1
PRINT P, , A, B(S), C(S)
NEXT
NEXT
END
XX:
IF E = E# THEN
IF B > A THEN
IF E MOD 2 THEN
U = B
V = A
DO WHILE V > 1
W = U MOD V
U = V
V = W
LOOP
IF U = 1 OR V = 1 THEN
C = (E * E + D * D) \ 2
Z = Z + 1
B(Z) = B
C(Z) = C
END IF
END IF
END IF
END IF
RETURN

------------------------------------------

Hmmm... An hour or two later, I have this MUCH faster version...

dow

-------------------------------------------------------

' pythagorean triples
DEFLNG A-Z
DEFDBL Q
L = 1000 ' limiting length of shortest side
DIM B(1 TO SQR(L)), C(1 TO SQR(L))
P = 0
CLS
FOR A = 3 TO L
Z = 0
IF A MOD 2 THEN
FOR D = 1 TO SQR(A) STEP 2
E# = A / D
E = E#
B = (E * E - D * D) \ 2
GOSUB XX
NEXT
END IF
FOR T = 1 TO 2 * SQR(A / 2)
Q1 = T / 2
Q2 = A / T
E# = Q1 + Q2
E = E#
D = Q2 - Q1
B = D * E
GOSUB XX
NEXT
FOR S = 1 TO Z
P = P + 1
PRINT P, , A, B(S), C(S)
NEXT
NEXT
END
XX:
IF B > A THEN
IF D MOD 2 AND E MOD 2 THEN
IF ABS(E# - E) < .000000000001# THEN
U = B
V = A
DO WHILE V > 1
W = U MOD V
U = V
V = W
LOOP
IF U = 1 OR V = 1 THEN
C = (E * E + D * D) \ 2
Z = Z + 1
B(Z) = B
C(Z) = C
FOR S = Z TO 2 STEP -1
IF B(S - 1) < B(S) THEN EXIT FOR
SWAP B(S - 1), B(S)
SWAP C(S - 1), C(S)
NEXT
END IF
END IF
END IF
END IF
RETURN

-------------------------------------------------

If anyone is interested, here is an even faster version of my program
that prints "Pythagorean Triples", A, B, C, such that A^2 + B^2 = C^2,
A < B < C, all in "lowest terms", and sorted on A and then on B.

On my machine, it will generate all the triples in which A <= 10000
(there are 19271 of them) in about half a minute. However, the PRINT
instruction has to be disabled. PRINTing adds a lot to the time.

dow

(david.williams(at)NOSPAMablelink.org)

------------------------------------------------

' pythagorean triples
DEFLNG A-Z
L = 10000 ' limiting length of shortest side
P = 0
CLS
FOR A = 3 TO L
IF A MOD 2 THEN
FOR D = INT(SQR(A) - 1) OR 1 TO 1 STEP -2
E# = A / D
E = E#
IF E = E# THEN
B = ((E + D) * (E - D)) \ 2
GOSUB XX
END IF
NEXT
ELSE
FOR Q1 = INT(SQR(A / 2)) TO 1 STEP -1
Q2# = A / (Q1 + Q1)
Q2 = Q2#
IF Q2 = Q2# THEN
E = Q1 + Q2
IF E MOD 2 THEN
D = Q2 - Q1
B = D * E
GOSUB XX
END IF
END IF
NEXT
END IF
NEXT
END
XX:
IF B > A THEN
U = B
V = A
DO WHILE V > 1
W = U MOD V
U = V
V = W
LOOP
IF V THEN
C = (E * E + D * D) \ 2
P = P + 1
PRINT P, , A, B, C
END IF
END IF
RETURN

--------------------------------------------------

Here is the latest version. If the PRINT line is disabled, it takes
only a few seconds to produce all the 19,271 triples in which the
smallest number (the shortest side of the triangle) is less than or
equal to 10,000. PRINTing the answers slows it down a lot, however.

dow

----------------------------------------------

' pythagorean triples
DEFLNG A-Z
Mn = 3 ' minimum length of shortest side
Mx = 10000 ' maximum length of shortest side
N = 0
Z# = SQR(SQR(2) - 1)
CLS
FOR A = Mn TO Mx
SELECT CASE A MOD 4
CASE 0
A2 = A \ 2
FOR F = INT(Z# * SQR(A2)) TO 1 STEP -1
T# = A2 / F
G = T#
IF G = T# THEN
E = F + G
IF E MOD 2 THEN
D = G - F
B = D * E
GOSUB XX
END IF
END IF
NEXT
CASE 2
CASE ELSE
FOR D = INT(Z# * SQR(A) - 1) OR 1 TO 1 STEP -2
T# = A / D
E = T#
IF E = T# THEN
B = ((E + D) * (E - D)) \ 2
GOSUB XX
END IF
NEXT
END SELECT
NEXT
END
XX:
U = B
V = A
DO WHILE V > 1
W = U MOD V
U = V
V = W
LOOP
IF V THEN
C = (E * E + D * D) \ 2
N = N + 1
PRINT N, , A, B, C
END IF
RETURN

--------------------------------------------

I wonder if it might be possible to generate these triples using a
forward difference technique along the following lines.

DEFLNG A-Z
CLS
LET N = 0
LET A1 = 1
LET B0 = 0
LET B1 = 0
DO WHILE A1 < 9999
LET A1 = A1 + 2
LET B0 = B0 + 4
LET B1 = B1 + B0
LET N = N + 1
PRINT N, , A1, B1, B1 + 1
LOOP
END

Of course the above only generates a particular subset of the triples,
rather than all of them and it doesn't remove common factors as your
method does. On the other hand it is extremely fast (a fraction of a
second when it doesn't have to print) and it may be possible to extend
the idea.

Cheers

Derek

The whole point of the 'exercise' is to generate primitive
pythagorean triples, since multiple forms, if needed, can
be trivially computed from these.


Derek's programme implements the standard mathematical
algorithm, i.e. there's really nothing new about it. Even
the variable naming is the same. Renaming of variables and
the subroutine would make it far more readable, IMHO. That
is, if measningful names are allowed in QBasic. Presumably
comments are allowed, too, so one really shouldn't hold
back on them. All these things are an invitation to people
to actually read and analyse the programme and, who knows,
maybe give some feedback. ;)

--
Mail sent to this email address is automatically deleted
(unread) on the server. Send replies to the newsgroup.

Oops, I really meant to say: David's programme ...

Sorry about that, David.

--
Mail sent to this email address is automatically deleted
(unread) on the server. Send replies to the newsgroup.

-> > Here is the latest version. If the PRINT line is disabled, it takes
-> > only a few seconds to produce all the 19,271 triples in which the
-> > smallest number (the shortest side of the triangle) is less than or
-> > equal to 10,000. PRINTing the answers slows it down a lot, however.

-> I wonder if it might be possible to generate these triples using a
-> forward difference technique along the following lines.

-> DEFLNG A-Z
-> CLS
-> LET N = 0
-> LET A1 = 1
-> LET B0 = 0
-> LET B1 = 0
-> DO WHILE A1 < 9999
-> LET A1 = A1 + 2
-> LET B0 = B0 + 4
-> LET B1 = B1 + B0
-> LET N = N + 1
-> PRINT N, , A1, B1, B1 + 1
-> LOOP
-> END

-> Of course the above only generates a particular subset of the triples,
-> rather than all of them and it doesn't remove common factors as your
-> method does. On the other hand it is extremely fast (a fraction of a
-> second when it doesn't have to print) and it may be possible to extend
-> the idea.

-> Cheers

-> Derek

Sure, it's possible to generate subsets of the triples much faster than
the complete set. It is also possible to generate triples extremely
fast by ignoring their order, e.g. just zipping through lots of values
of D and E and sticking them into the formulae:

A = D * E
B = (E^2 - D^2) / 2
C = (E^2 + D^2) / 2

If D and E are odd positive integers, and E>D, these formulae will
always produce a valid triple, though not necessarily in lowest terms.

But the rules of the game say that the triples must be in order, in
lowest terms, and with A < B < C, in the order in which they are
printed.

I keep fiddling with the program. Here's the latest version:

dow

------------------------------------------------------

' pythagorean triples
DEFLNG A-Z
Mn = 3 ' minimum length of shortest side
Mx = 10000 ' maximum length of shortest side
N = 0
Z# = SQR(2) - 1
CLS
FOR A = Mn TO Mx
SELECT CASE A MOD 4
CASE 0
A2 = A \ 2
FOR F = INT(SQR(Z# * A2)) TO 1 STEP -1
G = A2 \ F
IF G * F = A2 THEN
E = F + G
IF E MOD 2 THEN
D = G - F
B = D * E
GOSUB XX
END IF
END IF
NEXT
CASE IS <> 2
FOR D = INT(SQR(Z# * A)) - 1 OR 1 TO 1 STEP -2
E = A \ D
IF E * D = A THEN
B = ((E + D) * (E - D)) \ 2
GOSUB XX
END IF
NEXT
END SELECT
NEXT
END
XX:
U = B
V = A
DO WHILE V > 1
W = U MOD V
U = V
V = W
LOOP
IF V THEN
C = (E * E + D * D) \ 2
N = N + 1
PRINT N, , A, B, C
END IF
RETURN

-------------------------------------------------

-> > Derek's programme ...

-> Oops, I really meant to say: David's programme ...

-> Sorry about that, David.

I've been called a lot worse!

No problem!

dow

Here's another version, somewhat faster than the last one:

-------------------------------------------------------

' pythagorean triples
DEFLNG A-Z
Mn = 3 ' minimum length of shortest side
Mx = 10000 ' maximum length of shortest side
N = 0
Z# = SQR(2) - 1
CLS
FOR A = Mn TO Mx
T = A MOD 4
SELECT CASE T
CASE 0
S = A \ 2
FOR F = INT(SQR(Z# * S)) TO 1 STEP -1
G = S \ F
IF G * F = S THEN
IF (F XOR G) AND 1 THEN
E = F + G
D = G - F
GOSUB XX
END IF
END IF
NEXT
CASE IS <> 2
S = INT(SQR(Z# * A))
FOR D = S - 1 OR 1 TO 1 STEP -2
E = A \ D
IF E * D = A THEN GOSUB XX
NEXT
END SELECT
NEXT
END
XX:
U = E
V = D
DO WHILE V > 1
W = U MOD V
U = V
V = W
LOOP
IF V THEN
IF T THEN B = ((E + D) * (E - D)) \ 2 ELSE B = D * E
C = (E * E + D * D) \ 2
N = N + 1
PRINT N, , A, B, C
END IF
RETURN

-------------------------------------------------

While I'm at it, here's another little program that I found on Fidonet.
It generates "near square" Pythagorean triples, of the form A, B, C,
where B = A + 1, and C is the hypotenuse of the triangle. If the
numbers are large (the program generates triples up to 16 digits long),
these triangles resemble half a square, divided across a diagonal.

I haven't a clue how this thing works. Would anyone care to figure it
out?

-----------------------------------------------------

DEFDBL A-Z
b0 = 0: c0 = 1
b1 = 1: c1 = 1
FOR n = 1 TO 20
b2 = b1 * 6 - b0 - 2
c2 = c1 * 6 - c0
PRINT b2 - 1; TAB(25); b2; TAB(50); c2
b0 = b1: b1 = b2: c0 = c1: c1 = c2
NEXT

-------------------------------------------------

dow

This looks like the same method of generating a sequence of triples as
I posted earlier, except that it's been applied to the "most nearly
equilateral triangle" subset instead of the "most extremely pointy
triangle" subset. That's why the formulae inside the loop are
different. It's based on Newton's difference method of producing
functions for sequences which would require quite a long post to
describe fully.

Cheers

Derek

David Williams" < XXXX@XXXXX.COM > wrote in message
news: XXXX@XXXXX.COM ...




It was surprisingly easy to convert this program to XBasic - only 6 changes
were made

[1] add a prolog containing a PROGRAM and END PROGRAM statement, an IMPORT
statement to enable square root and a declaration for the main function
Entry()

[2] FUNCTION and END FUNCTION statements for Entry()

[3] Comment out DEFLNG A-Z(not needed), END (not needed) and CLS (not used
by XB)

[4] Change XX: to SUB XX and RETURN to END SUB

[5] Change 2 to 2# and SQR() to SQRT()

[6] Change IS <> 2 to 1,3


Here is the programme in XBasic

PROGRAM "Pythagorean Triples"
IMPORT "xma"
DECLARE FUNCTION VOID Entry()
FUNCTION VOID Entry()
'DEFLNG A-Z ' not needed
Mn = 3 ' minimum length of shortest side
Mx = 10000 ' maximum length of shortest side
N = 0
Z# = SQRT(2#) - 1 ' 2# and SQRT
'CLS ' not used

FOR A = Mn TO Mx
T = A MOD 4
SELECT CASE T
CASE 0
S = A \ 2
FOR F = INT(SQRT(Z# * S)) TO 1 STEP -1
G = S \ F
IF G * F = S THEN
IF (F XOR G) AND 1 THEN
E = F + G
D = G - F
GOSUB XX
END IF
END IF
NEXT F
CASE 1,3 ' IS <> 2
S = INT(SQRT(Z# * A))
FOR D = (S - 1) OR 1 TO 1 STEP -2
E = A \ D
IF E * D = A THEN GOSUB XX
NEXT
END SELECT
NEXT A
'END 'not needed

SUB XX 'XX:
U = E
V = D

DO WHILE V > 1
W = U MOD V
U = V
V = W
LOOP

IF V THEN
IF T THEN B = ((E + D) * (E - D)) \ 2 ELSE B = D * E
C = (E * E + D * D) \ 2
N = N + 1
' PRINT N, , A, B, C
END IF
END SUB 'RETURN
END FUNCTION
END PROGRAM


For Mx=10,000 this runs too fast to measure even in interpreted mode.
Changing the upper limit to 100,000 results in a run time of 220
milliseconds on a P4/2.6. The run time T(n) where n is the upper limit is
interesting - one would expect T(n) = n**1.5, and this seems to happen ;
here are a few values:

T(2**19) = 2234 ms
T(2**20) = 6157 ms
T(2**21) = 16937 ms
T(2**22) = 46766 ms
T(2**23) = 130218 ms

This is approximately

T(n) = 8.65 * 1E-6 * ( n ** 1.47 )




Vic

-> It was surprisingly easy to convert this program to XBasic - only 6 changes
-> were made
..................
-> For Mx=10,000 this runs too fast to measure even in interpreted mode.
-> Changing the upper limit to 100,000 results in a run time of 220
-> milliseconds on a P4/2.6. The run time T(n) where n is the upper limit is
-> interesting - one would expect T(n) = n**1.5, and this seems to happen ;
-> here are a few values:

-> T(2**19) = 2234 ms
-> T(2**20) = 6157 ms
-> T(2**21) = 16937 ms
-> T(2**22) = 46766 ms
-> T(2**23) = 130218 ms

-> This is approximately

-> T(n) = 8.65 * 1E-6 * ( n ** 1.47 )




-> Vic

Yes. I figured that (in QBasic notation) the run-time should be
approximately proportional to Mx^(3/2), where Mx is the greatest value
of the smallest number in the triple.

It appears that, if speed is measured by the rate at which the program
generates triples, then it is remarkably close to constant. As A (the
smallest of the 3 numbers) increases, greater numbers of triples are
generated for each value of A, on average. So, even though the rate of
change of A decreases, the rate at which triples are generated changes
very little.

dow

Some people suggested that I should make a tidier version of that very
fast Pythagorean Triples program (in QBasic), divided into modules, and
with some comments. So here it is. Tidiness comes at a price. This is
slightly, but only slightly, slower than the unstructured program. But
it is neater, and maybe easier to understand.

dow

---------------------------------------------------------

' Pythagorean Triples

' David O. Williams. 2004

' Calculates and prints integer triples, A, B, C, such that
' A < B < C, they have no common factors (>1), and A^2 + B^2 = C^2.
' List is printed in order of increasing A, then of B.
' A triple counter is also shown.

DECLARE FUNCTION NoComFacs& (X&, Y&)
DECLARE SUB PrintOut ()
DECLARE SUB Even ()
DECLARE SUB Odd ()

DEFLNG A-Z

Mn = 3' minimum value of smallest number in triple
Mx = 10000' maximum value of smallest number in triple

DIM SHARED A, B, C ' three numbers in triple (increasing order)
DIM SHARED N, Z# ' triple counter, and much-used constant

N = 0
Z# = SQR(2) - 1

CLS

FOR A = Mn TO Mx
SELECT CASE A MOD 4
CASE 0: Even
CASE 1, 3: Odd
' case 2 produces no valid triples
END SELECT
NEXT

END

SUB Even ' handles cases when A is even (and a multiple of 4)
S = A \ 2
FOR F = INT(SQR(Z# * S)) TO 1 STEP -1
G = S \ F
IF G * F = S THEN
IF (F XOR G) AND 1 THEN
IF NoComFacs(F, G) THEN
E = F + G
D = G - F
B = D * E
C = (E * E + D * D) \ 2
PrintOut
END IF
END IF
END IF
NEXT
END SUB

FUNCTION NoComFacs (X, Y) ' non-zero if X and Y have no common factors
U = Y
V = X
DO WHILE V > 1
W = U MOD V
U = V
V = W
LOOP
NoComFacs = V
END FUNCTION

SUB Odd ' handles cases when A is odd
FOR D = INT(SQR(Z# * A)) - 1 OR 1 TO 1 STEP -2
E = A \ D
IF E * D = A THEN
IF NoComFacs(D, E) THEN
B = ((E + D) * (E - D)) \ 2
C = (E * E + D * D) \ 2
PrintOut
END IF
END IF
NEXT
END SUB

SUB PrintOut
N = N + 1
PRINT N, , A, B, C
END SUB

----------------------------------------------

added a bunch of documentation to that Pythagorean Triple program of
mine, and streamlined it a bit further. If anyone out there is
interested in how this kind of thing works, I hope this will explain
it.

dow

-----------------------------------------------------

' PYTRIPLE.BAS - Pythagorean Triples

' David O. Williams. 2004
' XXXX@XXXXX.COM

' Calculates and prints integer triples, A, B, C, such that
' A < B < C, they have no common factors (>1), and A^2 + B^2 = C^2.
' The list is printed in order of increasing A, then of B.
' A counter of triples is also shown.

' A more detailed explanation is at the end of the main module.

DECLARE FUNCTION NoComFacs& (X&, Y&)
DECLARE SUB PrintOut ()
DECLARE SUB Even ()
DECLARE SUB Odd ()

DEFLNG A-Z

DIM SHARED A, B, C ' three numbers in triple (increasing order)
DIM SHARED N, Z#, Q, R ' counter, utility factor, and loop limits

CLS

DO
PRINT "Both values must be >= 0, and Maximum >= Minimum"
PRINT
INPUT "Minimum value of smallest number in triple"; Mn
INPUT "Maximum value of smallest number in triple"; Mx
PRINT
IF Mn >= 0 AND Mx >= 0 AND Mx >= Mn THEN EXIT DO
BEEP
PRINT "Illegal value/s!"
PRINT
LOOP

IF Mx > 2 AND (Mx > Mn OR (Mn <> 4 AND Mn MOD 4 <> 2)) THEN
PRINT "Count", , " A", " B", " C"
PRINT
ELSE
PRINT "There are no triples in this range!"
END
END IF

Z# = SQR(2) - 1
Q = INT(SQR(Z# * Mn)) - 1 OR 1
R = INT(SQR(Z# * Mn / 2))
N = 0

FOR A = Mn TO Mx
SELECT CASE A MOD 4
CASE 0: Even
CASE 1, 3: Odd
END SELECT
NEXT

END

'----------------------------------------------------------

' Brief explanation:

' Pythagorean triples, if they are in lowest terms with no common
' factors, can be written as: Odd#^2 + Even#^2 = Big#^2. Big# is the
' largest integer (corresponding to the hypotenuse of the right-angled
' triangle). However, Odd# may be smaller or larger than Even#.

' If the three above numbers have no common factors, it is possible
' to define two odd positive integers, D and E, with E > D, such that:

' Odd# = D * E
' Even# = (E^2 - D^2) / 2
' Big# = (E^2 + D^2) / 2

' These definitions satisfy Pythagoras's Theorem. It is possible to
' prove that all valid triples, in lowest terms, can be written this
' way, with odd-integer values of D and E. (They must both be odd, so
' that their product is Odd#.)

' If a triple is written as A, B, C, with A < B < C, C must correspond
' to Big#, but A can be either Odd# or Even#, and B will be the other.

' The program treats these two possibilities separately. If A is odd,
' so it must be Odd#, the program simply searches for two odd integers
' whose product is A. After confirming that they have no common
' factors, which would mean that the triple is not in lowest terms,
' the program calculates B and C from them, and prints them out.

' If A is even, it is useful to define two further numbers, F and G,
' with G > F, such that:

' D = G - F
' E = F + G

' This means that:

' F = (E - D) / 2
' G = (D + E) / 2

' Since D and E are both odd, their sum and difference are both even,
' so F and G are integers. However, the sum and difference of F and G
' are E and D, which are odd, which