[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]