# brianbec's WebLog

## August 2004 - Posts

##### Great Care is Needed

In response to my "Invitation to Flip me off," Jeffrey Sax pointed out that the Multiplicative Congruential (MC) Pseudo-Random Number Generator (PRNG) z <- (a*z)mod p has better statistics than I had thought (see feedback to the blog entry, and a=16807, p=2^31-1). Of course (!) this MC generates 31 bits, not 32 bits of output each cycle, and those bits are evenly distributed. The 48 to 52 ratio of 0's to 1's arises from counting bits in 32-bit registers, rather than in 31-bit registers. The expected number of 0 bits in a 31-bit register is 15.5, and miscounting them in a 32-bit register results in a ratio of 15.5 / 32, which is about 48 pct. Point scored for Jeffrey! The C# System.Random behaves equally better. Both also score a cumulative 21.3 on ComScire's tests when scaled up to 32 bits via double precision as follows: (z * (2^32-1) / (2^31-2)).  This is not nearly as good as Marsaglia's or Mersenne Twister's 24+ scores, but it is not completely hopeless. HOWEVER, it does demonstrate that the integer output of MC and System.Random cannot be used without conditioning & full understanding: even a non-naive, careful examination like mine can be fooled by the apparent but untrue.

A small mystery remains. The scaled-up MC and System.Random numbers still have a bit ratio of .4947, where as Marsaglia and Mersenne are 0.5000. Mathematica code is in the feedback of "Invitation..." for anyone who would like to play around. Ideas, anyone?

##### An Invitation to Flip Me Off

One of the simplest possible applications of a pseudo-random number generator (PRNG) is to simulate a coin flip. Get your PRNG to give you one bit at a time, map 1 to 'heads' and 0 to 'tails' (or the other way, it matters not), and flip away ad nauseum.

