Analysis
The problem I have chosen to solve is pseudorandom number generation (PRNG) using the LMC instruction set. I have chosen this problem because the LMC supports the two most basic mathematical operations- addition and subtraction of integers- which should be substantial enough to construct a basic PRNG algorithm. PRNG algorithms are incredibly important in computer security, being used to derive encryption and decryption keys and secret keys that are important to be unpredictable given data either side of the key or data encoded via the key.
A number of different algorithms are available for random number generation. The first is a more mathematical approach called Multiplicative Congruential Generators.
An LCG is an algorithm that applies a modulo function of divisor ‘m’ across a dividend consisting of a ‘seed’ that has been multiplied by a constant ‘a’. The value of this operation is then used as the next seed.
A varient of this known as a mixed congruential generator, where a certain constant ‘c’ is added on to the dividend. Both the mixed and multiplicative generators are examples of linear congruential generators (LCGs).
Another algorithm is called the middle-squares method. The middle-squares method works by squaring a number, and then stripping the least significant digit and as many most significant digits as required until you have a number the same digits as you started. This is then repeated on the output of the function to generate another number.
The method used fairly typically on Linux systems is ‘/dev/random’- this is a file that updates with random noise sampled from the real world, typically things like network timings, power fluctuations and other minute, inevitable and uncontrollable environmental variables.
Of these 3 methods, the LCG is most ideal for the LMC instruction set. Multiplication and remainder can both be performed using just addition and subtraction, whereas doing manipulation required for the middle-squares method will be considerably more difficult. The middle-squares method will also be very limited on the LMC due to the restrictions on the size of a single piece of data being low (-999 to 999). An implementation of /dev/random would not be possible, since the LMC does not allow access to monitoring required for noise collection.
The problem can be therefore broken down into 3 stages: multiplying the seed by the constant ‘a’, performing remainder division on the result of this and the constant ‘m’, and then putting the result of this back into the algorithm.
However, this first raises the question of what values to use for ‘m’ and ‘a’. Using a prime as ‘m’ is a good idea, because it allows every number less than ‘m’ to be selected given the right value of ‘a’. ‘a’ must be a primitive element of ‘m’ for the LCG to have this property. 997 is the largest prime number usable in the LCG, and 505 is the smallest primitive root of 997 that multiplies by any value (other than 1) to a value greater than 997. Therefore, our ‘m’ value will be 997 and our ‘a’ 505.
Given a seed of 1, we should expect a period of 996 values (every number appears once, then 1 appears again and the cycle repeats), in seemingly random order. This isn’t true randomness, and a clear flaw in this system is that numbers will only appear once per cycle, however this could be resolved using another remainder operation on the result to generate numbers inside a given range:
( (a * seed) MOD m ) MOD maximum
This would generate random numbers within the given range, but may express a bias
Any seed should yield a period of 996 values before the seed reappears. The order should be pseudorandom, and this could be tested by comparing the intervals of the numbers generated. Another more visual test however could be to map the results and watch for patterns in the path produced (the first result moves upward, the second across, the third upward and so on). This visual approach makes data that is practically random or non-random easily distinguishable.
It must be tested with different seeds to check the result always appears as random.
A successful random number generator will generate numbers with an unpredictable interval between the numbers and an unpredictable pattern for repetitions below 996 (or the period length)
Implementation
Implementation can be split into 3 parts of functionality: remainder, multiplication and then iteration.
A method to calculate the remainder involves subtracting a value until the result is non-positive, and then adding one of the result on:
DIVISOR ← INPUT
DIVIDEND ← INPUT
ASSERT DIVIDEND ≥ DIVISOR
diff ← DIVISOR
diff ← diff – DIVIDEND
diff ← diff + DIVIDEND
RETURN diff
In the real implementation, the DIVISOR will be a hard-coded constant of 997. The DIVIDEND will be the seed for generation (or the last randomly generated number). Assert is a condition test: if the assertion is false, then the algorithm will fail. The nature of the algorithm in execution ensures that the assertion will only fail when the seed is 1, which is a none-issue, however it’s important that it is clear that the dividend should be greater than the divisor.
This can be implemented in an LMC with the following logic:
start LDA mod_r
BRP end
BRA mod
end OUT
HLT
mod LDA mod_a1
mod_p SUB mod_a2
BRP mod_p
ADD mod_a2
STO mod_r
BRA start
mod_a1 DAT 997
mod_a2 DAT 5
mod_r DAT -1
This algorithm is identical to the pseudocode, however it is missing the assertion present in the pseudocode, since I believe it to be unnecessary for the implementation of the algorithm.
The next algorithm required is a method to multiply, since the LMC doesn’t natively implement this. Integer multiplication can be replicated by repeated addition of values, or the sum of a set of length ‘n’ containing only the value ‘m’ (equal n ✕ m). This can be replicated in pseudocode as below:
M ← INPUT
N ← INPUT
OUT ← 0
COUNTER ← M - 1
WHILE COUNTER ≥ 0:
OUT ← OUT + N
COUNTER ← COUNTER - 1
RETURN M
‘M’ and ‘N’ are both arguments, so the result returned is ‘M ✕ N’. I implemented this in a similar fashion to the remainder function, using a number of memory spaces as arguments to the function to make it easier to reuse:
start LDA mul_r
BRP cont
BRA mul
LDA mul_r
OUT
HLT
mul LDA mul_a2
STO mul_c2
LDA mul_a1
SUB mul_s
STO mul_c1
mul_p LDA mul_a2
ADD mul_c2
STO mul_c2
LDA mul_c1
SUB mul_s
STO mul_c1
BRP mul_p
LDA mul_c2
SUB mul_a2
STO mul_r
BRA start
mul_a1 DAT 4
mul_a2 DAT 10
mul_c1 DAT -1
mul_c2 DAT -1
mul_s DAT 1
The function is broken into 3 main parts, denoted by the bold statements. The first part is the preparation sequence. The arguments are loaded from memory and cloned into new positions mul_c2 and mul_c1. This allows us to edit the values in place and still have access to the original values for use later. 1 is subtracted from the first argument before it is stored.
The function then enters the iterative stage. The second argument is loaded into the accumulator and added once to the cloned first argument. Then, the cloned first argument is loaded and reduced by one using the constant value mul_s. This is then used to overwrite the previous value. Supposing that the cloned multiplier is still above or equal to zero, the function branches back to the top of the stage and repeats.
Once this is done, one more operation must be performed (since the counter ranges down to -1 not zero) – the second argument must be subtracted once to ensure that it’s N ✕ M and not N ✕ (M + 1). It then writes this to the mul_r return space and branches to the top of the program to exit the function.
During testing I noticed an error in the function and the potential for a simpler approach to the final stage:
start LDA mul_r
BRP cont
BRA mul
cont LDA mul_r
OUT
HLT
mul LDA mul_a2
STO mul_c2
LDA mul_a1
SUB mul_s
STO mul_c1
mul_p LDA mul_a2
ADD mul_c2
STO mul_c2
LDA mul_c1
SUB mul_s
STO mul_c1
BRZ mul_end
BRA mul_p
mul_end LDA mul_c2
STO mul_r
BRA start
mul_a1 DAT 4
mul_a2 DAT 10
mul_c1 DAT -1
mul_c2 DAT -1
mul_s DAT 1
mul_r DAT -1
The fix was that a cont (short for continue) tag at the top, for the program to jump to. Missing this tag caused an assembling error.
The other change was to switch the BRP tag for a combination of a BRZ and BRA. This enables the program to run without requiring the subtraction at the end of the program.
In writing this, another change was made to the end of the algorithm. It occurred that with the limited memory space on the LMC (and that each instruction takes a set time to complete), the program could be shortened and sped up by removing mul_c2 in favor of just having mul_r.
The algorithm now looks like so:
; -- snip --
mul_p LDA mul_a2
ADD mul_r
STO mul_r
LDA mul_c1
SUB mul_s
STO mul_c1
BRZ start
BRA mul_p
; -- snip --
It is faster and shorter than the original algorithm. It also removes the variable mul_c2, reducing memory usage.
I found another set of improvements to make to the algorithm. In the implementation of the PRNG algorithm, both arguments don’t need to be retained after the multiplication gets done (it’d be better if it wasn’t kept so that it can’t be followed in a debugger). The new version looks like the following (one argument is retained- the ‘a’ in the algorithm):
mul LDA mul_a2
STO mul_r
mul_p LDA mul_a1
SUB mul_s
STO mul_a1
BRZ start
LDA mul_a2
ADD mul_r
STO mul_r
BRA mul_p
I moved the BRZ up and rearranged the two processes that occur in the mul_p. This means that the preparation section can be shortened down to only 2 lines which consist of cloning the value from arg2 into the mul_r. mul_c1 is now not used.
Now is just the initialization code. Some pseudocode is available below:
rand ← 1
LOOP:
part1 ← MULTIPLY(rand, 505)
rand ← MOD(part1, 997)
YIELD rand
The yield keyword is for returning values as the function proceeds. In the real implementation, an OUT can be used.
Here’s the initialization code I wrote:
start LDA rand
STO mul_a1
BRA mul
mul_e LDA mul_r
STO mod_a1
BRA mod
mod_e LDA mod_r
OUT
STO rand
BRA start
To implement this, some small modifications were made to the other functions. mod now jumps to mod_e and mul jumps to mul_e.
The data values were also set correctly. mul_a2 was set as 505 and mod_a2 was set as 997. This produced valid output, but was too slow so I altered the values to 97 and 40. However, this produced in cases negative results. I therefore had to add a method to ensure numbers were 0 > x > m. I implemented the following:
start LDA rand
STO mul_a1
BRA mul
mul_e LDA mul_r
BRP cont
ADD big
cont STO mod_a1
BRA mod
mod_e LDA mod_r
OUT
STO rand
BRA start
big DAT 999
This doesn’t allow for every number to show up during the loop however. Tests were conducted to find the numbers that provided the best ranges. 17 and 11 provided a full set of numbers. I figured that the issue is due to the negative numbers. 17✕17 is less than 999. The root of 999 is 31.607, so therefore a prime less than this should provide a good range if given a good ‘a’ value. 31 is a prime, and 17 should be a sufficient ‘a’ value. This provides a good and completed range of values. Before proceeding to try and improve the program, I compacted the code down, copy-pasting the functions into the place that they sit in the main function:
start LDA mul_a2
STO mul_r
mul_p LDA rand
SUB mul_s
STO rand
BRZ mul_e
LDA mul_a2
ADD mul_r
STO mul_r
BRA mul_p
mul_e LDA mul_r
mod_p SUB mod_a2
BRP mod_p
ADD mod_a2
OUT
STO rand
BRA start
rand DAT 1
mod_a2 DAT 997
mod_r DAT -1
mul_a2 DAT 505
mul_s DAT 1
mul_r DAT -1
This is the shortest I managed to make the program. For testing, I set the mod_a2 as 997 and the mul_a2 as 505.
I then tested it with a number of different seed values (and replaced the BRA start with a HLT). Here’s a table of the results for the first 12 seeds.
Seed | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 9 | 10 | 11 | 12 |
Value | 505 | 13 | 518 | 21 | 526 | 34 | 539 | 42 | 547 | 55 | 560 | 63 |
Change | +504 | +11 | +515 | +17 | +521 | +28 | +532 | +36 | +538 | +45 | +549 | +51 |
These results show a clear pattern. Ignoring the hundreds, which alternate between 0 and 5, we get +8, +5, +3, +5 repeated between the results. However, in a real example the program would be reseeded with the last value it produced:
Seed | 1 | 505 | 155 | 314 | 649 | 911 | 285 | 2 | 13 | 568 | 983 | 663 |
Value | 505 | 155 | 314 | 649 | 911 | 285 | 2 | 13 | 568 | 983 | 663 | 982 |
Change | +504 | -350 | +159 | 335 | 262 | -626 | -283 | +11 | +555 | +415 | -320 | +319 |
These results appear much more random than the set that was seeded sequentially. The intervals are unpredictable between two numbers, and the average of the numbers generated is 463.92, approximately half of 997. I displayed the approximate ‘randomness’ of these results using a box plot, however evidently there is not enough data to form a conclusion, therefore a better result set is provided in the evaluation section below.
Evaluation
The algorithm can be rewritten in very little Python:
This defines almost the same algorithm as was implemented above. However, different outputs are produced due to how the LMC deals with numbers that are greater than 999 (in Python, 999 + 1 doesn’t become -999)
‘r’ is a function that takes a single argument ‘x’ that it uses to produce a random number. ‘m’ and ‘a’ are set at the top as constants, and ‘s’ is the seed/random number produced.
To produce graphs to visually represent the variance of these results, I used R and ggplot2, a graphing library available for R. I generated 4000 results and then used the modulo operator to set these between 0 and 10. This made it possible to produce histograms that show how evenly distributed the results are of the algorithm:
The results are shown to be random through how the boxes are all approximately in line with each other. I then performed a test for first-
gap correlation by lagging the results from each other:
Since the points form no pattern, it implies that there’s no correlation between a single result and the next result. The code used to produce these graphs was like so:
Exam Style Q:
Two different random number LCG algorithms with ‘m’ values of 997 produce the following sequences:
1 | 505 | 155 | 314 | 649 | 911 | 285 | 2 | 13 | 568 | 983 | 663 | 982 |
1 | 987 | 161 | 155 | 143 | 120 | 12 | 927 | 160 | 154 | 144 | 120 | 12 |
State the next 2 terms of the second sequence.
927, 160
The second algorithm is modified to produce the following sequence:
1 | 527 | 361 | 455 | 343 | 322 | 415 | 360 | 554 | 344 | 417 | 399 | 464 |
Based on the terms shown, which of these sequences suggests a better random algorithm? Give reasons why.
The first, since the IQR spans approximately half the range of the resuts. On the second algorithm, the IQR is much tighter.
Part of the algorithm responsible for the top sequence looks like so:
f1 LDA a1
f2 SUB a2
BRP f2
ADD a2
STO res
a1 DAT 505
a2 DAT 997
res DAT -1
What is the meaning of the opcode BRP?
Jump to the specified instruction if the accumulator is holding a positive value.
Describe the process this algorithm performs.
The algorithm calculates the remainder of a1 and a2, and then stores it in the position labelled ‘res’