## Generating Random Numbers On Decimal Machines

Porting a phrase generator from a binary (2's complement) machine to a decimal (0-9) machine started the following -
The randomizer was considered too slow and not optimized for decimal machines -

from Michael Albaugh - Dec 20, 2019 2:57 pm
attached file
FMakRand-f.txt
 An update on my experiments, including a surprise. ... Context: ``` > A Fortran version of MakRand() would also have multiply and divide by ten in place > of the moves, in the rotates: (typed in a hurry, please excuse typos) > > IMSDS = IR1/1000 ISOLATE R1 MSDS > ILSDS = IR1-(IMSDS*1000) ISOLATE R1 LSDS > IR1 = (ILSDS*1000)+IMSDS COMBINE TO COMPLETE ROTATE 3 PLACES RIGHT > IMSDS = IR2/100000 ISOLATE R2 MSDS > ILSDS = IR2-(IMSDS*100000) ISOLATE R2 LSDS > IR2 = (ILSDS*100000)+IMSDS COMBINE TO COMPLETE ROTATE 5 PLACES RIGHT > > Two divides, four multiplies. I probably missed some optimization, but still, I feel the > major issue is "You can't say that in 1401 Fortran" (or not concisely) ``` So I decided to actually implement a small test program. As predicted above, I badly muffed the hasty transcription, but the code attached runs, and produces the same first 20 numbers as the Autocoder version. Concentrating on my third question: > 3) Can MakRand() retain its efficiency edge implemented in 1401 FORTRAN. > (I suspect not, given that there are four multiplies and two divides inherent in the > Fortran version.) So I looked at the generated code snippet of first rotate, which I annotated before findingout that the "indexing" was actually a type flag (A-bit or B-bit but not AB -> Integer): C ISOLATE R1 MSDS 7 IMSDS = IR1/1000 4745: B 700 CALL Interp 4749: DCW @6PW#6MW/GUF@ : IMSDS+X2 # IR1+X2 / CONST1000+X1 4760: DCW @|@ RETURN? C ISOLATE R1 LSDS 8 ILSDS = IR1-(IMSDS*1000) 4761: B 700 CALL Interp 4765: DCW @6OW#6PW*GUFN&6MW@ : ILSDS+X2 # IMSDS+X, * CONST1000+X1 N & IR1 4781: DCW @|@ Return C COMBINE TO COMPLETE ROTATE 3 PLACES RIGHT 9 IR1 = (ILSDS*10000)&IMSDS 4782: B 700 CALL Interp 4786: DCW @6MW#6OW*GVA&6PW@ : IR1+X2 # ILSDS+X2 * CONST10K+X1 & IMSDS+X2 4801: DCW @|@ Return C DITTO FOR R2, 5 PLACES 10 IMSDS = IR2/100000 4802: B 700 The thing that leapt out at me was that integer arithmetic is also interpreted, although we (or at least I) had assumed only FP was. The operand type is encoded in the zone bits of the tens digit of addresses (I went down a rabbit hole trying to make sense of what appeared to be +X2) So that first string, after a bit of massage because I forgot to set FORTRAN character set, and to split it up as the interpreter would see it: 6PW integer at 4676, which the listing says is IMSDS = assignment 6MW int @ 4646, IR1 / GUF int constant 1000 | (record mark) end of expression I did not really expect the compiler to recognize mul/div by powers of ten and convert to Moves, but then I did not expect integer arithmetic to be interpreted either. Anyway, this confirms that my highest priority should be figuring out how to define user functions and patch the library to include them. If I can figure that out, it would be cool to add MakRand to our library. -Mike

from Ronald Mak/SJSU - Dec 02, 2019 11:11 pm
 Don?t go below the sample size of 30, otherwise the Central Limit Theorem no longer applies.

from Robert Garner Dec 02, 2019 11:02 pm
 Ron, > Generate 1,000,000 random values and you?ll get a much nicer bell curve. Here?s the histogram with 1,000,000 bins, where each bin is average of 100 Mak-Random values:

