HHC 2004 Contest Problem Revisited

Roger Hill rhill at siue.edu
Sat Nov 27 15:06:37 PST 2004


Hi all,

The Programming Contest that Joe Horn presented at HHC 2004 raises some
interesting theoretical questions in addition to being a fun programming
challenge.  A couple of weeks after the conference, while I should have been
doing some exam grading, I happened to noticed my entry still sitting in my
HP-49G+ and between thinking more about the problem and grading the exam,
guess which won out...  Anyway, I came up with a program just 77 bytes long.
Of course it is long since too late for contest purposes, but inspired by
Jeremy Smith's interesting article(**) documenting and discussing the
various contest entries, I thought I would write it up and share it with you
along with some further notes on some of the mathematical concepts involved.

To remind you of the basic requirements:  "The program must graph a large
circle (at least 60 pixels high for all HP48 & HP49 models), then plot
random pixels inside that circle until approximately half the pixels are
dark, then exit without user intervention, leaving the graph displayed. ...
Assume machine default flag settings (except HP49; assume RPL mode).
Altered flag settings must be returned to default status upon program
completion.  Stack contents before and after the program must be the same."

I'll start off by giving the program and then explain it.  Some of the
explanations apply to some of the other contest entries as well.

  << "4!2E3(0,0)3 6!0{#0#0}0" OBJ-> FREEZE PVIEW ARC
     START RAND SQRT 3 * -1 RAND ^ SQ * PIXON NEXT >>

(HP-49G+ checksum #3153h, 77 bytes).  Of course SQRT means the square-root
symbol.

The initial text string, unpacked by the OBJ-> command, gives the various
numerical parameters for the later commands.  0 FREEZE is to keep the
display preserved after the proram has terminated, {#0h,#0h} PVIEW displays
the graphic, and (0,0) 3 6! 0 ARC draws the circle as an arc centered about
(0,0) of radius 3 units (= 30 pixels with the calculators default settings)
from angle 720 to angle 0.  The calculator seems to be smart enough to not
waste time retracing itself after a complete circle has been drawn, so there
is no functional disadvantage in giving an angle larger than necessary;
using 6! saves some bytes and works in any of the angular modes.  If I were
actually going to *use* the program and not worry about byte count, I would
put ERASE at the beginning of the whole program to get rid of any previous
graphic, but this was not required by the contest rules.

The heart of the program is the algorithm between START and NEXT which turns
on random pixels in the circle.  The initial and final parameters 24 and
2000 for the loop give 1977 points in all -- more about this number later.
A point is plotted by giving its coordinates as a complex number and
executing PIXON.  Since the points must be in the circle, it is natural to
use polar coordinates (r,theta):

  (x,y) = x + i*y = r*cos(theta) + i*r*sin(theta) = r*exp(i*theta).

To get a uniform distribution of points in the circle, it is clear that the
angle theta should be uniformly distributed between 0 and 2*pi radians, but
what about the radius r?  If r were uniformly distributed between 0 and 3,
the points would end up more concentrated near the center, because for given
angles, points with smaller r are closer together than points with larger r.
If we want a uniform distribution, the average number of points with
coordinates between r and r + dr and between theta and theta + d(theta)
should be proportional to the area of that region, which is approximately a
rectangle of sides dr and r*d(theta) if dr and d(theta) are small:

  dN = K * d(area) = K * dr * r*d(theta) = K * r * dr * d(theta)

where K is a constant.  Now this can be written as

  dN = (K/2) * d(r^2) * d(theta).

This means that if r^2 is uniformly distributed between its minimum and
maximum values (0 and 9), and theta is uniformly distributed between its
minimum and maximum values (0 and 2*pi), then the points will be uniformly
distributed throughtout the area.  The command RAND gives a (pseudo-)random
number uniformly distributed between 0 and 1, so we will get the desired
distribution by letting

  r^2 = 9 * RAND   i.e   r = 3 * SQRT(RAND),  and
  theta = 2*pi * RAND .

The complex number describing the plotted point is then

  (x,y) = 3 * SQRT(RAND) * EXP(i*2*pi*RAND) .

The exponential can be greatly simplified (and more bytes saved) by using
another trick based on the identity EXP(i*pi) = -1:

  EXP(i*2*pi*RAND) = (EXP(i*pi*RAND))^2 = ((-1)^RAND)^2 .

The beauty of this is that taking -1 to any real fractional power will
automatically produce a complex number on the calculator without having to
build it up from real and imaginary parts.  This seems to work regardless of
whether complex numbers are enabled on the HP-49G+.  

Consider next how many points need to be plotted.  Jeremy gave a nice
discussion of this, to which I might add a little:  Let N be the number of
pixels in the region of interest, and let M be the number of points plotted.
If we want half of the pixels (N/2) to be dark on the average, we have to
plot more than N/2 pixels because there is a significant probability that
some pixels will get "hit" two or more times.  Let p be the ratio M/N, which
is the average number of times a given pixel gets "hit" in the plotting
process.  In the limit where both N and M are large, the probability that a
given pixel gets "hit" k times -- call it P(k) -- is given by a Poisson
distribution:

  P(k) = EXP(-p) * p^k / k!

Thus the probability that a given pixel does not get hit at all is P(0) =
EXP(-p), the probability that the pixel gets hit exactly once is EXP(-p)*p,
the probability that the pixel gets hit exactly twice is EXP(-p)*p^2/2!, and
so on -- the sum of these probabilities comes out to 1 as it should be.  If
we want half of the pixels to be dark on the average, that means we want the
probability of not getting hit at all to be 1/2, i.e.   P(0) = 1/2,   or
EXP(-p) = 1/2,   or   p = LN(2) = 0.693 .  This means   M/N = 0.693,   so we
actually have to plot M = 0.693*N pixels.  All this is of course on the
average; there's no guarantee that we will get *exactly* N/2 dark pixels.  I
assumed that the phrase "approximately half" in the contest rules allows for
this; otherwise one would have to have a much more elaborate program that
keeps track of the dark pixel count.

For a circle of radius 3 units = 30 pixels, the area in pixels is pi*30^2
which is about 2827 pixels (this is N), so the number of pixels plotted
should be around M = 0.693*2827 = 1960 pixels.  Any number around that
should "approximately" do the job, and I settled for 1977 which took the
fewest bytes.  I also ignored some pixel quantizaton effects that Jeremy
considered in his article, e.g. the fact that the circle itself is 1 pixel
thick, figuring that the "approximately" would cover such slop as well.

There is one thing that this program requires and assumes:  The
pseudo-random number operation RAND must produce random numbers strictly
between (but not including) 0 and 1.  I believe is the way most
pseudo-random-number generators work -- anyone know for sure about the
49G+'s algorithm?  If RAND were to produce either 0 or 1 exactly, then -1 to
that power would fail to give a complex number, causing an error with the
PIXON command, but I have not seen this happen myself.

At any rate, it all brings up a lot of interesting probability concepts, and
I hope at least some of you find it interesting and/or useful.

See you at HHC 2005...

-- Roger

(**)  Jeremy's article being referred to is at
http://www.peak.org/~jeremy/calculators/HHC%202004%20Conference%20User%20RPL
%20Programming%20Contest.pdf




More information about the HHC mailing list