Recently I’ve taken some breaks from coding for profit to code for glory. Google Code Jam Online Round 1 began on Friday with sub-round A. To compete in Online Round 1 contestants had to successfully solve a “small” and “large” input set in the Qualification Round on July 16th. The top 840 contestants from Online Round 1A will compete in Online Round 2 on August 2nd, along with the top 840 finishers from Online Round 1B and 1C. I managed to place in the top 840 and advance to Round 2, but this post is about the problem that stumped me.

The last problem, “Numbers”, was pretty tricky with only about 4% of contestants (93 out of 2394) solving it. The question is (deceptively) simple:

In this problem, you have to find the last three digits before the decimal point for the number (3 + √5)ⁿ.

For example, when n = 5, (3 + √5)⁵ = 3935.73982… The answer is 935.

For n = 2, (3 + √5)² = 27.4164079… The answer is 027.

The possible values of n for the “large input” case are

2 <= n <= 2000000000

meaning you’ll need more than high precision floating point math to solve the problem.¹

I wasn’t able to solve it myself, but I enjoyed making sense of the top solutions after the contest and will explain some of them here. All solutions submitted during the contest, including those I reference here, can be downloaded from the scoreboard. If you didn’t compete in Round 1A it might be fun to stop at this point and try to solve “Numbers” yourself.

A key observation

All the solutions I’ve deciphered come down to one important fact: it suffices to calculate

(*) (3 + sqrt(5))^n + (3 - sqrt(5))^n - 1 mod 1000

Before making sense of (*) let’s establish some notation:

rp = 3 + sqrt(5)
rm = 3 - sqrt(5)

So, with this notation, the problem is to find

floor(rp^n) % 1000

and (*) can be rewritten as

(*) rp^n + rm^n - 1 mod 1000

To see (*), note that

(**) rp^n + rm^n is integral for all n

by expanding both terms via the binomial theorem and canceling the odd powers of sqrt(5), and so (*) follows from (**) and the fact that

0 < (3 - sqrt(5))^n < 1 for all n

OK, so how do we calculate rp^n + rm^n? There are at least three types of solutions in the top 10, but they’re mostly related to linear recurrences.

My favorite solution

Before talking about recurrences though, I’ll explain my favorite solution, by the 7th place finisher Reid (who might be Reid Barton, a hell of a problem solver apparently). It’s very light on theory, and easy to code. We’ll get it as a special case of a general formula.

If

p = a + sqrt(b)
m = a - sqrt(b)
R(n) = p^n + m^n

then for e = 0 or 1 we get

R(n) R(n+e) = p^{2n+e} + m^{2n+e} + (p^e + m^e)    (pm)^n
            = R(2n+e)             + 2 (1 + e(a-1)) (a^2 - b)^n

If we think of R(n) as an exponentiation, this formula gives a generalization of the recursion

x^{2n+e} = x^n x^{n+e}

used for efficient exponentiation which comes up later. In our case

a = 3
b = 5

and so

R(2n+e) = R(n) R(n+e) - 2 (1 + 2e) 4^n

Reid’s Haskell implementation is already quite readable, but a simple Python implementation is

@memoize
def R(n):
    if n == 0:
        return 2
    elif n == 1:
        return 6
    else:
        q,e = divmod(n,2)
        return (R(q) * R(q+e) - pow(4,q,1000) * 2 * (1 + 2*e)) % 1000

def solve():
    N = input()
    return '%03i' % ((R(N) - 1) % 1000)

Linear recurrences

Recall that given a linear recurrence

a_{n+2} = c_1 a_{n-1} + c_0 a_n

you can express a_n as a linear combination of roots of the
recurrence’s characteristic polynomial

c(x) = x^2 - (c_1 x + c_0)

if all the roots are distinct. This fact usually comes up in relation to the Fibonacci numbers in an algorithms or linear algebra class.

Since

(x - rp)(x - rm) = x^2 - 6x + 4

we see that the recurrence

a_{n+2} = 6 a_{n+1} - 4 a_n

has a solution

a_n = p rp^n + m rm^n

for some numbers p and m. Since we really want to calculate

rp^n + rm^n

we take p = m = 1 and get base cases

a_0 = p + m = 2
a_1 = p rp + m rm = 6

