Pythagorean triples: Difference between revisions
Content added Content deleted
(Add PHP) |
(Added Pythagorean triples for EDSAC) |
||
Line 1,017: | Line 1,017: | ||
Up to 10000000000: 14557915466 triples, 702304875 primitives. |
Up to 10000000000: 14557915466 triples, 702304875 primitives. |
||
Up to 100000000000: 161750315680 triples, 7023049293 primitives.</pre> |
Up to 100000000000: 161750315680 triples, 7023049293 primitives.</pre> |
||
=={{header|EDSAC order code}}== |
|||
Not much optimization is done in this code, which is best run on a fast simulator |
|||
if the maximum perimeter is more than 10^5 or so. |
|||
A maximum perimeter of 10^7 takes 340 million EDSAC orders |
|||
(about 6 days on the original EDSAC). |
|||
The number of primitive triples divided by the maximum perimeter seems to tend to a limit, |
|||
which looks very much like (ln 2)/pi^2 = 0.07023049277 (see especially the FreeBASIC output below). |
|||
<lang edsac> |
|||
[Pythagorean triples for Rosetta code. |
|||
Counts (1) all Pythagorean triples (2) primitive Pythagorean triples, |
|||
with perimeter not greater than a given value. |
|||
Library subroutine M3, Prints header and is then overwritten. |
|||
Here, the last character sets the teleprinter to figures.] |
|||
..PZ [simulate blank tape] |
|||
PFGKIFAFRDLFUFOFE@A6FG@E8FEZPF |
|||
@&*!MAX!PERIM!!!!!TOTAL!!!!!!PRIM@&#. |
|||
..PZ |
|||
[Library subroutine P7, prints long strictly positive integer; |
|||
10 characters, right justified, padded left with spaces. |
|||
Closed, even; 35 storage locations; working position 4D.] |
|||
T 56 K |
|||
GKA3FT26@H28#@NDYFLDT4DS27@TFH8@S8@T1FV4DAFG31@SFLDUFOFFFSFL4F |
|||
T4DA1FA27@G11@XFT28#ZPFT27ZP1024FP610D@524D!FO30@SFL8FE22@ |
|||
[Subroutine for positive integer division. |
|||
Input: 4D = dividend, 6D = divisor. |
|||
Output: 4D = remainder, 6D = quotient. |
|||
37 locations; working locations 0D, 8D.] |
|||
T 100 K |
|||
GKA3FT35@A6DU8DTDA4DRDSDG13@T36@ADLDE4@T36@T6DA4DSDG23@ |
|||
T4DA6DYFYFT6DT36@A8DSDE35@T36@ADRDTDA6DLDT6DE15@EFPF |
|||
[Subroutine to return GCD of two non-negative 35-bit integers. |
|||
Input: Integers at 4D, 6D. |
|||
Output: GCD at 4D; changes 6D. |
|||
41 locations; working location 0D.] |
|||
T 200 K |
|||
GKA3FT39@S4DE37@T40@A4DTDA6DRDSDG15@T40@ADLDE6@T40@A6DSDG20@T6D |
|||
T40@A4DSDE29@T40@ADRDTDE16@S6DE39@TDA4DT6DSDT4DE5@A6DT4DEFPF |
|||
[************************ ROSETTA CODE TASK ************************* |
|||
Subroutine to count Pythagorean triples with given maximum perimeter. |
|||
Input: 0D = maximum perimeter. |
|||
Output: 4D = number of triples, 6D = number of primitive. |
|||
0D is changed. |
|||
Must be loaded at an even address. |
|||
Uses the well-known fact that a primitive Pythagorean triple is of the form |
|||
(m^2 - n^2, 2*m*n, m^2 + n^2) where m, n are coprime and of opposite parity.] |
|||
T 300 K |
|||
G K |
|||
A 3 F [make link] |
|||
E 16 @ [jump over variables and constants] |
|||
[Double values are put here to ensure even address] |
|||
[Variables] |
|||
[2] P F P F [maximum perimeter] |
|||
[4] P F P F [total number of Pythagorean triples] |
|||
[6] P F P F [number of primitive Pythagorean triples] |
|||
[8] P F P F [m] |
|||
[10] P F P F [n] |
|||
[Constants] |
|||
[12] P D P F [double-value 1] |
|||
[14] P1F P F [double-value 2] |
|||
[Continue with code] |
|||
[16] T 69 @ [plant link for return] |
|||
A D [load maximum perimeter] |
|||
T 2#@ [store locally] |
|||
T 4#@ [initialize counts of triangles to 0] |
|||
T 6#@ |
|||
A 12#@ [load 1] |
|||
T 8#@ [m := 1] |
|||
[Next m, inc by 1] |
|||
[23] T F [clear acc] |
|||
A 8#@ [load m] |
|||
A 12#@ [add 1] |
|||
T 8#@ [update m] |
|||
H 8#@ [mult reg := m] |
|||
C 12#@ [acc := m AND 1] |
|||
A 12#@ [add 1] |
|||
T 10#@ [n := 1 if m even, 2 if m odd] |
|||
[Here to count triangles arising from m, n. |
|||
It's assumed m and n are known coprime.] |
|||
[31] A 31 @ [call the count subroutine,] |
|||
G 70 @ [result is in 6D] |
|||
S 6 D [load negative count] |
|||
G 40 @ [jump if count > 0] |
|||
[No triangles found for this n. |
|||
If n = 1 or 2 then whole thing is finished. |
|||
Else move on to next m.] |
|||
T F [clear acc] |
|||
A 14#@ [load 2] |
|||
S 10#@ [2 - n] |
|||
G 23 @ [if n > 2, go to next m] |
|||
E 64 @ [if n <= 2, exit] |
|||
[Found triangles, count is in 6D] |
|||
[40] T F [clear acc] |
|||
A 4#@ [load total count] |
|||
A 6 D [add count just found] |
|||
T 4#@ [update total count] |
|||
A 6#@ [load primitive count] |
|||
A 12#@ [add 1] |
|||
T 6#@ [update primitive count] |
|||
[47] T F [clear acc] |
|||
A 10#@ [load n] |
|||
A 14#@ [add 2] |
|||
U 10#@ [update n] |
|||
S 8#@ [is n > m?] |
|||
E 23 @ [if so, loop back for next m] |
|||
[Test whether m and n are coprime.] |
|||
T F [clear acc] |
|||
A 8#@ [load m] |
|||
T 4 D [to 4D for GCD routine] |
|||
A 10#@ [load n] |
|||
T 6 D [to 6D for GCD routine] |
|||
A 58 @ [call GCD routine,] |
|||
G 200 F [GCD is returned in 4D] |
|||
A 4 D [load GCD] |
|||
S 14#@ [is GCD = 1? (test by subtracting 2)] |
|||
E 47 @ [no, go straight to next n] |
|||
G 31 @ [yes, count triangles, then next n] |
|||
[64] T F [exit, clear acc] |
|||
A 4#@ [load total number of triples] |
|||
T 4 D [return in 4D] |
|||
A 6#@ [load number of primitive triples] |
|||
T 6 D [return in 6D] |
|||
[69] E F |
|||
[2nd-level subroutine to count triangles arising from m, n. |
|||
Assumes m, n are coprime and of opposite parity, |
|||
and m is in the multiplier register. |
|||
Result is returned in 6D.] |
|||
[70] A 3 F [make and plant link for return] |
|||
T 91 @ |
|||
A 2#@ [acc := maximum perimeter] |
|||
T 4 D [to 4D for division routine] |
|||
A 8#@ [load m] |
|||
A 10#@ [add n] |
|||
T D [m + n to 0D] |
|||
V D [acc := m*(m + n)] |
|||
[Need to shift product 34 left to restore integer scaling. |
|||
Since we want 2*m*(m+n), shift 35 left.] |
|||
L F [13 left (maximum possible)] |
|||
L F [13 more] |
|||
L 128 F [9 more] |
|||
T 6 D [perimeter to 6D for division routine] |
|||
A 4 D [load maximum perimeter] |
|||
S 6 D [is perimeter > maximum?] |
|||
G 89 @ [quick exit if so] |
|||
T F [clear acc] |
|||
A 86 @ [call division routine,] |
|||
G 100 F [leaves count in 6D] |
|||
E 91 @ [jump to exit] |
|||
[89] T F [acc := 0] |
|||
T 6 D [return count = 0] |
|||
[91] E F |
|||
[Main routine. Load at an even address.] |
|||
T 500 K |
|||
G K |
|||
[The initial maximum perimeter is repeatedly multiplied by 10] |
|||
[0] P50F PF [initial maximum perimeter <---------- EDIT HERE] |
|||
[2] P 3 F [number of values to calculate <---------- EDIT HERE] |
|||
[3] P D [1] |
|||
[4] P F P F [maximum perimeter] |
|||
[6] P F P F [total number of triples] |
|||
[8] P F P F [number of primitive triples] |
|||
[10] P F [negative count of values] |
|||
[11] # F [figures shift] |
|||
[12] @ F [carriage return] |
|||
[13] & F [line feed] |
|||
[14] K 4096 F [null char] |
|||
[Enter with acc = 0] |
|||
[15] S 2 @ [initialize a negative counter] |
|||
T 10 @ [(standard EDSAC practice)] |
|||
A #@ [initialize maximum perimeter] |
|||
T 4#@ |
|||
[19] T F [clear acc] |
|||
A 4#@ [load maximum perimeter] |
|||
T D [to 0D for subroutine] |
|||
A 22 @ [call subroutine to count triples] |
|||
G 300 F |
|||
A 4 D [returns total number in 4D] |
|||
T 6#@ [save locally] |
|||
A 6 D [returns number of primitive in 6D] |
|||
T 8#@ [save locally] |
|||
[Print the result] |
|||
A 4#@ [load maximum perimeter] |
|||
T D [to 0D for print subroutine] |
|||
A 30 @ [call print subroutine] |
|||
G 56 F |
|||
A 6#@ [repeat for total number of triples] |
|||
T D |
|||
A 34 @ |
|||
G 56 F |
|||
A 8#@ [repeat for number of primitive triples] |
|||
T D |
|||
A 38 @ |
|||
G 56 F |
|||
O 12 @ |
|||
O 13 @ |
|||
A 10 @ [load negative count] |
|||
A 3 @ [add 1] |
|||
E 53 @ [out if reached 0] |
|||
T 10 @ [else update count] |
|||
A 4#@ [load max perimeter] |
|||
U D [temp store] |
|||
L 1 F [times 4] |
|||
A D [times 5] |
|||
L D [times 10] |
|||
T 4#@ [update] |
|||
E 19 @ [loop back] |
|||
[53] O 14 @ [done; print null to flush printer buffer] |
|||
Z F [stop] |
|||
E 15 Z [define entry point] |
|||
P F [acc = 0 on entry] |
|||
</lang> |
|||
{{out}} |
|||
<pre> |
|||
MAX PERIM TOTAL PRIM |
|||
100 17 7 |
|||
1000 325 70 |
|||
10000 4858 703 |
|||
100000 64741 7026 |
|||
1000000 808950 70229 |
|||
10000000 9706567 702309 |
|||
</pre> |
|||
=={{header|Eiffel}}== |
=={{header|Eiffel}}== |