Pythagorean triples: Difference between revisions
Added Pythagorean triples for EDSAC
(Add PHP) |
(Added Pythagorean triples for EDSAC) |
||
Line 1,017:
Up to 10000000000: 14557915466 triples, 702304875 primitives.
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}}==
| |||