So, we’re in business if we can calculate a_n for n as large as 2 billion.

A naive solution

A naive looping implementation of the recurrence in C calculates the 2 billionth term in about 45 seconds on my 1.66 GHz Core Duo. So, it’s easy to calculate any particular term in the sequence in the allotted 8 minutes, but calculating all 100 terms in the “large input” set separately is about 10 times to slow. However, there are some obvious workarounds:

1. store all 2 billion terms on disk (about 4GB) and have an intelligent way to index into them
2. sort the input and then iterate through the input in parallel with calculating the sequence, storing only the terms you need (or more generally, have some way to quickly check if the current term in the sequence is one of the input terms: it’s much too slow to simply iterate over the whole list of 100 inputs each time you calculate a new sequence term)

A decent C programmer (not me) implementing (2.) would sort the input in C, but here’s a C version where I’ve used Python to pre-sort the input

#include <iostream>

int length = 100;
/* {input, case #, place to store output} */
int I[100][3] = {{103,14,0},{104,70,0},{105,85,0},{203,41,0},{204,78,0},{205,40,0},{1023,55,0},{32769,74,0},{75225234,23,0},{209877672,26,0},{276560520,27,0},{308200077,12,0},{330356541,45,0},{365921634,60,0},{368818964,36,0},{410638565,73,0},{444322222,4,0},{526886908,100,0},{546776723,39,0},{587767441,69,0},{590464207,77,0},{736875403,91,0},{800561938,13,0},{861924427,61,0},{903179180,88,0},{910062006,1,0},{936105644,59,0},{955677912,62,0},{1000875506,18,0},{1011977351,33,0},{1016518317,8,0},{1051461600,64,0},{1064277706,29,0},{1065829044,49,0},{1073741823,65,0},{1073741824,46,0},{1095316023,76,0},{1109957771,31,0},{1123049112,48,0},{1155744626,21,0},{1177572005,11,0},{1184485437,20,0},{1191967198,34,0},{1193741318,7,0},{1225400731,30,0},{1230910728,63,0},{1260829455,52,0},{1284145748,51,0},{1317056443,50,0},{1329765563,53,0},{1365380858,35,0},{1376932921,16,0},{1449613541,2,0},{1477014138,44,0},{1480436753,43,0},{1514881215,66,0},{1526692805,96,0},{1535334186,17,0},{1540548799,38,0},{1543660610,93,0},{1544963801,58,0},{1562234411,25,0},{1570752061,9,0},{1598846451,75,0},{1605814285,5,0},{1611433491,37,0},{1661296617,54,0},{1731301533,80,0},{1736284618,90,0},{1737629019,84,0},{1746961562,98,0},{1791276241,92,0},{1795024971,56,0},{1809415844,22,0},{1813833154,89,0},{1865363233,42,0},{1880802833,94,0},{1894509466,32,0},{1918517182,71,0},{1995422805,24,0},{1999999981,83,0},{1999999982,86,0},{1999999983,81,0},{1999999984,57,0},{1999999985,3,0},{1999999986,67,0},{1999999987,97,0},{1999999988,28,0},{1999999989,15,0},{1999999990,68,0},{1999999991,19,0},{1999999992,6,0},{1999999993,79,0},{1999999994,47,0},{1999999995,10,0},{1999999996,82,0},{1999999997,95,0},{1999999998,99,0},{1999999999,72,0},{2000000000,87,0}};

