Text 161, 255 rader
Skriven 2004-08-31 21:54:48 av DAVID WILLIAMS (1:250/514)
Kommentar till en text av MILES MAXTED
Ärende: Re: Pythagorean Triples
===============================
-> DW> I keep trying to streamline this thing further, and I do sometimes
-> DW> manage to shave just a little bit off the execution time. Here's the
-> DW> latest, and fastest, version. Can anyone do better?
-> Mmmm - this one runs here between 0 - 10000 in 4 seconds -
-> probably because of the shift from integers. It's getting to be
-> an impressive case-study, nevertheless, and is safely filed away
-> against some future use....
I have the impression that sometimes changes that make the program
faster on my machine make it slower on yours, and vice versa. Could
be...
Anyway, I've been trying to "modularize" the thing, to satisfy people
who like programs that way. This version is *nearly* as fast, on my
machine, as the previous fastest one. Who knows? Maybe on your machine,
it will be even faster.
Incidentally, I haven't shifted away from integers, except in that one
version where I tried to make it handle *huge* numbers. That version
was *really* slow!
dow
--------------------------------------------------
' PYTRIPLE.BAS - Pythagorean Triples
' David O. Williams. 2004
' david.williams@ablelink.org
' 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 Odd (A&)
DECLARE SUB Even (A&)
DEFLNG A-Z
CLS
DO
PRINT "Legal ranges: 0 <= Minimum <= Maximum <= 2^16 (65536)"
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 ^ 16 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: CALL Even(A)
CASE 1, 3: CALL Odd(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.
' 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 MOD 4 = 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 Odd and Even 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 long-integer format for all of the variables
' (except Z#), which sets a limit on the largest numbers that can be
' handled. If A exceeds 2^16 (65536), overflows occur when the
' values of B and C are calculated, since the largest values of B and
' C can pass 2^31, which is more than the maximum value for a long
' integer. For this reason, an upper limit of 65536 is placed on the
' value of A.
' = end =
SUB Even (A)
STATIC H, J, K, L
S = A \ 2
IF NOT H THEN
H = -1
J = INT(SQR(S / Z#))
K = J + 1
L = INT(K * K * Z#)
ELSEIF S > L THEN
J = K
K = J + 1
L = INT(K * K * Z#)
END IF
FOR F = J TO 1 STEP -1
IF S MOD F = 0 THEN
G = S \ F
IF (F XOR G) AND 1 THEN
IF NoComFacs(F, G) THEN
X = G * G
Y = F * F
B = X - Y
C = X + Y
PrintOut A, B, C
END IF
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 Odd (A)
STATIC H, J, K, L
IF NOT H THEN
H = -1
J = INT(SQR(A / Z#)) - 1 OR 1
K = J + 2
L = INT(K * K * Z#)
ELSEIF A > L THEN
J = K
K = J + 2
L = INT(K * K * Z#)
END IF
FOR D = J TO 1 STEP -2
IF A MOD D = 0 THEN
E = A \ D
IF NoComFacs(D, E) THEN
B = ((E - D) \ 2) * (E + D)
C = D * D + B
PrintOut A, B, C
END IF
END IF
NEXT
END SUB
SUB PrintOut (A, B, C)
STATIC N
N = N + 1
PRINT N; TAB(15); A; TAB(32); B; TAB(55); C
END SUB
------------------------------------------------------
--- Platinum Xpress/Win/WINServer v3.0pr5
* Origin: The Bayman BBS,Toronto, (416)698-6573 - 1:250/514 (1:250/514)
|