228 lines
9.1 KiB
Plaintext
228 lines
9.1 KiB
Plaintext
[Farey sequence for Rosetta Code website.
|
|
Get number of terms by using Euler's totient function.
|
|
EDSAC program, Initial Orders 2.]
|
|
|
|
[Euler's totient function for each n = 1..1000 is calculated here as follows.
|
|
A wheel is defined for each prime p < sqrt(1000), i.e. for p <= 31.
|
|
When n = 0, the wheels are all 0. When n is incremented:
|
|
(i) the totient is initialized to n
|
|
(ii) the wheel for each prime p is incremented modulo p.
|
|
A prime p therefore divides n iff the wheel for p is 0. In this case:
|
|
(1) the totient is multiplied by (1 - 1/p)
|
|
(2) n is reduced by dividing it by p as many times as possible.
|
|
When all primes p have been tested, the reduced n must be either:
|
|
(a) 1, in which case the totient is finished; or
|
|
(b) a prime q > 31, in which case the totient is multiplied by (1 - 1/q).]
|
|
|
|
[Library subroutine M3, prints header, terminated by blank row of tape.]
|
|
PFGKIFAFRDLFUFOFE@A6FG@E8FEZPF
|
|
*!!!!!ORDER!!!!!TERMS@&#..
|
|
[PZ]
|
|
T 56 K
|
|
[Library subroutine P7, prints double-word integer > 0.
|
|
10 characters, right justified, padded left with spaces.
|
|
Closed, even; 35 storage locations; working position 4D.]
|
|
GKA3FT26@H28#@NDYFLDT4DS27@TFH8@S8@T1FV4DAFG31@SFLDUFOFFFSF
|
|
L4FT4DA1FA27@G11@XFT28#ZPFT27ZP1024FP610D@524D!FO30@SFL8FE22@
|
|
|
|
[Subroutine (not from library) for integer short division.
|
|
Input: dividend at 4F, divisor at 6F
|
|
Output: remainder at 4F, quotient at 6F
|
|
Working location 0D. 37 locations.]
|
|
T 100 K
|
|
GKA3FT34@A6FUFT35@A4FRDS35@G13@T1FA35@LDE4@T1FT6FA4FS35@G22@
|
|
T4FA6FA36@T6FT1FAFS35@E34@T1FA35@RDT35@A6FLDT6FE15@EFPFPD
|
|
|
|
[Put address of primes at 53.
|
|
Primes are therefore referred to by code letter B.]
|
|
T 53 K
|
|
P 160 F
|
|
T 160 K
|
|
P 11 F [number of primes (as address)]
|
|
P1FP1DP2DP3DP5DP6DP8DP9DP11DP14DP15D
|
|
|
|
[Put address of wheels at 54.
|
|
Wheels are therefore referred to by code letter C.
|
|
Number of wheels = number of primes, at the moment 11]
|
|
T 54 K
|
|
P 180 F
|
|
|
|
[Main routine]
|
|
T 200 K
|
|
G K
|
|
[Long variable]
|
|
[0] P F P F [sum of Euler's totient function over all n]
|
|
[Short variables]
|
|
[2] P F [n = order of Farey sequence]
|
|
[3] P F [reduced n, as prime factors are taken out]
|
|
[4] P F [partial totient; initially n, finally Euler's phi(n)]
|
|
[5] P F [current prime p]
|
|
[6] P F [residue of n by prime p]
|
|
[7] P F [negative counter for steps]
|
|
[8] P F [negative counter within step]
|
|
[9] P F [negative counter for primes]
|
|
[Short constants]
|
|
[10] P D [1]
|
|
[11] P 100 F [step, as an address (for convenience)]
|
|
[12] P 10 F [number of steps, as an address]
|
|
[13] # F [figure shift]
|
|
[14] @ F [carriage return]
|
|
[15] & F [line feed]
|
|
[16] K 4096 F [teleprinter null]
|
|
[17] A C [order to read first wheel]
|
|
[18] T C [order to write first wheel]
|
|
[19] A 1 B [order to read first prime]
|
|
|
|
[Subroutine to multiply partial totient by (1 - 1/p)]
|
|
[20] A 3 F
|
|
T 31 @
|
|
A 4 @ [load partial totient]
|
|
T 4 F [to dividend]
|
|
A 5 @ [load prime p]
|
|
T 6 F [to divisor]
|
|
[26] A 26 @ [for return from next]
|
|
G 100 F [call division routine]
|
|
A 4 @ [partial totient again]
|
|
S 6 F [subtract quotient]
|
|
T 4 @ [update partial totient]
|
|
[31] E F [exit with acc = 0]
|
|
|
|
[Enter with accumulator = 0]
|
|
[Reset all wheels to 0, working from 31 down to 2.]
|
|
[32] A B [load number of wheels]
|
|
S 2 F [dec by 1]
|
|
[34] A 18 @ [make order 'T m C' for address m]
|
|
T 36 @ [plant order]
|
|
[36] T C [reset this wheel]
|
|
A 36 @ [get order again]
|
|
S 2 F [dec address by 1]
|
|
S 18 @ [compare with order for first wheel]
|
|
E 34 @ [loop back till done]
|
|
[Initialize sum to 1]
|
|
T F [clear acc]
|
|
T #@ [clear sum (both words + sandwich bit)]
|
|
A 10 @ [load 1 (single word)]
|
|
T @ [to sum (low word)]
|
|
T 2 @ [order of Farey sequence := 0]
|
|
S 12 @ [load negative number of steps (typically -10)]
|
|
[Here acc = negative step count]
|
|
[47] T 7 @ [update negative step count]
|
|
S 11 @ [load negative step size (typically -100)]
|
|
[Here acc = negative count within a step]
|
|
[49] T 8 @
|
|
A 2 @ [inc n]
|
|
A 10 @
|
|
U 3 @ [initialize reduced n := n]
|
|
U 4 @ [initialize partial totient := n]
|
|
T 2 @ [update n]
|
|
|
|
[Loop through primes p. Inc wheel for prime p by 1 mod p.
|
|
If wheel = 0, then p divides n.
|
|
If so, update partial totient and reduced n.]
|
|
S B
|
|
T 9 @ [initialize count]
|
|
A 19 @ [order to read first prime]
|
|
T 64 @ [plant in code]
|
|
A 17 @ [order to read first wheel]
|
|
T 66 @ [plant in code]
|
|
A 18 @ [order to write first wheel]
|
|
T 88 @ [plant in code]
|
|
[63] T F
|
|
[64] A F [load prime]
|
|
T 5 @ [store]
|
|
[66] A C [read wheel (residue of n mod p)]
|
|
A 10 @ [inc]
|
|
U 6 @ [store locally]
|
|
S 5 @ [reached p yet?]
|
|
G 86 @ [skip if not]
|
|
|
|
[Here if p divides n.
|
|
Need to multiply partial totient by (1 - 1/p)
|
|
and divide reduced n by highest possible power of p.]
|
|
T F [acc := 0]
|
|
T 6 @ [wrap residue from p to 0]
|
|
[Update partial totient, multiply by (1 - 1/p)]
|
|
[73] A 73 @ [for return from next]
|
|
G 20 @ [call subroutine]
|
|
[Divide reduced n by p as many times as possible
|
|
(it must be divisible by p at least once)]
|
|
[75] A 3 @ [load reduced n]
|
|
T 4 F [to dividend]
|
|
A 5 @ [load prime p]
|
|
T 6 F [to divisor]
|
|
[79] A 79 @ [for return]
|
|
G 100 F [call division routine; clears acc]
|
|
S 4 F [load negative of remainder]
|
|
G 86 @ [stop dividing if remainder > 0]
|
|
A 6 F [quotient from division]
|
|
T 3 @ [update reduced n]
|
|
E 75 @ [try another division]
|
|
[86] T F [clear acc]
|
|
A 6 @ [get residue for this prime]
|
|
[88] T C [write back]
|
|
A 9 @ [load negative prime count]
|
|
A 2 F [inc count]
|
|
E 103 @ [out if done all primes]
|
|
T 9 @ [else update count]
|
|
A 64 @ [inc addresses in the above code]
|
|
A 2 F
|
|
T 64 @
|
|
A 66 @
|
|
A 2 F
|
|
T 66 @
|
|
A 88 @
|
|
A 2 F
|
|
T 88 @
|
|
E 63 @ [loop for next prime]
|
|
|
|
[Tested all primes up to 31 for this n.
|
|
Reduced n is now either 1 or a prime > 31]
|
|
[103] T F
|
|
A 3 @ [get reduced]
|
|
S 2 F [subtract 2]
|
|
G 111 @ [skip if reduced n = 1]
|
|
A 2 F [else restore value]
|
|
T 5 @ [copy to prime p]
|
|
[109] A 109 @ [for return from next]
|
|
G 20 @ [call routine to update partial totient]
|
|
[Update sum of Euler's totient over 1..n.
|
|
Note sum is double word, while totient is single word.
|
|
Totient is converted to double before adding to sum.]
|
|
[111] T F [clear acc]
|
|
T D [clear 0D (i.e. 0F, 1F and sandwich bit)]
|
|
A 4 @ [load totient (single word)]
|
|
T F [to 0F]
|
|
A D [load totient from 0D as double word]
|
|
A #@ [add to sum]
|
|
T #@ [update sum]
|
|
|
|
[On to next n]
|
|
A 8 @ [load negative count]
|
|
A 2 F [add 1]
|
|
G 49 @ [loop until count = 0]
|
|
|
|
[Here when finished this step.
|
|
Typically, n has increased by 100.
|
|
Show n and the sum of Euler's totient.
|
|
Note accumulator = 0 here.]
|
|
T D [clear 0D (i.e. 0F, 1F and sandwich bit)]
|
|
A 2 @ [load n (single word)]
|
|
T F [to 0F; now 0D = n for printing]
|
|
[124] A 124 @ [for return from next]
|
|
G 56 F [call library subroutine to print n]
|
|
A #@ [load sum (double word)]
|
|
T D [to 0D for printing]
|
|
[128] A 128 @ [for return from next]
|
|
G 56 F [call library subroutine to print sum]
|
|
O 14 @ [print new line (CR, LF)]
|
|
O 15 @
|
|
[On to next step]
|
|
A 7 @ [load negative step count]
|
|
A 2 F [add 1]
|
|
G 47 @ [loop until count = 0]
|
|
[Here when finished whole thing]
|
|
[135] O 16 @ [output null to flush teleprinter buffer]
|
|
Z F [stop]
|
|
E 32 Z [define start of execution]
|
|
P F [start with accumulator = 0]
|