from Michael Albaugh Dec 02, 2019 4:53 pm
 > On Dec 2, 2019, at 11:56 AM, Ronald Mak/SJSU wrote: > > Generate 1,000,000 random values and you?ll get a much nicer bell curve. How long will that take on the 1401? Hmmm. Using my rough memory of instruction timing (this is called "covering my ass about wild errors") I get about 209 clock cycles (92 for instruction fetch, 117 for data) to _generate_ one value (counting the Return, but not the Call), so at 11.5 uSec per clock, that would be 2403.5 uSec per value. So a million values would take 2403.5 seconds roughly 40 minutes. Of course the time taken for the analysis of those values is not included. If someone wants to port Python to the 1401, I'm not gonna stop you, but I'm not holding your beer. :-) Generation code snippet with times (I and D for Instruction and Data cycles): ``` * PRNG PROPER ORG 500 MAKRND SBR MRRTN&3 SAVE RETURN ADDRESS * I 5 D 3 * ROTATE R1 RIGHT BY THREE CHARACTERS. MCW R1-3,ACC MOVE NUM MSDS TO ACC LSDS * I 8 D 8 MCW R1 MOVE NUM LSDS TO ACC MSDS (CHAIN) * I 5 D 6 MCW ACC,R1 * I 8 D 14 * ROTATE R2 RIGHT BY FIVE CHARACTERS MCW R2-5,ACC MOVE NUM MSDS TO ACC LSDS * I 8 D 4 MCW R2 MOVE NUM LSDS ACC MSDS (CHAIN) * I 5 D 10 MCW ACC,R2 * I 8 D 14 MCW R2,SUM * I 8 D 14 A R1,SUM * I 8 D 14 MZ @0@,SUM-6 IGNORE OVERFLOW * I 8 D 2 MCW R2,R1 * I 8 D 14 MCW SUM,R2 * I 8 D 14 MRRTN B MRRTN RETURN TO CALLER * I 5 ```

Ronald Mak/SJSU - Dec 02, 2019 11:56 am
 Generate 1,000,000 random values and you?ll get a much nicer bell curve. How long will that take on the 1401?

Robert Garner wrote December 1, 2019 at 10:32 PM
 Ron, I slightly modified your program, attached, to run under Python 3.2 (my MacBookPro?s default), and to use the standard matplotlib histogram routine (not the prettier seaborn plotting library, which I need to learn. :) Here?s the resulting histogram of your sum of 100 ?Mak Random? variables ? same as yours. :) My plan is to try higher accumulated sums and then next run the Diehard random generator test suite. :)) This is fun exploring your efficient random number generator for 1400-class computers. :) - Robert

An e-mail from Ron Mak dated Nov 29, 2019
Subject: Normally distributed random numbers on the 1401
From: Ronald Mak/SJSU (Add as Preferred Sender)
Date: Fri, Nov 29, 2019 2:09 pm
To: Robert Garner
Cc: Marc Verdiell , Ed Thelen , Carl Claunch , Ken Shirriff , W Van Snyder , Mike Albaugh , Cay Horstmann , Steven Moffitt , Luca Severini
Hi, random fans.

I hope you all had a great Thanksgiving! Please enjoy some of my leftover turkey:

I hacked together a 1401 Autocoder program that generates, prints, and graphs 10,000 normally distributed random numbers (i.e., they fit a bell curve). As in my previous program, I use my double-rotation algorithm to generate uniformly distributed random numbers. In this program, I keep a running sum of the last 100 uniform random numbers. I accomplish this by maintaining a circular buffer of the last 100 numbers. Each time I generate a new uniform random number, I subtract the oldest number in the buffer from the sum, add the new number to the sum, and enter the new number into the buffer to replace the oldest. According to the Central Limit Theorem of statistics, all the averages of the individual sets of 100 buffered numbers will approach a normal distribution. I chose sample sizes of 100 for the buffer because on the 1401, you can divide by 100 is by simply ignoring the last two digits of the sum.

I?ve attached the Autocoder program and the Python program for comparison. The outputs and graphs of both programs match, so either I got them both right or they?re both wrong. To run the Python program in a terminal window:

ipython -i MakNormG.py

The ?G? is for ?graph?, which opens in a separate window.

Confused by all this? No, I?m not a statistician. I don?t even play one on TV. I only teach the subject ...

-- Ron

P.S. I noticed a problem with ROPE, which I believe some of you have encountered, and perhaps my former top student Luca has fixed. If you accidentally send the 1401 simulator into an infinite loop, and you attempt to recover by aborting ROPE, the simulator continues to run! ROPE forks the simulator as a separate Java thread, and normally, it shuts down the thread when your 1401 program completes. But if you kill ROPE with an abort signal (control-C, force kill, or whatever), you will have to kill the simulator process manually (e.g., kill -9 i1401r). Otherwise, the zombie process holds onto its files, and when you restart ROPE, you will have problems accessing those files.

Ronald Mak
 http://www.cs.sjsu.edu/~mak/ Department of Computer Engineering Department of Computer Science Department of Applied Data Science San Jos? State University One Washington Square San Jos?, CA USA 95192 ?On two occasions I have been asked, ?Pray, Mr. Babbage, if you put into the machine wrong figures, will the right answers come out?? I am not able rightly to apprehend the kind of confusion of ideas that could provoke such a question.? -- Charles Babbage, 1791-1871 World?s first computer scientist ?The purpose of computing is insight, not numbers.? -- Richard Hamming, 1915-1998, Mathematician
Attachments included:
MakNormG.s, MakNormG.lst, MakNormG.out, MakNormG.py, Screen Shot 2019-11-29 at 2.08.31 PM.png

