December 28, 2009

On-Demand Prime Generation

Many Project Euler problems are made easier if you have a good facility for generating and detecting prime numbers. Usually when I see such a problem, I would immediately switch to Haskell, but the other day a problem took about a minute to solve (way too long), so I ported the functionality to C.

In both, I'm doing a simple trial division algorithm: a number is prime if no other number divides it. A simple glance at Wikipedia shows you only have to test up to and including the square root of the number you're testing, and the nature of primes means you really only need to test factorization by other primes.

---

In Haskell, I used a former exercise from my CS173 class that highlighted laziness as a language feature to achieve this. First, we declare our list of primes:


primes :: [Int]
primes = filter isPrime [1..]


And then we define the predicate to determine if a number is prime or not, on which our list depends:


isPrime :: Int -> Bool
isPrime 1 = False
isPrime 2 = True
isPrime n = checkFactors (takeWhile (\ x -> x <= (floor $ sqrt (fromIntegral x))) primes) n

checkFactors :: [Int] -> Int -> Bool
checkFactors [] _ = True
checkFactors (x:xs) num = (num `mod` x) /= 0 && checkFactors xs num


For those who aren't Haskellites (even those who are, since my Haskell probably isn't very pretty or idiomatic), the code is doing this:

  • The list primes is the list of all integers beginning at 1, with all the non-primes filtered out, where non-primes are determined by the function isPrime.

  • The function isPrime tests a number x for primality by testing divisibility for all prime numbers below the square root of x, where the prime numbers come from the list primes.


Huh? You might have read over that a few times trying to piece it together and couldn't, as primes depends on isPrime and isPrime depends on primes. But not to worry, this is laziness at work :)

Laziness means that Haskell will only perform the computations it needs as it needs them, so as long as isPrime only needs as many numbers as primes has already computed, this works. Alternatively, so long as primes only needs to generate numbers isPrime can check, we're golden.

Imagine a snake that is three feet long eating a foot of itself from its tail, but growing a foot and a half longer as a result. It can do this continuously and grow to indefinite size!

More concretely, suppose I want to check if 9 is prime. A call to isPrime 9 causes the following to happen:

  • The takeWhile in isPrime will take all numbers from primes that are less than or equal to 3.

  • The language will compute values for primes by running each integer through isPrime, and will do this while those values are less than or equal to 3.
  • isPrime 1 will fail, isPrime 2 will pass, and isPrime 3...

  • ... we're not sure, so we repeat the steps above for all numbers in primes less than the square root of three. This gives us the empty list, which by definition of isPrime tells us 3 is prime.

  • We return this result, and stop taking values from primes since this is less than or equal to 3. We divide 9 by 3. This fails the predicate, so 9 is not prime.


Next to shell languages, Haskell is the only language I can think of in wide use (if you consider Haskell in wide use...) that employs laziness by default, and it allows you a powerful abstraction like this. The laziness allows us to define the list of all primes to be used freely in the code as any other list, the program will only calculate as many primes as it needs, and it can calculate arbitrarily many primes. Furthermore, the primes are comparatively fast to compute, since there are no 'wasted' divisibility tests.

---

And so I solved a number of Project Euler problems in Haskell just so I can use this facility. That being said, the speed of the language itself left me wanting, so I ported a somewhat similar abstraction to C. We'll do this in the opposite order, first defining an isPrime:

BOOL
isPrime(int num)
{
  int curr_prime, index = 0;
  double limit = sqrt(num);
  while ((curr_prime = takePrime(&index)) <= limit) {
    if (num % curr_prime == 0) return FALSE;
  }
  return TRUE;
}

This relies on a takePrime function, which can be called continuously to fetch the next prime from an index value. It looks like:

unsigned int
takePrime(int* indx)
{
  unsigned int val = prime_ptr[*indx];
  if (val != UNCOMPUTED) {
    ++(*indx);
    return val;
  }
  else {
    unsigned int last_prime = prime_ptr[(*indx)-1];
    unsigned int next_prime;
    for (next_prime = last_prime + 1; ; ++next_prime){
      if (isPrime(next_prime)) break;
    }
    prime_ptr[*indx] = next_prime;
    ++(*indx);
    return next_prime;
  }
}

