fast pythagorean triple generator

Here's one for people who like interesting algorithms. It generates
triples of integers, A, B, and C, such that A^2 + B^2 = C^2. However,
it does it without using SQR, which is slow. The full algorithm is
explained in the program.

I remember thinking very hard about this when I wrote it about 3 years
ago. Now, I'd have to think hard again to understand it...

Enjoy!

dow

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

' PYTRIPLE.BAS - Pythagorean Triples

' David O. Williams. 2004

' Version date: 2004 Sep 05

' 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& (P&, Q&)
DECLARE SUB PrintOut (A&, B#, C#)
DECLARE SUB DoOdd (A&)
DECLARE SUB DoEven (A&)

DEFLNG A-Z

CLS

DO
PRINT "Legal ranges: 0 <= Minimum <= Maximum <= 2^27 (134217728)"
PRINT
INPUT "Minimum value of smallest number in triple"; Mn
INPUT "Maximum value of smallest number (or ENTER for same)"; Mx
PRINT
IF Mx = 0 THEN Mx = Mn
IF Mn >= 0 AND Mx >= Mn AND Mx <= 2 ^ 27 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"; TAB(16); "A"; TAB(33); "B"; TAB(56); "C"
PRINT
ELSE
PRINT "There are no triples in this range!"
END
END IF

CONST Z# = 2.414213562373095# ' SQR(2) + 1

' main loop
FOR A = Mn TO Mx
SELECT CASE A AND 3
CASE 0: DoEven A
CASE 1, 3: DoOdd A
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. In
' calculating C, the equation C = B + D^2 is used, which is easy to
' prove since B is Even#.

' 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 implies that the parities of F
' and G must be opposites, so one is odd and the other even.

' Since, in this case, A is the same as Even#, it is given by:

' A = (E^2 - D^2) / 2

' Writing G - F for D and G + F for E, and simplifying, this gives:

' A / 2 = F * G

' Since F and G are of opposite parity, this means that A / 2 must be
' even, so A must be a multiple of 4. There are no valid triples when
' (A AND 3) = 2.

' When A is a multiple of 4, the program looks for factors of A / 2.
' It checks that they are of opposite parity, and have no common
' factors. It then uses these values of F and G to calculate B and
' C, using the easily-proved formulae:

' B = G^2 - F^2
' C = G^2 + F^2

' (Strictly, the parity check could be omitted. Since they are factors
' of an even number, at least one of the found values of F and G must
' be even. The possibility that they are both even would be detected
' by the common-factor test. However, the parity check is much simpler
' and faster than the common-factor routine, so having it in the
' program saves some time.)

' The FOR ... NEXT loops in the DoOdd and DoEven SUBroutines search
' for the factors D and E, or F and G, respectively, The loops run
' from high to low values of the counting variables, since this makes
' the triples appear in the desired order. Also, the ranges of the
' loops are limited so that only triples in which B > A appear. The
' variable J governs these ranges. It is initialized when each
' SUB is first called, and is incremented in the SUBroutine as
' the value of A increases. The number L is a pre-calculated
' limit. When A (or S, which is just A / 2 when A is even) passes
' this limit, J is incremented appropriately, and a new value
' of L is calculated. On most of the occasions when the
' subroutines are used, it is not necessary to increment J, so
' all that has to be done is a simple comparison, e.g IF A > L THEN,
' which turns out not to be true. This saves a lot of time. Also,
' this method does not use the slow operation SQR repeatedly. A
' full mathematical treatment of the situation uses the number
' (SQR(2) + 1) several times. This number is therefore treated as a
' CONSTant in the program, named Z#.

' This version uses double-precision floating-point variables for B
' and C, and in the calculations that lead to them. This allows
' triples with very large numbers to be calculated. The upper limit
' for A is set at 2^27, which is the limit beyond which QBasic cannot
' do the calculations completely accurately.

' A simpler version of the program, capable of going up only to A =
' 2^16 (65536), can be made by making B, C, and X into long-integer
' variables (dropping the "#" suffixes), and removing the CDBL
' commands. This requires small changes in the SUBs DoOdd, DoEven,
' and PrintOut, and, in the main module, to the DECLARE SUB PrintOut
' statement and the statements of the maximum limit of A. This simpler
' version shows the algorithm more clearly.

' = end =


SUB DoEven (A)

STATIC H, J, K, L

S = A \ 2

IF NOT H THEN ' initialize loop limit
J = INT(SQR(S / Z#)) ' next integer below SQR(S/Z#)
K = J + 1
L = INT(K * K * Z#) ' when S passes this value, increment J
H = -1 ' flag shows initialization done
ELSEIF S > L THEN ' increment loop limit
J = K
K = J + 1
L = INT(K * K * Z#)
END IF

FOR F = J TO 1 STEP -1 ' search for factors of S (which is A/2)
IF S MOD F = 0 THEN
G = S \ F
IF (F XOR G) AND 1 THEN ' parity check
IF NoComFacs(F, G) THEN ' common-factor check
X# = CDBL(G) * G
Y = F * F
B# = X# - Y
C# = X# + Y
PrintOut A, B#, C#
END IF
END IF
END IF
NEXT

END SUB

SUB DoOdd (A)

STATIC H, J, K, L

IF NOT H THEN ' initialize loop limit
J = INT(SQR(A / Z#)) - 1 OR 1 ' next *odd* integer < SQR(A/Z#)
K = J + 2 ' loop steps in 2's: odd numbers only
L = INT(K * K * Z#) ' when A passes L, J must be incremented
H = -1 ' flag shows initialization done
ELSEIF A > L THEN ' increment loop limit
J = K
K = J + 2
L = INT(K * K * Z#)
END IF

FOR D = J TO 1 STEP -2 ' search for factors of A (must be odd)
IF A MOD D = 0 THEN
E = A \ D
IF NoComFacs(D, E) THEN ' common-factor check
B# = ((E - D) \ 2) * CDBL(E + D)' divide soon: overflow late
C# = D * D + B#
PrintOut A, B#, C#
END IF
END IF
NEXT

END SUB

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

SUB PrintOut (A, B#, C#)
STATIC N
N = N + 1
PRINT N; TAB(15); A; TAB(32); B#; TAB(55); C#
END SUB

-------------------------------------------------------
.

Relevant Pages

  • Pythagorean triples
    ... Here's a program that finds and prints out "Pythagorean Triples", ... DECLARE SUB DoOdd ... ' main loop ... ' to define two odd positive integers, D and E, with E> D, such that: ...
    (comp.lang.basic.misc)
  • Re: printing double sided reports
    ... Here is some code that allows you to have different margins on even and odd ... Dim MoveMargin As Integer ... Private Sub ReportHeader_Format(Cancel As Integer, ...
    (microsoft.public.access.reports)
  • Re: Odd/even page different detail section
    ... works overlay the odd set over the even set in the detail section. ... Dim MoveMargin As Integer ... Code the Report Header Format event: ... Private Sub ReportHeader_Format(Cancel As Integer, ...
    (microsoft.public.access.reports)
  • Re: conditional formatting w/ check for beginning odd or even
    ... ' first row is even ... Sub abcEven() ... Dim even As Long, odd As Long ... One macro gives the ...
    (microsoft.public.excel.programming)
  • Re: Integer right triangles
    ... Since the difference between any two adjacent squares is an odd ... If K is a divisor of A^2 and K and A are of same parity etc... ... This is a valuable method to get triples (by scanning all divisors ... All primitive triples =1) are generated once each by: ...
    (rec.puzzles)