from Mike Albaugh, Nov 11, 2019 5:52pm
 > Mike, > > Did you get a chance to corner/ask Knuth last weekend if he knew of an efficient random number > generator for a decimal/character-oriented computer like our 1401? Yes and no. I mentioned the project and that I had cracked open his book for some tips, but of course blew my credibility when I mis-remembered which volume. He correctly said "Seminumerical Algorithms", which made sense, but of course my memory of taking it off the shelf said "Fundamental Algorithms". When I got home I looked at the book itself sitting on my desk, and Surprise, Surprise he was right. Anyway, I suspect that section on Linear Congruential in the _right_ book (that is, shelved to the right of the one I mis-remembered :-) has plenty of tips, although it wouldn't be a textbook without having the student work exercises leading to enlightenment. I see you have included Van, which is great, because I suspect the real slog-work is going to be the process of adding a RNDOM() subroutine to the 1401 library, creating a new compiler tape that includes it. MY plan, such as it is, is to try to do _any_ subroutine and call it from FORTRAN. Compared to that, the RNDOM subroutine itself is fairly simple. The one in the example code is a classic Linear Congruential, of the form (expressed in C: ``` unsigned int seed; unsigned int rndom() { seed = ((seed * A) + C) % M); return seed; } ``` with some faffing around to convert the range from 0 to 65535 to the range 0.0 to 1.0-epsilon, for compatibility with the 7090 version. For my original version, the constant M is 65536, because the original version was in fact written in C (and simultaneously 68000 assembly). One of the few things C guarantees about arithmetic is that "unsigned int" will be evaluated as an actual two's complement number, so M = 65536 can simply be a mask with 0xFFFF. IT is common with Linear Congruential PRNGs to pick an "easy" M like this. But of course on a decimal machine, that is not an easy modulus. A and C must be chosen (primarily having no prime factors in common with M) to avoid issues with sequence length. Seminumerical Algorithms has a section that is a great start, along with the exercises. But of course when attempting the choice of A and C (Given M) one should probably start with some way to evaluate the results. Given the vast amount of memory and speed that modern computers need to properly render Cat Videos, they are well equipped to actually do these tests relatively quickly, even for a nine-decimal digit (to contain the equivalent of the 7090's mantissa, and then some. So, I'd love to see someone come up with appropriate tests, either in general or "good enough" that after reducing to 0..11, the number of occurrences of each value are reasonably equal, and that sequences (R1, R2, R3) from the output are also reasonably equally distributed. Beware that human generated attempts at random numbers often fail this test, because they are light on "can't happen" sequences (e.g. 1,2,3 or 7,7,7) that are just as likely as any of the other 1000 (or 1728 :-) Note that my existing generator is a bit shy of 6 digits (there's a reason many computers had 18-bit integers until the fetish for 8< > * > > I?ve asked Ron Mak if he has any ideas. > (As he had authored The Java Programmer?s Guide to Numerical Computing). > Ron ?> Perhaps you might assign the question to one of your more engaged students. :)) Thanks!! :))) Hey, I'd love so help, at least with testing, and if someone is ambitious enough to find some good A and C, it could be fun. Meanwhile, I need to get back to trying to understand the 1401 Fortran subroutine library/linking process. BTW: anybody getting this only now, I swear my original html did not produce "Code Blue". No idea, other than it has been quite some time since it was OK to feed code through a modern email agent and try to compile and run it. (What we got back, didn't live too long") > - Robert > > * This is your current random number generator (in 1401 FORTRAN). > When it was running last week, it felt like it took about 4 seconds per invocation/monostich. > 206 ISEED = (ISEED * 13849) + 23861 > C TAKE MODULO 65536 TO MIMIC ORIGINAL BINARY > ITOP = ISEED/65536 > ISEED = ISEED - (ITOP * 65536) > C CONVERT TO FLOAT AND MAP 0..65535. to 0..1.0 > RNDNO = ISEED > RNDNO = RNDNO/65535. Note that the slowness is almost certainly two divides, one floating point. Using a "machine friendly" M would not help this, as Fortran of this era did not have a shift operator (in keeping with the attempt to be radix-agnostic). That's why any decimal version will have to be in assembly language to gain much. Thanks to all, Mike

from Robert Garner, Nov 11
 Mike, Did you get a chance to corner/ask Knuth last weekend if he knew of an efficient random number generator for a decimal/character-oriented computer like our 1401?* I?ve asked Ron Mak if he has any ideas. (As he had authored The Java Programmer?s Guide to Numerical Computing). Ron ?> Perhaps you might assign the question to one of your more engaged students. :)) Thanks!! :))) Cheers, - Robert * This is your current random number generator (in 1401 FORTRAN). When it was running last week, it felt like it took about 4 seconds per invocation/monostich. ``` 206 ISEED = (ISEED * 13849) + 23861 C TAKE MODULO 65536 TO MIMIC ORIGINAL BINARY ITOP = ISEED/65536 ISEED = ISEED - (ITOP * 65536) C CONVERT TO FLOAT AND MAP 0..65535. to 0..1.0 RNDNO = ISEED RNDNO = RNDNO/65535. ```