Where prime_ptr is an array of integers, memset to some value UNCOMPUTED (I've left out most of the header information, as well as the init() and finished() calls that make this all work).

This takePrime preserves previously computed values (this is the purpose of the first condition), so if you ever need to check the i-th index of a value you've already computed, you simply return it. Otherwise, you take the last value you've computed previously, and increment all integers above it until you find your next prime. When you do, store it in the array and return it.

The C abstraction has a number of shortcomings over the Haskell version; namely, you can't pass in an arbitrary index and receive a prime number (last_prime only checks the prime immediately behind the one you are trying to compute. If you've only computed 6 primes and you ask for the 9th, you Segfault). You also lose the list abstraction (prime_ptr is static and this set of functions is #included, so I choose not to 'export' it), and you can't calculate arbitrarily many primes (the array has a fixed size).

All that being said, changes to correct those are pretty easy to implement; I've never needed them. But this allowed me to brute-force a few problems that might have had more elegant solutions. While I hate on C pretty frequently, gcc/g++ produce some pretty slick executables.

December 14, 2009

Prologgin'

Exam season is a royal pain, but I find ways to entertain myself. My most recent delight has been Project Euler: I discovered it late last semester, and took a break until about November to tackle some Facebook puzzles. But with the gaps between commitments getting smaller and smaller, and with a desire to use some more unorthodox languages (Facebook doesn't accept Scheme or SML submissions, for example), tackling 15 or so of the easier Project Euler problems proved fitting.

I recently solved Problem 54, an unusually straightforward implementation problem (most problems require some mathematical insight, this was a straight-up write-it-out problem). In Problem 54, they provide a file with 1000 games of poker, and you must determine how many games the first player wins.

So here I saw an opportunity: Prolog! Wouldn't it be more fun to just declare the rules of Poker and say "go," rather than hard-code every individual evaluation possibility?

-----
For those who don't know what Prolog is, here's a quick (very, very brief) primer: Prolog is a language that attempts to satisfy truth clauses given a set of relations, and rules that govern them. So if I were to define the following relation R:

(a R b) is true if and only if a is the reverse of b.


Then we know that ("paul" R "luap") is a true statement, but ("paul" R "robert") is false.

In Prolog, you supply the language with any relations you like (such as R, above), and when you provide data to those relations it can tell you if they are true or not.

Yawn. Here's the kicker: In Prolog, you can supply variables to the relations, and the language will bind those variables to values that satisfy its truth, without having to specify exactly how to find such a value. Using the example above, if I fed Prolog

x R "cantankerous"


Prolog returns

x = "suoreknatnac"


The key to this is that I never told Prolog how to reverse a word, I simply declared that R is true when one side is the reverse of the other. This becomes more clear as I go through the poker example: when you program in Prolog, you aren't thinking "what instructions can I give a machine primarily using destructive memory updates such that I can compute my goal?" (generally what happens in imperative programs), or "what functional abstractions can I define such that one's output evaluates to my goal?" (functional programs), but rather "what is my goal?" (logical programs, of which Prolog is the most popular language).

----

I'll throw a little syntax in before we dive into the example. In Prolog, square brackets define a list, so [] is an empty list, [paul, robert, annalisa] is a list with my siblings and I, and [[galosh, wader],[souvlaki, moussaka, gyro],gawker] is a list of two lists and "gawker."

In my poker program, each card is a list of the cards value and its suit (e.g. [ace,spades] or [8,diamonds]). A list of five cards is a hand, and a list of two hands is a game between two players.

Throwing it all together, here's are some of the rules I wrote that determine a player's hand in poker:


determine_hand([[_,X],[_,X],[_,X],[_,X],[_,X]], flush).


In English, this says: the relation determine_hand is true if two conditions are met. The first value is a hand of cards, for which the second value on every pair (the suit) is the same value X. Second, the second value of determine_hand is the value "flush." The underscore in the place of the values of the cards tells Prolog "we don't care what goes there," since getting a flush is only dependent on the suits of the cards. Here is another:


determine_hand([[10,X],[jack,X],[queen,X],[king,X],[ace,X]], royal_flush).


This clause says: the relation determine_hand is true if two things occur: the first, its left side is a 5-tuple of pairs. The values of the represented cards must be 10, jack, queen, king, and ace; the second value (the suit) for each card must all be the same value X. Secondly, the right side of determine_hand must be the value "royal_flush."

So if I then prompted Prolog with:


determine_hand([[10,clubs],[3,clubs],[8,clubs],[queen,clubs],[6,clubs]], HandType).


Prolog would search the possible values for HandType (variables begin with capital letters) until it found some value to make it true given the rules I've provided above. We see that all suit values are the same ("clubs"), so Prolog replies:*


HandType = flush.

-----
Like any evaluator of conditional statements, the relation rules can be chained together with standard boolean operators. The following mean:

  • left:-right means left is true if right is true.

  • Commas mean logical AND.

  • Semicolons mean logical OR.

  • And, of course, you can group with parenthesis (AND gets higher precedence than OR).


This should be enough (coming at you very fast!) to give you a flavor for how the program worked. Here is the top-level relation:


winner(H1, H2, Winner) :-
   sort_hand(H1, Sorted_Hand1),
   sort_hand(H2, Sorted_Hand2),
   determine_hand(Sorted_Hand1, X1),
   determine_hand(Sorted_Hand2, X2),
   beats(X1, X2, Verdict),
   (Verdict = X1, Winner = H1;
    Verdict = X2, Winner = H2;
    Verdict = tie,
    tiebreak(X1, Sorted_Hand1, Sorted_Hand2, SortedWinner),
     (SortedWinner = left, Winner = H1 ;
      SortedWinner = right, Winner = H2)).


It goes something like this in English: The winner relation is true if the following are true, for two poker hands H1 and H2, and some value Winner:

  1. Sorted_Hand1 is bound to the first value that makes sort_hand(H1, Sorted_Hand1) true. In this case, it's true when a hand is sorted by ascending value, so if

    H1 = [[4,spades],[king,clubs],[9,hearts],[3,diamonds],[9,spades]].

    the predicate becomes true if

    SortedHand = [[3,diamonds],[4,spades],[9,spades],[9,hearts],[king,clubs]].

    Which Prolog will find for us ^_^

  2. Find the same value for the second Hand.

  3. X1 and X2 are the values of determine_hand for the Sorted_Hands, and these are bound to the appropriate types of hand in poker (e.g. "two pair" and "full house").

  4. beats is true when the third value (in this case, Verdict) is the winner of the first two hands. If the two hands are the same, Verdict becomes "tie."

  5. One of three things is true: Either Verdict is X1, and Winner is H1 (the left player has a higher hand), OR Verdict is X2, and Winner is H2, OR there was a tie, in which case tiebreak becomes true if it's fourth value is the winner of two hands of the same type. (the problem guarantees there's always a clear winner, so we won't have actual tie games, such as two royal flushes).


----

It's a pretty radical departure from more traditional ways of programming. For those interested in logic programming, there are some great chapters near the end of PLAI on the subject. Also, The Reasoned Schemer is essential, and the 'logical Scheme' they use is an interesting counterpoint to Prolog.

I firmly believe you should base your language choice on the problem you're trying to solve, and you shouldn't contort your problem to fit the language (this is why Lisp is so much fun). I don't run into too many problems where Prolog is the answer, but when I always greatly look forward to when I do ^_^

--
Prolog links:

* = If, given the definitions, you then passed determine_hand(X,Y), where X was a royal flush, Y would be bound to "flush." Why? Because the first predicate we defined was successful (a royal flush is just a more specialized flush, and Prolog saw that determine_hand was true for the flush first). How do we get around this? You can either by ordering your clauses appropriately, or explicitly stating that one clause can only be true when the others are false.

November 14, 2009

Reading List

Almost as contentious (maybe more) than asking a programmer what their favorite editor/language is would be to ask them what the most important books are for a budding programmer to read. Every smart programmer I know has a small list of books that are essential (they claim) to know anything important about programming, software, etc.

Taking the union of all such books would yield a stack taller than you are. I'm doing my best to ferociously play catch-up.

I tend to read books 2 or 3 at a time: one is usually about a specific technology, one higher-level about programming/coding as a whole, and one usually more for fun (if I'm lucky, I get to read fiction! That normally happens most in the summer).

Because broadcasting life is more fun than living it, here's what I've been reading, and what's next in queue.

Read last year or so...

These are what I've read since starting school, Fall 2009.

If you have any suggestions, let me know!

The list of books in my room that I haven't started yet (or aren't actively finishing up) is far larger, and I'll get to those next. Wish me luck!

November 11, 2009

New Blog

Hi all,

In a financial confusion (new cards invalidate old ones, internet companies), the former host of morepaul.com dropped me and my site after a single, spam-collected e-mail.

Thanks to Google Cache, I was able to recover most of my posts (though little of the uploaded content), and have switched to Blogger for now.

In many ways, the change is very good: like when I moved from LiveJournal to More Paul, I think it's time for a change in tone. In this iteration, I hope to include:

  • Thoughts on Code: You wouldn't know it from my previous blog, but I think about code a lot. I'll hopefully include code snippets of puzzles I'm solving, projects I'm working on, etc.

  • Thoughts on books, especially relating to code: Like any good programmer, I'm always reading. I hope to keep a list as to what's working for me and what isn't.


This hopefully won't come at the expense of my previous tone, which was "here's what I'm doing with my life right now." But like the author of this great A List Apart essay, no technical detail, no matter how arcane, will be as uninteresting as my previous blog if I care about it enough.

Besides: it's a blog. It's an invitation for vanity. So if I tell you how to manually curry functions in Scheme with macros, I will, dammit!

More later, let this just be a stopgap. In the meantime, welcome!

(for those who were wondering about the title, it's a reference to The Little Schemer, a most wonderful book for anyone (of any level) to learn to think about problems recursively. Near the end you explore some pretty hip theoretical topics, but you wouldn't know it if I hadn't told you ^_^).