int main(int argc, char ** argv) {
  int a = 2, b = 6, c = 0;
  for (int i = 1; i <= 2000000000; ++i) {
    int b_ = (6*b - 4*a + 4000) % 1000;
    a = b; /* a is the ith term in the sequence */
    b = b_;
    while (I[c][0] == i && c < length)
      I[c++][2] = (a + 999) % 1000;
  }
  for (c = 0; c < length; ++c)
    printf(”Case #%d: %03d\n”, I[c][1], I[c][2]);
}

and the unix sort to post-sort the output

$ time ./a.out | sort -t # -k2,2g > C-large.out
./a.out  49.76s user 0.00s system 99% cpu 49.784 total
sort -t # -k2,2g > C-large.out  0.00s user 0.00s system 0% cpu 49.784 total

So, with a little bookkeeping, we calculate all 100 terms needed in the time it takes to calculate the last of them.

Periodicity

I’ve see two approaches used to speed up the calculation employed by the top 10 finishers and they’re both more sophisticated than my suggestions above. The 5th place finisher, Plagapong, has a simple solution: since we’re modding out by 1000 there are only 1000^2 possible values for (a_{n-1}, a_{n-2}) in

a_n  = 6 a_{n-1} - 4 a_{n-2}

and so the mod 1000 sequence must be periodic with some period at most 1000000. Plagapong wrote a helper program to calculate the period (it’s 100) and where the sequence starts being periodic (at n = 3).

Fast matrix exponentiation

The 1st and 2nd place finishers, Bohua and yuhch123, took a more heavy weight approach. Anyone that’s familiar with the relationship between exponentiation and linear recurrences probably also knows that exponentiation comes from viewing the recurrence as powering a matrix. In this case, the matrix equation is

(a_n a_{n-1}) = m^n (6 2)

for

m = 6 -4
    1  0

(yuhch123 uses (6 2), which corresponds to the base cases (a_1 a_0). Bohua, OTOH, uses (28 6), which corresponds to (a_2 a_1)). Doing exponentiation “efficiently”, namely finding x^n with O(lg n) multiplications, is a common trick, and both Bohua and yuhch123 employ it.

Here’s a Python version of their solutions:

def multmod(m1,m2):
    m3 = [[0,0],[0,0]]
    for r in range(len(m1)):
        for c in range(len(m2[0])):
            m3[r][c] = sum(m1[r][k]*m2[k][c]
                           for k in range(len(m1[0]))) % 1000
    return m3

@memoize
def powmod(n):
    if n == 1:
        return [[6,-4],[1,0]]
    else:
        q,r = divmod(n,2)
        m = multmod(powmod(q),powmod(q+r))
        return m

def solve():
    N = input()
    return ‘%03i’ % ((multmod(powmod(N),[[6],[2]])[1][0] - 1) % 1000)

where memoize is a memoization decorator defined by

def memoize(f):
    cache = {}
    def g(*args):
        if args not in cache:
            cache[args] = f(*args)
        return cache[args]
    return g

Fibonacci numbers

The next type of solution, employed by the 3rd and 4th place finishers neal.wu and newman, uses the Fibonacci numbers. Let

pp = (1 + sqrt(5))/2
pm = (1 - sqrt(5))/2

be the golden ratios and let F_n denote the nth Fibonacci number. Then

pp^2 = 1 + pp = rp/2
pm^2 = 1 + pm = rm/2
pp pm = -1

and so by the closed form expression for F_n

       pp^n - pm^n
F_n =  -----------
         sqrt(5)

we have

          (rp^n + rm^n) - 2^{n+1} (-1)^n
(F_n)^2 = ------------------------------
                      5 2^n

and so

rp^n + rm^n = 5 2^n (F_n)^2 - (-2)^{n+1}

With that you’re mostly done, but you still need to calculate F_n. Both of them use fast matrix exponentiation of the Fibonacci matrix

1 1
1 0

There’s still more to learn

I’m still working on the 6th place solution by Ahyangyi, which is another matrix power solution. Ahyangyi’s matrix isn’t a recurrence matrix, but it is similar to the recurrence matrix we saw above, meaning powering it gives linear combinations of rp^n and rm^n. I can understand that Ahyangyi’s solution works using matrix similarity together with the matrix identity

a k*b   a' k*b'   a'' k*b''
      *         =
b   a   b'   a'   b''   a''

(where a, a’, b, b’, and k are arbitrary numbers, and a” and b” depend on a, a’, b, and b’), but I have no idea how to derive it in a plausible way. Can anybody derive it?

---

Footnotes

1. A floating point library won’t solve the “large input” set, but it’s enough for the “small input” set, where

2 <= n <= 30

There are some Python solutions using the decimal library, and
yuhch123, the 2nd place finisher whose solution to the “large input” set is analyzed above, used the unix calculator bc to solve the “small input” set. yuhch123 used the bc script

scale=50
for (i = 2; i <= 30; i ++)
   (3 + sqrt(5))^i
quit

and some C code for IO.

I recently accepted an invitation from MetaEzra–blog on all things related to Cornell University–to participate in a panel discussion of the 08-09 New Student Reading Project selection. Incoming students to Cornell read a selected book over the summer, so what they meet in the fall, deep, intellectual conversations replace heavy-petting and makeout sessions.

This time around, they’re reading Lincoln at Gettysburg, a choice I presume meant to tease out comparisons to the present political moment. The panel consists of twenty-something Cornell graduates involved in national politics. One guy is the national volunteer organizer for Obama; another guy is managing editor for the American Spectator and commentator on Fox News. I rep the federal judiciary and Sonya Labs, ya heard?

Editor-in-Chief Matt Nagowski said something about wanting from me a perspective slightly detached from the prevailing winds of the campaign year. The judiciary is often portrayed as the politically neutral third branch of government. Judges have lifetime appointment. And the application of the law is portrayed as a neutral, scientific process.

Law is often expressed as the accumulated wisdom of civilization; therefore, progress comes slow. This, of course, attempts to justify the process of law. On the other hand, reason provides the foundation for law as an inquiry, therefore reason is expected to prevail in the long run. Often, a new set of facts will show how courts stretched an old legal rule beyond the rationale that justifies it. The new decision will then wipe away the old rule by reasoning from underlying principles.

As you know, faithful reader, sanity has not prevailed in the industry of legal research… but it will, naturally. I’ve read some hand-wringing in the blogosphere concerned with the conservatism of law and lawyers in adopting new tools. Using the established tools lends credibility, authority, and a fall guy. No one ever got fired for using IBM!

However, law is done by judges who stick around forever, doing law the same way. Legal research, on the other hand, is done by young associates and law clerks who cycle in and out quickly. In other words, in law, the judge himself is an entrenched barrier to change. Legal research lacks a comparable wall to climb.

See, in law school, the use of precedent, the understanding of the structure of law, etc, are all learned, but they also jive with some sense of fairness and rational governance that we have in our heads, just walking around. In law school, too, lawyers-to-be learn about how legal research is done, but that process doesn’t cohere with any intuition about how research ought to be done.

The new case of legal research presents a novel set of facts. In the old days, Westlaw and Lexis sat firmly on the notion that research generally was a gigantic pain and a logistical nightmare. Young lawyers today have used Google, so they know what’s possible in terms of the ease with which they can retrieve information. Therefore, we must wipe away the old rule that said Wexis was broken but still acceptable, because that was all we knew. The world has changed; irrationality cannot persist for long.

From Sonya Labs: Can Startups Be a Pittsburgh Thing? (rustvalley.com):

I’m not sure why the Sonya Labs blog attracted the attention of so many of my friends and colleagues, but it turns out that Shimon was only the first of three people who pointed me at that blog yesterday and today. Interesting. Sonya Labs has a higher profile among geeks who know me than they probably thought they did.

Now, if only they could figure out how to turn that brand-awareness in that very small and very odd market into money…

Rodrigo’s piece about Pittsburgh has sparked some healthy debate in the blogosphere. Check out the ongoing conversation:

• Hacker News | Forget Silicon: How to Be Steel Valley (news.ycombinator.com)
• Why you SHOULDN’T start a tech company in Silicon Valley (Tim Marman’s Loosely Coupled)
• Sonya Labs: Can Startups Be a Pittsburgh Thing? (rustvalley.com)
• NY Times: 36 Hours in Pittsburgh (rustvalley.com)

…and the Fame-O-Meter continues to rise.

Of all places, I never would’ve expected to build my startup in Pittsburgh. I moved to the ‘burgh from Chicago when Sonya Labs got a seed-stage investment from AlphaLab. It is not so unfathomable that I’m here, though, it actually makes quite a bit of sense. Even Paul Graham, a Pittsburgh native, who is famous for advising that startups go to Boston or Sillicon Valley says so!

Pittsburgh has the opposite problem: plenty of nerds, but no rich people. The top US Computer Science departments are said to be MIT, Stanford, Berkeley, and Carnegie-Mellon. MIT yielded Route 128. Stanford and Berkeley yielded Silicon Valley. But Carnegie-Mellon? The record skips at that point. Lower down the list, the University of Washington yielded a high-tech community in Seattle, and the University of Texas at Austin yielded one in Austin. But what happened in Pittsburgh? And in Ithaca, home of Cornell, which is also high on the list?

…

Rich people don’t want to live in Pittsburgh or Ithaca. So while there are plenty of hackers who could start startups, there’s no one to invest in them.

Well, maybe Paul doesn’t quite say so, but he does try to explain why there isn’t a startup hub here. It simply makes sense for one of the top CS programs in the country to turn its home into one. His explanation is a bit flawed, though, no rich people in one of the steel capitals of the past century? Maybe there is no young money, but there has got to be money somewhere. Not only that, but Pittsburgh is close to all of that east coast money.

On the other hand, if all the hackers that study here are happy to leave and there is nothing to attract hackers to come, money doesn’t count for much. Pittsburgh seems to have more of that problem than a lack of money.

The puzzle goes a bit deeper. For a web startup made of cheap servers, ramen and young people, the main expenses to get started are rent and food. Pittsburgh then, being a lot cheaper than Boston or the Valley, makes even more sense. Here our AlphaLab money will last us at least six months and possibly even eight, whereas it would only last three or four in the other places.

Of course, this puzzle is a lot more complex than solvable with a single answer, but if I was a CMU hacker working on something startup worthy, there are just about enough incentives to stay in Pittsburgh; and if I was a mere mortal hacker there even seem to be enough incentives to come here.

Mike Madison hits closer to what I think is one of the main problems:

This puts the lie, I think, to the most famous and durable of Pittsburgh stereotypes, that this town and region are noteworthy for their honesty and work ethic. If you have a job, that stereotype certainly seems to fit — but a big part of that job seems to be keeping it intact, and keeping others at arms’ length. Is there a “Not Welcome” sign posted in the region’s employment markets? It sure seems that way.

I think that more than a “Not Welcome” sign, Pittsburgh has no sign at all. At least when it comes to startup founders. Since landing here, I’ve noticed the vibe of reinvention just about everywhere I’ve set foot. The city is and has been trying to renew itself into being a startup hub among other things.

Whereas I could always feel the academic attitude of Boston, the entrepreneurial one of the bay area and the alterno-hipster culture of Portland from afar, it was only after a few weeks here that I felt anything about Pittsburgh. I’ve been thinking about startups for the last three years and never once had Pittsburgh crossed my mind as a potential place, but places like Austin, Seattle, Chicago (before I moved there) and even Portland had. There isn’t enough noise!

How to go about hanging the right sign at the door is also tricky business. The first part would be a bit of a media campaign in the adequate parts of the blogosphere. I mean nerds and hackers blogging to nerds and hackers about nerdiness, hackery, and Pittsburgh. A bit of a catch-22, if they aren’t coming or staying here in the first place, though. Maybe we can help. We’re nerdy looking and we mumble “linux, “open source” and “python” enough to pass for hackers.

Programs like AlphaLab are also a good start. If the program is successful in a couple of years there will be a solid network of alumni founders in the city which will make their noise and attract more people in turn. Even more important than helping founders start is helping people who are already here continue, though. Ahem.

Self interest aside, another important piece to attracting founders is getting VCs to invest smaller amounts in more companies. Convincing VCs to invest angel-like amounts and angels to behave like VCs will be part of the next wave of web startups wherever it happens. The trend about the costs of startups are clear by now, so I won’t go into why most startups that are a few months old don’t need several million, but only several hundred thousand to move forward. If you want to give us a few million at the right valuation we won’t object, though.

Finally whatever “Not Welcome” sign needs to be done away with. Duh. There is no room for anti-immigration attitudes. Letting Mike make my point:

Moreover, there are sizable communities in the Pittsburgh region that see potential increases in immigration rates as undesirable — either because immigration of lower-skilled workers threatens existing blue-collar employment and depresses wages, or because in-bound higher-skilled workers compete for positions with people who already live here, or both. Somewhere in Pittsburgh, someone is asking why Sycor wants to raise the H1-B visa cap rather than hire skilled people who already live in Pittsburgh.

…

The problem, in other words, is that immigration is perceived by many as a threat to the pie that we already have, rather than as part of a process of growing the pie.

That one is much more difficult so I’ll let people like Mike work on that problem.

Nathan Collins, CEO

From LexisNexis Versus Westlaw - Update (Legal Research Plus)

The law school library respondents preferred Westlaw at a rate of nearly four-to-one, and seventy-eight percent of the federal court/government respondents preferred Westlaw. And, sixty percent of the law firm respondents preferred Westlaw over LexisNexis. Of the library communities of respondents, no single group preferred LexisNexis over Westlaw.

Paul Lomio and Erika Wayne over at Stanford Law School have completed a truly awe-inspiring work on differences between Lexis and Westlaw.

In reading the survey comments in the article, I’m struck by how age often seems to account for the difference between Lexis and Westlaw. One librarian writes that she simply “trust[s] Shepard’s more” because she’s “a dinosaur.” Others write about how older attorneys to all citation checking as “Shepardization” which younger attorneys have to understand means either Shepardize or KeyCite a case.

Some of the librarians trace the differences between Shepard’s and KeyCite back to how they work in paper form. I’m impressed by this knowledge because I don’t fully understand it and consider it to be somewhat arcane. In law school, I remember getting trained on the books, but of course when I got electronic access I never looked back. I imagine that distinguishing Westlaw from Lexis on the basis of how West book digests work is a legacy of aging librarians, God bless ‘em.

If it’s as I say, that differences in perception of Lexis versus Westlaw are on account of age, I think it’s because, if at one time they had substantial dissimilarities, they have evolved to the point of being nearly identical. Most recent graduates that I know consider them to be completely interchangeable, and will simply use whichever service their firm sticks in front of their face. In law school, if you told me which service was offering better sandwiches at their lunchtime talks and I could tell you which service I couldn’t live without.

From Part II: A Review of PreCYdent.com: A “Google-Styled” Search Engine for Lawyers and Laymen (Res Ipsa Blog)

While users do not have to register to use PreCYdent, there are a number of unique features that are available only to registered users. For instance, registered users can rate cases on a five-star system, mark cases as relevant, add keyword “tags” to cases, as well as upload opinions, and statutes. There are also discussion pages where registered users can post comments or questions about cases or statutes. Additionally, PreCYdent has a professional networking feature that allows registered users to connect with other attorneys.

Benson at Res Ipsa Blog has a nice writeup about PreCYdent. While I’ve heard the above before when distinguishing PreCYdent, I have yet to read anyone tell me why those features are useful.

For instance, Netflix allows users to rate individual movies, so that Netflix can later guess how much I will like a particular movie and make me recommendations based on that. The main reason for this is because films in the abstract lack any networking structure. Metacritic, on the one hand, takes advantage of the existing network of critics who rate films. Netflix, on the other hand, makes use of its user base to rate films and come up with recommendations.

But case law already has a network structure ripe for exploitation. The main question on my mind, then, is why would I want some random user’s rating of a case when judges, arbiters of what the law is, have already “rated” a case by citing or not citing it in an opinion?

Movie preferences are inherently subjective. In other words, I may care more about whether Joe in Idaho likes WALL-E than if Ebert likes WALL-E because Joe and I share more similar taste in movies. However, as an attorney, I don’t care what Joe in Idaho thinks about New York Times v. Sullivan, I care about what a judge thinks. And judges have already indicated their thoughts through the network citation structure of law.

In this situation, there exist two types of networks that can be exploited: social networks and the network structure of law. However, the latter is essentially a social network where judges like Learned Hand and Richard Posner are friends who share citations and post viral videos on each other’s wall. The former is a social network with people like me. Would you prefer my take on Law & Economics or Judge Posner’s? Just sayin’.

From Most Corporate Blogs Are Unimaginative Failures (blogs.wsj.com)

Many businesses have launched corporate blogs in an effort to better communicate with customers and capture a little Web-2.0 mojo. But Huffington Post they ain’t: Not only are these corporate blogs boring as paint, but the businesses behind admit they don’t have much value.

Unfortunately, we aren’t rich enough yet to get talked about in the Wall Street Journal. But one day, god willing, the Wall Street Journal will diss our blog, too.

One of Sonya Labs’ fans asked me whether we had a typo on our About page. I told him that the legal research industry has got to get to Web .20 before it can even think about Web 2.0. Some of our competitors have hilarious “Web 2.0″ features that I won’t even get into right now. Jumping from 1990 to 2008 without stopping through the years in between decontextualizes “Web 2.0″ in the worst way. A question as a point of clarification: if, in 1990, you picked as your domain name a normal word with all of the vowels removed, would you get rich?