I have now tested in TWO independent ways the Park-Miller minimum standard PRNG, which has a single, 32-bit integer state variable, z, and computes the new z as ((16807 * z) mod (2^31-1)), which is also the result -- that is, the new PRN. No matter how hard I try (and I've been blogging on this for the last month :) I can't get it to do anything other than the following:

Heads: 48 percent of the time
Tails: 52 percent of the time

over runs of up to 639 million 32-bit z's, an appreciable fraction of the period of 2-billion some odd. Ok, so now I'm inviting y'all to check this. Please, please, prove me wrong! A small warning -- if you do not have BIGNUMS, you will overflow 32-bit integers under a naive implementation: your implementation will not be correct. My blog has lots of information about how to avoid this (the Schrage trick), and you may also see Numerical Recipes and The Art of Computer Programming, both of which vet the technique. Or you can just use a programming system that has bignums.

##### Why, Why, Y, O Y?

I've worked through the derivation of the Y combinator several times, and always found it to be a slippery concept. I'll "get it" for a while, but I move on to other things, and, by the time I come back, I've got to think about it again. One reason it's so hard to grip is the amount of syntactic noise necessitated by peripheral issues like applicative order and normal order evaluation, even in very clean programming languages like Scheme. Here, I will work an exercise (9.9, in fact) from Paul Hudak's lovely book, "The Haskell School of Expression," which results in the cleanest derivation of the Y combinator I've yet encountered. I haven't plumbed all the depths of Haskell, yet, so I am not entirely sure how it sails around the applicative vs. normal issues that make Y slippery in Scheme, though implicit currying might be a large part of it. However, I'll just "take it," because it does sail around. First, here is the normal, recursive definition of mod, expressed in Euclid's antenaresis procedure (repeated subtraction). You should be able to read this even if you don't know a thing about Haskell:

antenaresis a b = if a < b then a
else antenaresis (a-b) b

Now, here is the miraculously small definition of Y in Haskell:

y f = f (y f)

and, the answer to Paul's exercise,

ant = y (\g a b -> if a < b then a else g (a-b) b)

Painful, ain't it?  Here's why it works. Let f be (\g a b -> ...), so (y f) is ((\g a b ...) (y (\g a b ...))). This is a function of two arguments, because the first argument, g, got sucked up by the fixed-point calculation y (\g a b...). Those two arguments get bound to the variables a and b, and g gets bound to the fixed point, which is a function of two arguments. This is hardcore black magic, but it works. And it seems less slippery to me: I might be able to hold on to it longer.

For you Scheme-heads amongst us, here is the Y-combinator I use in regression tests for my Scheme implementations. I don't even want to look at this for very long, let alone dissect it like the one above, but knock yourselves out if you like:

(define appY
(lambda (f)
((lambda (x) (f (lambda (a) ((x x) a))))
(lambda (x) (f (lambda (a) ((x x) a)))))))

Posted Wednesday, August 25, 2004 4:16 PM by brianbec | with no comments

##### Nailed it! (slaps his forehead)

Well, I had to get one little theorem from the literature, but I nailed the coin-counting problem without backtracking. This is pretty nice. The theorem (in Wilf's cited lecture, and in the Combinatorica book, and implied in SICP (Structure and Interpretation of Computer Programs)) requires us to look at not just p(n), the number of partitions of n, but p(n,k), the number of partitions where the maximum species is k. The theorem is that p(n,k) = p(n-k,k) + p(n,k-1). That's a little difficult to parse mathematically and a little tricky to prove, but it translates, programming-wise, into a statement so clear and blatantly obvious it had me slapping my forehead. It's nice that insights from programming sometimes feedback into math, isn't it? Ok, here it is:

Let's call a coin (or species) of value k a k-coin. So a quarter is a 25-coin.

When the highest coin species is a k-coin, the partitions are of two kinds: those where there are k-coins and those where there are no k-coins (DUUUUH). But, this is key, and permits us to write a recursion and to PROVE the result: the partitions where there are k-coins are for amounts less than n. In fact, we can write the list of partitions where there ARE k-coins as ones with 1 k-coin appended to the list of partitions for amount n-k.  Pretty cool, it works, and it allowed me to write and verify the following Haskell program in just a few minutes:

-- This is the type of the program named "pay"
pay {- amt species -} :: Amount -> Species -> Multipliers

-- These are the auxilliary types; [square brackets enclose Lisp-like lists]
type Amount      = Integer
type Species     = [Integer]
type Multipliers = [[Integer]]

-- This function lets us put zeros in front of all the lists in a list
prepends anInt []    = []
prepends anInt aList = map (\l -> (anInt : l)) aList

-- This function lets us increment the first element in a list ...
bumpFirst []     = []
bumpFirst (l:ls) = ((l+1):ls)

-- And this one lets us increment the first elements of all the lists in our list of partitions
bumpFirsts []    = []
bumpFirsts aList = map bumpFirst aList

-- Here we go. The list of partitions when the list of species is empty is a list of the empty partition.
pay amt [] = [[]]
-- If we have just ONE species at the end of a list of species, just divide the amount by the species
-- with error checking of course
pay amt (spc : []  )
| (amt `mod` spc) == 0  =  [[ amt `div` spc ]]
| otherwise             =  error "can't make change; sorry"
-- Finally, here is the meat of the little theorem. (spc:spcs) is a list comprehension that matches
-- the highest species (our k-coin) to the variable spc, and the rest of them to spcs.  Gotta love it.
pay amt (spc : spcs)
| (amt-spc) >= 0 = (bumpFirsts ((pay (amt-spc) (spc : spcs))))
++ (prepends 0 (pay amt spcs))
| otherwise      = prepends 0 (pay amt spcs)

##### A Breed Apart

So here's a venerable little problem, simple to state, tough to get
right, even tougher to program. I give you an amount, a, and a list of
currency denominations, l; you give me back

1. the different configurations or partitions, that is, the
multipliers or coin counts of each of the denominations that add up
to the given amount

2. the most efficient partition, meaning that with the least number of
coins

3. an indication whether the most efficient partition is unique

For instance, I give you a = 40 cents, and l = [25,10,5,1], that is,
quarters, dimes, nickels, and pennies; you give me back the following
31 partitions:

{[1,1,1,0], [1,1,0,5], [1,0,3,0], [1,0,2,5], [1,0,1,10], [1,0,0,15],
[0,4,0,0], [0,3,2,0], [0,3,1,5], [0,3,0,10], [0,2,4,0], [0,2,3,5],
[0,2,2,10], [0,2,1,15], [0,2,0,20], [0,1,6,0], [0,1,5,5], [0,1,4,10],
[0,1,3,15], [0,1,2,20], [0,1,1,25], [0,1,0,30], [0,0,8,0], [0,0,7,5],
[0,0,6,10], [0,0,5,15], [0,0,4,20], [0,0,3,25], [0,0,2,30], [0,0,1,35],
[0,0,0,40]}

with the first one, consisting of 3 coins, being the most efficient and
obviously unique, by inspection.

If we're allowed quarters only, that is, the denominations are the
list [25], we can't do it, so the partitions are the empty set. If the
denominations are [25,10], that is, quarters and dimes, the only
partition is [0,4]. If the denominations are [25,10,5], adding nickels
to the mix, we get the following partitions: {[1,1,1], [1,0,3],
[0,4,0], [0,3,2], [0,2,4], [0,1,6], [0,0,8]}, for seven distinct
possibilities.

Ok, so what's a general rule or procedure for computing these
configurations or partitions?

To come up with the most efficient partition, so long as there are
pennies in the mix, I came up with the following pretty little Haskell

makeChange amt [] = []
makeChange amt (den : dens) =
(amt `div` den) : makeChange (amt `mod` den) dens

This just says that when the denominations are nil, the result is nil,
otherwise, the first element of the partition is the integer quotient
of the amount with the largest denomination, and the rest of the
partition is, recursively, the same procedure applied to the remainder
of the integer division and the rest of the denominations.  Neat and
clean. Trouble is, it's not robust. It MUST have pennies in the mix to
pick up all the 'slop' -- the remainders not divisible by the next
highest denomination from pennies.  If not, it will fail when a = 40
and l = [25,10], for instance.

Well, I thought, the general solution will fall quickly to my years of
accumulated experience and treachery, to my rapier-like cleverness and
finesse, to a relentless assault of blunt-force trauma, blah, blah,
blah.

Being a physicist, the first thing to occur to me was DOT PRODUCT!
Since the inner product of the denominations and each configuration
must equal the given amount, each solution started to look like a
vector in a lattice space, where the basis vectors were "lengths"
equal to each denomination. The difficulty with this approach is that
the needed length of the vector is a "Manhattan" metric rather than a
Euclidean. So much for that idea (though, there might still be
something there).

After spending some time with the rapier, the bludgeon, and outright
treason, I concluded that backtracking in an n-ary tree-shaped data
structure was necessary, since every attempt to generate the "next"
from a "prior" resulted in duplicates. I could stuff the duplicates in
a hash table and simply refused to generate them a second time, but
that was too ugly, even for me.

I decided to sneak a peek at history, and found that this problem has
an extensive literature. Even determining the existence of a solution
is NP-hard. However, there do exist recursive procedures, for
instance, those found at http://theory.cs.uvic.ca/~cos/. They are
based on some theorems that can be found in Herbert Wilf's lecture on
integer partitions:

http://www.math.upenn.edu/%7Ewilf/PIMS/PIMSLectures.pdf

As a bonus enticement, let me say that this subject was attacked by
the famous Ramanujan, who closed it out with one of those stupendously
rococo formulas of his involving 24-th roots of unity in the complex
plane, pi times square roots of two, and derivatives of hyperbolic sin
functions!  It was one of his crowning achievements, and must be seen
to be believed. He was surely a breed apart, especially when attacking
integer partitions.

##### Negative Mods Revisited

A few blogs back, I noted that Excel and C# differ when negative numbers are given to the "mod" function. The language Haskell has support for BOTH styles, and the clarity is refreshing. Haskell has one pair, 'div' and 'mod', which truncate toward 0 (the Excel style), and another pair, 'quot' and 'rem', which truncate toward to -Infinity (the C# style).  I had blogged that both styles were acceptable, and it's nice to see some validation of that opinion. Check this out (pasted out from a Hugs98 session -- see Haskell link above)

Main> zipWith mod [30, -30, 30, -30] [21, 21, -21, -21]
[9,12,-12,-9] :: [Integer]
(307 reductions, 472 cells)
Main> zipWith rem [30, -30, 30, -30] [21, 21, -21, -21]
[9,-9,9,-9] :: [Integer]
(99 reductions, 172 cells)

Posted Friday, August 20, 2004 5:27 PM by brianbec | with no comments

##### Black Box, White Box, Grey Box

I implemented the Mersenne Twister algorithm and ran a sample of 80 million bits through both DIEHARD and ComScire's RNGMeter test suite. It passes with results comparable to Marsalia's certified bits. This was a black-box test because I don't understand the Mersenne Twister in any detail at all. But I shall be using it from now on for all casual randoms work. It's fast enough and I don't have to spend a week worrying about it every time I want to use it for ordinary work. I'll probably stitch it into my Luna/Scheme environment, which I'd like to get back to sooner rather than later.

In the white-box corner, I ran my Multiplicative Congruential (m = 2^31-1, a = 16807) through ComScire with a resounding "thud" (-- that's a "fail" ), providing supporting evidence for yesterday's DIEHARD failure.  That's a white-box test because not only did I write every byte of the code, I understand every aspect of it in excruciating detail.

In the grey box is C# "Random" class. It also fails both DIEHARD and ComScire, and the statistics look very similar to those of my MC, so I am guessing it's a comparable MC or LC (Linear Congruential).

DIEHARD continues to choke with an end-of-file error for the last few tests on my files created via MC and the C# system Random class. This is all the more mysterious, since I used EXACTLY the same sample sizes and "Console.WriteLine" code to print out the results as I used for the Mersenne Twister and I mimicked Marsaglia's file format for his "bits" files. I am suspecting some internal mungle up due to the statistics' being so bad. Larger samples don't help. Casual perusal with emacs, less, and he (hex edit) reveal no pertinent format or structure differences between the Mersenne Twister file, Marsaglia's bits files, and my MC and C# Random outputs. ComScire gets through all of them without choking, even when they fail execrably.

I may play around with the NIST tests, but their software has no "intelligent defaults," unlike the others, so I will have to study their methodology so as to do the runs.

##### Oh, the Brutality!

I found a well regarded battery of statistical tests for random-number generators. The battery is called "DIEHARD" (nyuk nyuk), and encapsulates the accumulated experience and wisdom of George Marsaglia, a recognized expert. Follow this link for his web page http://stat.fsu.edu/~geo/diehard.html.

This is a crushing, brutal series of tests that looks for patterns in dozens of ways. Anything passing these tests will look pretty random for a wide variety of applications. Marsaglia provides some "certified" random bits at http://stat.fsu.edu/pub/diehard/cdrom/. I ran DIEHARD on several of these "certified" files, and they pass. Here are the results. The basic idea is that none of the p-values should be too close to 1.000000 or 0.000000. Capisce?

Ok, so I applied it to the C# Random class output and to the Park-Miller minimum multiplicative PRN that I vetted on theoretical grounds in the last few blogs: Oh, oh, oh, the brutality! p-values of 1.000000 All Over The Place. Pain Here and Suffering There.

For some reasons I could not diagnose, DIEHARD refused to run its last handful of tests on MY data sets, encountering some kind of end-of-file error, even though my files were of the same format (so far as I could tell) as Marsaglia's. Probably some ^Z issue. I may not bother looking into it: the evidence is IN:

C#'s Random is probably multiplicative congruential: I know mine is. But neither one is anywhere near adequate for anything that requires DIEHARD-strength randoms.

##### Whiling Away the Highway Miles

If you're going to play the game of factoring license-plate numbers to
while away the highway miles with your child, here are some ideas to
make you quicker, because if your child is bright, it's going to be a
very short time before s/he is kicking your behind. First one to
report the factors correctly wins a round!

Remember, to do well at this game, we only need to learn how to divide
quickly by the 11 primes between 2 and 31, inclusive.

First, demolish factors of two by dividing them out so long as the
result stays even. Keep count, of course, because you're going to have
to report these factors at the end.

The next easiest test is divisibility by 5. If and only if the result
ends in a 5 digit, it's divisible by 5, assuming we've gotten rid of
any factors of 2.  If you've got a factor of 5, get rid of it, keeping
track, of course.

Next easiest is divisibility by 3. Keeping adding up all the digits
until you have only one digit, and if the leftover is 3, 6, or 9, then
the original number is divisible by 3.  Example: 579 -> 5+7+9 -> 12+9
-> 3+9 -> 12 -> 3, so 579 is divisible by 3. But 529 is not because
its digits add up to 7. Can you see why this technique works? Write
the original number as 5*100 + 7*10 + 9. Since 10 mod 3 is 1, and 100
mod 3 is 1, and, in fact, every power of ten mod 3 is 1, taking the
entire mess mod 3 amounts to just adding up the digits. So, if and
only if the entire mess mod 3 is zero then the original mod 3
zero. This is just a variation of the ancient technique of "casting
out nines" that people have been using for millenia to check sums.

I do not know a quick trick for divisibility by seven, so just
practice dividing by seven, it's not too hard.

There is a trick for testing divisiblity by 11. Any multiples of 11
less than 100 are obvious -- the two digits must be the same. But
between 100 and 1000, the number is divisible by 11 if and only if the
middle digit is the sum of the two outer digits. Example, 253 is
divisible by 11; in fact it equals 11*23. 781 = 11*71.  Without this
divisibility test, you might be forgiven for staring at 253 or 781 for
a few moments wondering whether they're prime.  Can you see why this
works? 11 is 10+1. The factor of ten shifts the victim over by one and
the factor of one adds the top digit of the victim to the shifted
victim, forcing the middle digit of the result to equal the sum of the
outer two digits.

Those are the easy tests that I know. If you know some
'if-and-only-ifs' for multiples of 13, 17, 19, 23, 29, and 31, please
feel free to post them in "feedback."  I only know some "ifs" to help
along.

At this point, I assume we've eliminated any factors of 2, 3, 5, 7,
and 11.

First, I test by 13 and 23 together, because the last digit of any
multiple of 13 or 23 would be the same. In other words, suppose the
last digit of the number is 1. This could be 13 or 23 times something
that ends in 7. Or, if the answer ends in 9, it could be 13 or 23
times something that ends in 3. If that answer ends in 3, it could be
13 or 23 times something that ends in 1. Finally, if it ends in 7, it
could be 13 or 23 times something that ends in 9. In practice, there
aren't very many possibilities below 1000. You could memorize them all
or just work them out in your head as you go along. If you taught your
child the multiplication table of primes including 13 and 23, you'd
better learn it too, or you're going to lose. But I'm lazy, so ...

I use 'bridge-hand' numbers to remind me of the multiples of 13. It's
quick to recall that there are 52 cards in a pack, and that makes four
hands of 13, or three of 39 plus the dummy. 39 is a great number to do
trial multiplications with, since it's one less than 40. To multiply
by 39, just multiply by 40 then subtract the multiplier. In practice,
doing trial divisions and multiplications by 13 is easy.

23 is not much more difficult because it's almost 25 and 25 is really
easy. So the multiples of 23 go like 25-2, 50-4, 75-6, 100-8, and so
on.

Likewise, I work on 19 and 29 together. They're both really easy
because they're one less than 20 and 30, respectively, so their
multiples are easy. I throw 31 in this mix, because it's one more than
30, and, if I'm testing 29, I might as well do 31 at the same time
since it's so close. Plus you can do super-quick tricks like 29*31 =
899, because it's (30^2-1): remember that (x+1)(x-1) = (x^2-1)

The toughest one is 17, but a nice factoid is that 3*17 = 51, which is
just one more than 50, which is really easy.  So, in practice, 17 is
not too hard.

To prove that any number less than 1000 is prime, you must prove it's
not divisible by any of the primes above (in fact, the proof will work
up to 1369 = 37*37 as stated last time, but we only need three digits
for the license plate game).  So just do the quick tests for 2, 3, 5,
and 11, then do the hard work for 7.  Go for 13 and 23 together, then
19, 29, and 31 together, then tackle 17.

If you really want to set some speed records, you can just memorize
the following list of primes less than 1000.  Have it on hand, anyway,
just to verify results because mistakes can be made. But don't let
your eight-year old look at it too long or s/he'll get it memorized
before you reach the next rest stop.

There are some easy standouts like 199, 499, and 599; and 277, 577,
677, 877, and 977. These are hard to FORGET, actually, but I find that
lazily working out the trial multiples makes the miles go by quickly
and keeps me awake behind the wheel.

{2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61,
67, 71, 73, 79, 83, 89, 97, 101, 103, 107, 109, 113, 127, 131, 137,
139, 149, 151, 157, 163, 167, 173, 179, 181, 191, 193, 197, 199, 211,
223, 227, 229, 233, 239, 241, 251, 257, 263, 269, 271, 277, 281, 283,
293, 307, 311, 313, 317, 331, 337, 347, 349, 353, 359, 367, 373, 379,
383, 389, 397, 401, 409, 419, 421, 431, 433, 439, 443, 449, 457, 461,
463, 467, 479, 487, 491, 499, 503, 509, 521, 523, 541, 547, 557, 563,
569, 571, 577, 587, 593, 599, 601, 607, 613, 617, 619, 631, 641, 643,
647, 653, 659, 661, 673, 677, 683, 691, 701, 709, 719, 727, 733, 739,
743, 751, 757, 761, 769, 773, 787, 797, 809, 811, 821, 823, 827, 829,
839, 853, 857, 859, 863, 877, 881, 883, 887, 907, 911, 919, 929, 937,
941, 947, 953, 967, 971, 977, 983, 991, 997}

##### Children's Games

Here is a little game I learned from an eight-year old, who devised it herself while riding around bored in the back seat of mom's car: factor license-plate numbers in your head.

If you have a child who's just learned the multiplication table 1 through 12, challenge the child to learn the multiplication table of the 12 numbers {1, 2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31}. There are just 12 of these numbers, so the new table is no bigger nor harder to learn than the old one. And all but the first number is prime, so the table gives all the information necessary to factor numbers less than 37*37 = 1369 (That's the first composite all of whose factors are greater than 31). In fact, knowing this multiplication table by rote makes it easy to factor most license-plate numbers quickly with practice, because most license-plate numbers are < 1000. All this is well within the mathematical capacity of a bright 8-year old.

The second game is to multiply challenge pairs of two-digit numbers using the formula (a+/-b)*(c+/-d) cleverly. For instance, to get 37*37 (the little number above), just rewrite it as (40-3)*(40-3), which easy to do in one's head (it's 40*40=1600 - 3*40*2=240 + 3*3=9.  Pick the a's and b's to be multiples of 10 and the rest is easy.

Your child will be amazed at his/her new powers, or at the very least, the parent will be amazed at the child's new powers.

Posted Tuesday, August 17, 2004 11:00 AM by brianbec | with no comments