Pi day ra-pi-dly approaching!
A while ago I discussed randomness with a friend. We had a discussion on what is random. There are so many different ways to generate random values, but in the end, they sort of turn out to not be that random after all. I though this was one of the reasons to my faulty pi in 1987.
Anyway. In this post, I just wanted to refresh the minds of you oldies who have worked with an ABC 80. And I am not an expert in arithmetics – as you can see below – but I found some interesting behavior of some graphs from a simple MATLAB code.
In some digital computers PRBS generators are used, a shift-register with a feedback polynomial guaranteeing very low correlation between values popping out from the PRBS (pseudo-random binary sequence) generator. Others pick up data from external sources and use that to generate random numbers.
Some computers employed other techniques like investigating the threshold of a diode, and more. However, even that can be a technique with flaws. For example, if the word length of the internal representation is too short, it cannot resolute a number into randomness. The quantization steps are too large.
Back in my youth I was given access to the ABC80 machines at my school. My teacher told me to go and calculate the value of pi using a method employing random numbers (the Monte Carlo method). He thought it was a nice exercise since he could not provide me with more exercises. (Or he just wanted to get rid of me).
- Generate two numbers: p and q along their corresponding axes.
- Check the result of p^2 + q^2 = r^2
- If r^2 < n^2, then the tupel x and y describe a point (p,q) inside a circle.
- Count the number of points inside the circle – it will result in a sum that is proportional to the area of the circle, i.e., proportional to pi. If we count the number of times it is inside and divide it by the total number of random tupels – and multiply by four it should “equal” pi.
The figure below says it all. It shows two vectors point at the tupels. One being inside (“Bra” = good) and one being outside (“Blä” = bad).
The problem, however, was that even though I ran many, many times and refined my code (well, “refine”… it was very BASIC – pun intended), I didn’t really get that close to pi. Why not? Well, we did not really solve the problem very accurately back then – n00bs as we were. But nevertheless it has nagged me since.
Since it popped up as a good story the other day, I decided to revisit the problem again. Unfortunately I did not find a good emulator for ABC 80. (TBD!) I also thought of the recent work by the Siliconeer on test signal generators for DACs. He has been working on deriving the requirements on number of bits and arithmetic to generate sinusoids with high spectral purity. Not too far from what the monte-carlo pi is about.
Further on, someone posted on the internet that:
“I am not 100% sure but i believe the ABC80 (remember?)
had a “real” hardware-based random number generator
built in based on a reversed-biased diode.”
This contradicted my previous assumption (my scapegoat! ) that they had a limited-length PRBS generator which could not provide me with enough accuracy. Must be something else. Something more trivial. Rudimentary? And it probably was. Maybe this internet quote could be a clue:
“I had one. I was 12 years back then.
If you asked for sqr(4) it said 2.002!”
See the graph below. I have indicated a discrete set of points that can be obtained inside the area given a certain resolution. In this case I have chosen eight levels, i.e., three bits to illustrate the available resolution. Let us ignore word and byte lengths in the ABC 80 for a moment and just look at the picture in front of us.
For this low number of levels – why bother running the Monte Carlo analysis at all? The purpose of Monte Carlo is that we want to converge towards a value with less number of simulations than would be required by exploring the entire design space, i.e., without running exhaustive testing. In the continuous-time, continuous-amplitud case, this would take forever, i.e., Monte Carlo would make sense. In the digital domain it is theoretically doable to test all cases.
For example, we can count the number of points inside: 52 out of 64. According to theory, it should then correspond to pi / 4. And indeed if we take 52 / 64 x 4 = 52 / 16 we get 3.25. not exactly pi, but it could be good enough in most cases. In some cases they could even be politics or important dates.
Massaging this simple algorithm and assuming that we have developed the “best” algorithm. An algorithm that ignores things like positive and negative numbers, extended word lengths to handle the squaring, summing, and comparison (we will soon see that they can simply be added to the list of parameters). Then one can conclude that the following must hold:
pi_approx = floor(x*pi/4)*4 / x;
where x = 2^W is the number of levels available. It is the best choice given that we can run our Monte Carlo forever (or just exhaustive testing of 2^(2*W) points… It is the best choice.
The flooring (to get an integer) comes from the inequality in the original equation of the algorithm and that we simply count the number of discrete points inside the circle – don’t worry we will also check what happens if we round upwards to pick the one slightly outside the circle.
The following MATLAB (octave) code does the job to present the approximation of pi as a function of the number of levels x = 2^W:
lambda = pi/4;
W = 32; x = 2.^(1:W);
pi_approx = floor(lambda*x)*4 ./ x;
k(1) = semilogy(log2(x), ( abs(pi - pi_approx) ),'bx-'); hold on;
pi_approx_hi = ceil(lambda*x)*4 ./ x;
k(2) = semilogy(log2(x), ( abs(pi - pi_approx_hi) ),'rx-'); hold off;
t(1) = title('Deviation from pi for different algorithms and number of bits');
t(2) = xlabel('Deviation'); t(3) = ylabel('Number of bits');
The plot below is the result and it shows the absolute deviation from the “real” pi (whatever internal accuracy my computer allows – but it is better than an ABC80’s). We see the ceiling actually does a better job now and then due to the rounding effects. There are also stair-case effects due to the same reason.
Anyways, the list of approximated pi’s are:
2.00000000000000 3.00000000000000 3.00000000000000 3.00000000000000 3.12500000000000 3.12500000000000 3.12500000000000 3.14062500000000 3.14062500000000 3.14062500000000 3.14062500000000 3.14062500000000 3.14111328125000 3.14135742187500 3.14147949218750 3.14154052734375 3.14157104492188 3.14158630371094 3.14158630371094 3.14159011840820 3.14159202575684 3.14159202575684 3.14159250259399 3.14159250259399 3.14159262180328 3.14159262180328 3.14159265160561 3.14159265160561 3.14159265160561 3.14159265160561 3.14159265346825 3.14159265346825
If, for sake of argument, we assume that the Zilog in ABC80 had 16-bits available for all arithmetic
operations, we should be able to get a pi of 3.14154052734375 according to the results above. Pretty good. However, notice that the 16 bits have to be taken into account for the sum it self. We have to look at the case when p^2 + q^2 = 2^W overflows. Once it overflows, the point is outside the range and thus not inside the circle.
Using that methodology we still have to guarantee that the product p^2 itself does not overflow, i.e., the number of levels that effectively can be used by the parts of the sum quickly goes down to W/2 – and even further than that since we need to cope with that two-thingie when p and q are equal. Thus we can only play around with 2^((W-1)/2) level if we are 13-year olds newbies in the mid 1980s with an ABC 80 in front of us. An original 16-bit setup thus quickly boils down to only 7.5 bits, or 7. From the graph and the list, the best we can get is 3.125 from that end of the spectrum.
pi = 3.125 ...
That reminds me of:
a way to redefine pi politically.
In addition: we have another side effect due to the limited length with which we can record: 2^16 samples. This also implies a truncation at ~1.53e-05 for the pi/4, i.e., a truncation of ~6.1e-05 for the pi which in turn would not give us a better pi than 3.1416.
Concludingly – do not try Monte Carlo, and if you do, optimize your code, and even if you do that, do not let it follow you 29 years. There are much better ways to calculate pi, not as fun though, and a lot of effort should have been spent on optimizing the code rather than just a set of for loops.
Like the Arabs saved early European history, some Finns have kept the ABC 80 knowledge – deep down in their archives:
Quite cute code. Rest of the story is googleable.
- Download an emulator and re-run the ABC 80 code.
- Prepare for March 14 (the Pi day).
- Revisit arithmetic and algorithms