[Normal distribution for Rosetta Code EDSAC program, Initial Orders 2] [================================================================================== Uses an accept-reject method, which requires only logarithms (no trig functions). Let u, v be independent uniform variates in (0, 1). Let x = -ln(u), y = -ln(v). Accept x iff (x - 1)^2 <= 2y. If x is accepted, negate it with probability 1/2. Then x is normally distributed with mean 0 and standard deviation 1. The algorithm is modified for this EDSAC version: (1) Uses EDSAC library subroutine L1 to calculate (1/32)log_2() rather than ln() (2) Since real numbers on EDSAC are restricted to the interval [-1, 1), scales so that the standard deviation is 1/4, and reject values >= 4 s.d. from the mean. In the histogram, counts are divided by 16, with rounding. On the EDSAC PC simulator, takes 2.75 EDSAC hours to find 4096 normal variates. [=================================================================================] [Arrange the storage] T46K P70F [N parameter: library subroutine P1 to print +ve fraction] T47K P200F [M parameter: main routine and dependent subroutines] T49K P92F [L parameter: library subroutine L1 to calculate log_2] T51K P130F [G parameter: generator for pseudo-random numbers] T52K P168F [A parameter: library subroutine S2 for square root] T54K P190F [C parameter: constants read in by library subroutine R9] [Library subroutine R9, reads non-negative integers at load time. Fractions can be read by converting to integers (multiply by 2^34). 15 locations, must be loaded at location 56.] T56KGKT20FVDL8FA40DUDTFI40FA40FS39FG@S2FG23FA5@T5@E4@ [Library subroutine M3, prints header at load time. M3 and header are then overwritten.] PFGKIFAFRDLFUFOFE@A6FG@E8FEZPF *MEAN#!K0BL!!*SD#!K0M25BL@&#..PZ [========================== C parameter ==========================] [Tell R9 where to store integers read from tape] E69K T#C ['T m D' in WWG, but this also works] [List of integers; separated by F; end list with #TZ] 123456789F5F1549082005F11908177887F858993#TZ [0] [seed for random number generator] [2] [minimum argument for library subroutine L1] [4] [1/(16*ln(2)) for accept-reject] [6] [ln(2) for scaling standard deviation] [8] [0.00005 for rounding in print routine] [========================== M parameter ==========================] E25K TM GK [0] P4096F [number of data, in address field; code below assumes 4096] [1] PF [negative count of data] [2] PF [worlspace, low word] [3] PF [workspace, high word] [4] PF PF [auxiliary for accept-reject algorithm] [6] PF PF [sum of value/count, for mean] [8] PF PF [sum of value^2/count, for variance] [10] PFPFPFPFPFPFPFPFPFPFPFPFPFPFPFPF [16 bins for histogram] [26] CF [11110...0 binary, to isolate bin index; also prints colon] [27] A10@ [A order for bin{0}] [28] AF [A order for bin{16}] [29] P8F [8 in address field] [30] MF [subtract from A order to make T order; also prints dot] [31] WF [1/8] [32] FF [-15/16, mid value of histogram bin{0}] [33] PF [mid value of current bin] [Teleprinter] [34] #F [set figures mode] [35] PF [36] AF [minus (in figures mode)] [37] !F [space] [38] @F [carriage return] [39] &F [line feed] [Enter with acc = 0] [40] S@ T1@ [store negated data count in address field] A#C T4D [copy seed to 4D for PRNG] [44] A44@ GG [initialize PRNG] T6#@ T8#@ [clear sum and sum of squares] O34@ [set teleprinter to figures] [Start of loop to generate normal variates] [49] TF [clear acc] [Calculate and store logarithms of uniform variates] [50] A50@ G176@ SD T2#@ [corresponds to x = -ln(u)] [54] A54@ G176@ SD T4#@ [corresponds to y = -ln(v)] [Accept or reject] A4#C RD [acc := (1/32*ln(2))] S2#@ [subtract x] TD HD ND [acc := negated square] H4#C V4#@ [add 2*(auxiliary value)] G49@ [reject if result < 0] [First log is accepted; multiply by 8*ln(2) so that s.d. becomes 0.25] [Reject values outside (-1,1) (that is >= 4 s.d. from mean, unlikely)] T6F [clear acc] H6#C V2#@ [times ln(2)] S31@ E49@ [if product >= 1/8 (unlikely) reject and try again] A31@ [restore acc after test] L2F [shift 3 left to complete scaling] YF T2#@ [round and store back] [Finally change sign with probability 1/2] T6F T4D A2F T4F [pass range = 2 to PRNG] [80] A80@ G1G [call PRNG, 0 or 1 returned in 0D] SD E87@ [don't change sign if 0D = 0] T6F S2#@ T2#@ [change sign, so value < 0] [87] [Here with acc = 0.] [To print the values, replace the following jump with a no-op (X F)] E94@ A2#@ TD [pass abs(value) to print routine] [90] A90@ G191@ O38@ O39@ [print CR, LF] [94] A2#@ R1024F YF [load value, divide by count, round] A6#@ T6#@ [update sum] H2#@ V2#@ R1024F YF A8#@ T8#@ [update sum of squares] [Inc count in appropriate bin] H26@ [mult reg := mask to isolate bin index] C3@ [acc := top 4 bits of value] R1024F [12 right, get bin -8..7 in address field] A29@ [add 8 to address, bin index now 0..15] A27@ [make A order for bin] U113@ [plant in code] S30@ [convert to T order] T115@ [plant in code] [113] AF [(planted) acc := count in bin] A2F [inc by 1 in address field] [115] TF [(planted) store updated bin count] A1@ A2F U1@ [inc negative count of variates] G49@ [loop till got required number] [Print mean and standard deviation] A6#@ TD [pass mean to print subroutine] [122] A122@ G191@ O37@ O37@ O37@ [print spaces] A8#@ H6#@ N6#@ T4D [calc variance, pass to square root s/r] [131] A131@ GA [4D := standard deviation] A4D U8#@ TD [pass to print subroutine] [136] A136@ G191@ [print s.d.] O38@ O39@ O38@ O39@ [print CR, LF twice] [Print histogram] A29@ LD A27@ [make A order for exclusive end bin] T28@ [store as test for end] A32@ T33@ [initialize mid-value of bin] A27@ [A order for bin{0}] [149] T158@ [loop: plant A order for current bin] TD A33@ U1F [0D = mid-value of bin, extended to 35 bits] A31@ T33@ [update mid-value for next time] [155] A155@ G191@ [print mid-value] O26@ [print colon] [158] AF [(planted) load number of hits in address field] A29@ [add 8 for rounding] R4F [divide by 16] E163@ [jump to middle of loop] [162] O175@ [print plus sign] [163] S2F E162@ [loop till printed enough plus signs] O38@ O39@ [print CR, LF] TF [clear acc] A158@ A2F [make A order for next bin] S28@ [compare with end order] E174@ [exit if no more bins] A28@ [restore acc after test] G149@ [loop back (A order is negative)] [174] O34@ [dummy character to flush print buffer] [175] ZF [halt program; also serves as plus sign] [Subroutine of main routine. Sets 0D := (1/32)log_2 of uniform variate] [176] A3F T188@ [plant return link as usual] [178] T4D [pass range = 0 to PRNG] [179] A179@ G1G [0D := uniform variate] AD S2#C [test for too small (unlikely)] G189@ [jump to try again if so] A2#C T6D [OK, pass uniform variate to logarithm subroutine] [186] A186@ GL [0D := logarithm] [188] ZF [(planted) return to caller] [189] TF E178@ [if failed, clear acc and try again] [Subroutine of main routine; prints signed fraction in 0D to 4 decimals.] [Wrapper for library subroutine P1; adds sign, rounding, layout.] [191] A3F T207@ [plant return link as usual] AD G197@ [acc := value, jump if < 0] O37@ E200@ [print space, jump to common code] [197] O36@ [value < 0, print minus] TD SD [acc := abs(valus)] [200] A8#C TD [add 0.00005 for rounding, pass to P1 in 0D] O35@ O30@ [print '0.'] [204] A204@ GN P4F [call P1 to print 5 decimals] [207] ZF [(planted) jump back to caller] [========================== G parameter ==========================] [Linear congruential generator, same algorithm as Delphi 7 LCG.] [38 storage locations.] [Initialize: Call 0G with 35-bit seed in 4D.] [Input: Call 1G with 35-bit range in 4D.] [Output: if range > 0, 0D := random 35-bit integer in 0..(range - 1).] [if range = 0, 0D := random 35-bit real in [0, 1)] E25K TG GKG10@G15@T2#ZPFT2ZI514DP257FT4#ZPFT4ZPDPFT6#ZPFT6ZPFRFA6#@S4#@ T6#@E25FE8ZPFT8ZPFPFA3FT14@A4DT8#@ZFA3FT37@H2#@V8#@L512FL512FL1024F A4#@T8#@H6#@C8#@T8#@S4DG32@TDA8#@E35@H4DTDV8#@L1FTDZF [==================== Library subroutines ============================] E25K TL [L1: Logarithm to base 2.] [0D := (1/32)*log_2(6D), provided 6D > 2^(-32)] [38 storage locations; working positions 4D and 8F.] GKA3FT33@E11@IFP1024FP512FA3@LDT6DADS4@TDS3@A6DG6@T8FS5@T4D H6DV6DS3@E34@A3@LDYFT6DA4DADTDA4DRDYFG17@EFA3@YFT6DE29@ E25K TN [P1: Print non-negative fraction in 0D, without layout or rounding] [21 storage locations.] GKA18@U17@S20@T5@H19@PFT5@VDUFOFFFSFL4FTDA5@A2FG6@EFU3FJFM1F E25K TA [S2: square root.] [Closed: 22 storage locations, working position 0D.] [Forms sqrt( C(4D)) where C(4D) > 0, and places result in 4D.] GKA3FT20@A4DS9@A6@UDHDR1FS21@TDN4DRDA4DYFT4DVDTDVDYFG5@EFSF [======================= M parameter again ======================] E25K TM GK E40Z [define entry point] PF [acc = 0 on entry]