From Formulas to Algorithms

Before computers, computation went like this: a scientist invented a formula; the user plugged in suitable values and evaluated the resulting variable-free expression. The user typically did not have the expertise to judge the correctness of the formula. Though the compilers of formula books tried to cover all important cases, it was common that a user failed to find one for the situation at hand. In such compilations one can find formulas for Annuity Whose Present Value Is One, Return on Investment, Moment of Inertia (for circular sheet, hollow circular cylinder, and dozens of other shapes), Curved Surface of Cone, Net Present Value, and many others (but maybe not the one you need right now).

Initially, computers reflected the age of the formula: FORTRAN comes from FORmula TRANslator. But FORTRAN was found more useful for writing algorithms in which formulas played a minor role.  Used in the right way, an algorithm can be self-explanatory in a way that the old-style formulas were not. The formulaic paradigm was authoritarian in the sense that the user did not get an explanation of the formula found in the book. Thus the algorithmic paradigm is less authoritarian: access to the source code implies access to an explanation, not of a formula, but of the number you are getting. Yet relics from the age of formulas linger. In this article I discuss two of these and use them to illustrate the power of the algorithmic approach.

Calendrical calculations are, more than those in other areas, encrusted in ancient lore and antiquated methods. Such anachronisms were the motivation for Dershowitz and Reingold to write Calendrical Calculations. They set out to make the subject computational and use an algorithmic language. Yet they remain to some extent under the thrall of the formulaic paradigm. For example, for the standard (“Gregorian”) calendar, they derive from more complex formulas, after considerable algebraic manipulation, that the number of days in months 1,…,m-1 is given by the floor of (367m-362)/12. To prevent all this algebra becoming more complicated, the assumption was made that February always has 30 days, so this still needs to be corrected later.

Let us compare this approach with the algorithmic alternative. To start with, the formula for the number of days in the completed months before a given date was needed as stepping stone to the serial number of the day in the year, starting with January 1 as day 1. In the algorithmic approach we might as well tackle that directly, for example with:

int dayNum(int d, int m, int y) {
  while (m > 1) { m--; d += monthLen(m, y); }
  return d;
}

My claim is that this is easier to explain than the formula. By “explain” I don’t mean a mathematical correctness proof. Many pages of formalism have been expended on justifying such code mathematically, with formally defined Variants and Invariants. The formalism tends to be more forbidding than the algebra of Calendrical Calculations. No, the power of the algorithmic approach is that the thinking behind the algorithm does not need to be formalized. For a programmer the justification does not even become conscious as the code is written. If the reader needs help, it can be given in the form of prose in the vein of:

The sum of d and the number of days in the months 1,…,m-1 is constant. Before entry into the loop this constant is the number that we want the function to compute. The loop maintains the value of the constant while driving m down. As the loop exits we have m equal to one, so that d is equal to the constant.

Of course the algorithm needs a function for the lengths of the months. Our grandparents (as did their grandparents before them) remembered these by:

Thirty days hath September,
April, June, and November;
All the rest have thirty-and-one,
Save February, with twenty-eight days clear,
And twenty-nine in each leap year,

which is almost C:

    int monthLen(int m, int y) {
      switch (m) {
        case 4: case 6: case 9: case 11: return 30;
        case 2: return leap(y) ? 29 : 28;
        default: return 31;
      }
    }

Of course we still need a function for telling leap years. Even after deciding to use C rather than algebra, we have a rearguard action to fight. In textbooks one often finds something like

int leap(int y) {
  return (y%400 == 0) || (y%4 == 0 && y%100 != 0);
}

Whether it is helpful to express it this way depends on how one remembers the rule. A common way to remember it is as a rule with exceptions, which in turn have exceptions. To be specific: a year is normally not a leap year, except when its number is divisible by 4, but it’s normal again when divisible by 100, except when divisible by 400. The following is a direct translation of this way of remembering the Gregorian calendar:

   int leap(int year) {
     if (year%400 == 0) return 1;
     if (year%100 == 0) return 0;
     if (year%4 == 0) return 1;
     return 0;
   }

This completes what I want to say about an algorithmic alternative to that particular formula of Calendrical Calculations. That formula at least has an easily accessible justification. In my next example this is not the case, though the temptation to use a formula is more understandable. Googling “Zeller’s Congruence” gives evidence of the use in programs of a formula as a quick and easy way to determine the day of the week for a given date in <day,month,year> format.

In the 19th century Christian Zeller [G] invented such a formula to compute the day of the week with y the year, m the month, and d the day of the month:

    ((d+(((m+1)*26)/10)+k+(k/4)+(j/4))-(2*j))%7.

Here k = y%100 and j = y/100; all divisions are integer divisions. The value of the formula represents Saturday as 0, and the other days of the week accordingly. (Actually, the month and year may need to be tweaked a bit; see the code further on.) As the formula maps an infinity of dates to a finite domain of seven values, it can be thought of as establishing an equivalence relation among dates. It’s not clear where “congruence” comes from, but that is what it seems to be called.

Zeller’s formula represents the formulaic paradigm: it solves the problem with a minimum of computation, it is not easy to prove correct, and, I guess, it required considerable ingenuity to invent. In the formulaic paradigm this effort is locked up in a black box to be preserved as intellectual capital for future use. Its use requires mindless and obedient substitution of values followed by the prescribed arithmetical operations.

I propose to use algorithms in two ways, minor and major. The minor one is to prove the correctness of Zeller’s formula; the major one is to do away with the formula and substitute an algorithm for it. This substitution has the advantage that the algorithm is more intelligible than the formula. It will turn out to have the disadvantage that more computation is needed to obtain the answer. But the difference in computational requirement is not noticeable when a computer is used to execute the algorithm.

How can I hope to prove the correctness of Zeller’s formula when I despair of even understanding it? It turns out to be easy to show that only a finite number of dates need to considered and that it is, for contemporary computers, a negligible computational task to evaluate the formula for each of these dates.

To see that the checking of the formula can be reduced to a finite number of cases, it useful to view the rule for leap years in the Gregorian calendar as a system of nested cycles rather than nested exceptions:

  • The innermost cycle consists of three normal years followed by a leap year; let’s call this a quad. The length of the quad is 4*365+1, which is 1461 days.
  • Wrapped around the quads is the next cycle: the century. It consists of 24 normal quads followed by one that is exceptional in being one day shorter. Let’s call this cycle a cent. The length of the cent is 25*|quad|-1, which is 36524 days.
  • The outermost cycle consists of three normal cents followed by an exceptional one, which is one day longer. I call this cycle a clav, in honour of Clavius [G], the astronomer of Pope Gregory XIII. Clavius was, I suspect, the brains behind the Gregorian calendar. The length of the clav is 4*|cent|+1, which is 146097 days.

Let me note, if I may indulge in a short digression, that this gives 146097/400 = 365.2425 days as the length of a year averaged over a clav, and, by the cyclical nature of the calendar, this is also the average over all eternity. As the average length of the year is actually 365.2422 days, we see that Clavius, with his simple rule, managed to create a calendar that is, on average, only off by 0.0003 days per year: it takes around three thousand years to be off by a single day.

After this digression I am ready to state

Theorem 1: The formula of Zeller gives the day of week for every day of the Gregorian calendar.

Proof:
We first reduce the theorem, which makes an assertion for an infinite number of dates, to a finite number of special cases. This can be done by observing that the length of a clav is divisible by 7. Hence two dates fall on the same day of the week if they differ only in the year and if this is in such a way that the difference in year numbers is divisible by 400 (you can call this Van Emden’s Congruence, if you must). Thus the correctness of Zeller’s Congruence only needs to be checked for the 146097 days of a single clav.

In mathematics it is a well-known stratagem to make a proof easier by splitting into special cases, like n<0, n=0, and n>0. But mathematicians take a dim view of, say, a dozen special cases. When Appel and Haken [G] split the Four-Color Theorem into 1936 special cases (each to be checked with a computer program), there were lots of people who objected. You can imagine what these people will think of puny Theorem 1 and the need for 146097 cases to be checked by a computer.

After noting this objection, the proof: as the program below terminates normally with as output a single line containing “146097”, the Theorem is proved.

#include<iostream>

int Zeller(int day, int month, int year) {
  int d = day, m = month, y = year;
  assert(!(m<3 && year == 0));
  if (m < 3) { m += 12; y -= 1; }
  int k = y % 100, j = y/100;
// Zeller's formula: h = 0 -> Saturday etc
  int h =
    ((d+(((m+1)*26)/10)+k+(k/4)+(j/4))-(2*j))%7;
/* From Wikipedia article "Zeller's Congruence"
   observed September 15, 2010.
   Convert from h to 1 ~ Monday, ..., 7 ~ Sunday
*/
  return (h+5)%7 + 1;
}
int proof() {
/* Returns the number of days for which the Zeller formula
   gives the true value as obtained by counting backwards
   from a known day. Fails at first case, if any, in which
   the formula does not give the true value.
*/
  int year0 = 1, year1 = 400;
  int dayOfWk = 7; // Dec 31, 400 was a Sunday
  int count = 0;
  for(int year=year1; year>=year0; year--)
    for(int month=12; month>=1; month--)
      for(int day = monthLen(month, year);
          day >= 1; day--
         ) { assert(Zeller(day,month,year) == dayOfWk);
             dayOfWk--; if (dayOfWk < 1) dayOfWk += 7;
             count++;
           }
  return count;
}
int main() {
  printf("%d\n", proof());
  return 0;
}

My reason for preferring the algorithmic to the formulaic paradigm is that I find writing the proof a straightforward exercise, whereas understanding Zeller’s formula seems a daunting task to me. I trust that my readers will not find my code unduly difficult to understand and will find it less mysterious than Zeller’s formula. Alan Robinson has called for the need of  “intelligibility engineering”; I hope he accepts this as a contribution.

So far the minor application of the algorithmic paradigm in matters calendrical. Minor, because the day of the week for distant dates belongs in trivia-land, like the storming of the Bastille happening on a Tuesday. Perhaps a bit more significant is that Zeller’s formula gives Sunday for December 7, 1941. The formula may also be used to assure someone born on December 13, 1973 that this day was not a Friday, so that subsequent mishaps have to be attributed to other causes.

My major example of demonstrating the algorithmic alternative to formulas is to replace Zeller’s formula by an algorithm. The reason for doing so is that the algorithm, although occupying many more times the number of lines of the formula, is intelligible: it can be read straight through compared to the need to puzzle over the single line of the formula.

What is an algorithmic alternative to a single-line formula for determining the day of the week? A promising starting point is the essence of a calendar: a system to name the days. In an algorithmic approach, a calendar uses the fact that the days form a sequence. This suggests to use integers for the names. In this way the calender does not only name the days, but does also count them. And counts not only the days themselves, but the number of days between any two dates.

I follow this approach. The result is a building block that can be used for many calendrical computations. The algorithm for numbering the days will be more intelligible than Zeller’s formula. The algorithm will also more widely applicable, while the day of week is just one of many uses of converting days to integers. Thus the algorithmic alternative to Zeller’s formula will take the form of a function that converts a date in <day,month,year> format to an integer. As a side effect, knowing the day of the week of just one date (like, say, today) implies the equivalent of Zeller’s formula via the remainder on division by seven.

I follow Dershowitz and Reingold in calling the integer name of a day rata diem and I abbreviate it to rd. These authors assign the rd of 1 to January 1 of year 1. Perhaps the idea is that the symmetry that the integers have with respect to 0 can be used in the rd naming scheme. This would be a promising approach, were it not for the fact that the Gregorian calendar, with its positive and negative years, is not symmetric: for positive years time moves away from the origin; for negative years it moves in the same absolute direction, which is towards the origin. As computer integers have a limited range anyway, a computerized calendar has to have an earliest date as well as a latest date. Accordingly, I use unsigned integers for the rata diem. As 32 bits is a commonly available size of unsigned integers and gives a generous number of years into past and future, I choose this size. As there is an even number of integers, there are two midpoints: 7fffffff and 80000000 in hexadecimal. I choose the later one of these as the rd of January 1 of year 1. I call this date the epoch, adopting another bit of calendrical terminology from [DR97].

Here is the first bit of the program:

typedef unsigned uint;
const uint epoch = 0x80000000; // Rata Diem of January 1, 1 C.E.

const int yLen = 365;               // length of normal year
const int quadLen = 4*yLen + 1;     // length of normal quad
const int centLen = 25*quadLen - 1; // length of normal century
const int clavLen = 4*centLen + 1;  // length of clav

Notice that I have switched to the C++ compiler. Not because I want to do any object-oriented programming, but because of the minor civilities that this improved form of C affords. Remember that I relied on my shaky pencil-on-back-of-the-envelope calculations to come to the conclusion that a clav has 146097 days? A man shouldn’t have to do that in this day ‘n age, but C forbids anything except a literal in the initialization of a constant. (I didn’t actually mind using a pencil. The real reason is that I want to replace magic numbers by their provenance, as documentation).

A promising strategy for converting Rata Diem to <day,month,year> is to split the rd into the days coming from the completed years (if any) in the past and the days coming from the current year. The latter number is computed by our old friend:

int dayNum(int d, int m, int y) {
  while (m > 1) { m--; d += monthLen(m, y); }
  return d;
}

Function dayNum() was written with positive years in mind. For negative years the completed years are greater than the current year and the days in the current year are not counted forward from the beginning, but backwards from the end of the year. Rather than tweak dayNum() to a form that works in negative years as well as positive years, I use the fact that the Gregorian calendar is the same for every clav. So I shift every date, irrespective of whether the year is positive, negative, or in between, to the clav that begins on January 1, 1 C.E. and ends on December 31, 400 C.E. I implement this shift by:

 while (year < 1) {
    year += 400; rd -= clavLen;
  }
  while (year > 400) {
    year -= 400; rd += clavLen;
  }

This takes linear time. Yes, it could be done in constant time. But the rule for modulo a negative number is less clear for me, so I prefer the slower alternative. After all, the goal is intelligibility.

See the Appendix for a complete set of C++ functions.

Zeller’s formula, though unintelligible, has the advantage of being computationally efficient, more so than my algorithmic alternative. In my experience the most extreme example of an unintelligible formula is Matijasevic’s polynomial [G] that is claimed to be a representation of the prime numbers. Not only is it a mystery to all except the few initiated why the polynomial should do so, but it is also computationally inefficient in the extreme. The polynomial is of degree 25 and has 26 variables. It is claimed to represent the prime numbers in the sense that, as the 26 variables range over the non-negative integers, all positive values of the polynomial are prime and vice versa. Contrast this, for example, with Dijkstra’s prime-number generator [DDH72], which is intelligible as well as reasonably efficient.

In this article I have concentrated on examples where algorithms are preferable to formulas. I’m not forgetting that there is a whole world of programming that exists on the strength of the opposite: in functional programming one is forced to abstain from algorithms and to express everything in formulas. Of course, the formulas allowed in functional programming form a less restrictive language than the one to which Zeller was limited. A modern language of formulas is enriched with conditionals and with recursion. These facilities introduce possibilities for explanation that were not open to Zeller.

References

[1] Information on this was found with Google.

[2] Calendrical Calculations by Nachum Dershowitz and Edward Reingold; Cambridge University Press 1997.

[3] Structured Programming by O.-J. Dahl, E.W. Dijkstra, and C.A.R. Hoare; Academic Press 1974.

[4] The ability to tell the day of week of distant dates occurs in several calculating prodigies; see The Great Mental Calculators by Steven Smith, Columbia University Press 1983. Another study is an essay by Oliver Sacks (“The Twins”, in The Man Who Mistook His Wife for a Hat, Picador, 1986) describing twin brothers with this faculty. Sacks reports (page 187) that the men had an IQ of sixty, could not do addition or subtraction with any accuracy, and could not comprehend the meaning of multiplication or division. This rules out Zeller’s formula as their secret weapon. In the Postscript on page 200 Sacks quotes one Israel Rosenfeld, apparently consulted for his supposedly superior knowledge of matters numerical:

Their ability to determine the days of the week within an eighty-thousand-year period suggests a rather simple algorithm. One divides the total number of days between ‘now’ and ‘then’ by seven. …

One wonders whether Mr Rosenfeld has tried out his “rather simple algorithm”.

Appendix

#include<iostream> // for assert and printf

typedef unsigned uint;
typedef struct{int day; int month; int year;} Date;
void printDate(Date d) {
  printf("%u, %u, %d\n", d.day, d.month, d.year);
}
Date mDate(int d, int m, int y) {
  Date date;
  date.day = d; date.month = m; date.year = y;
  return date;
}

const int yLen = 365;               // length of normal year
const int quadLen = 4*yLen + 1;     // length of normal quad
const int centLen = 25*quadLen - 1; // length of normal century
const int clavLen = 4*centLen + 1;  // length of clav
const uint epoch = 0x80000000; // Rata Diem of January 1, 1 C.E.

int leap(int year) {
  if (year%400 == 0) return 1;
  if (year%100 == 0) return 0;
  if (year%4 == 0) return 1;
  return 0;
}
int monthLen(int m, int y) {
  switch (m) {
    case 4: case 6: case 9: case 11: return 30;
    case 2: return leap(y) ? 29 : 28;
    default: return 31;
  }
}
int dayNum(int d, int m, int y) {
  while (m > 1) { m--; d += monthLen(m, y); }
  return d;
}
uint Date2RD(Date argDate) {
  int day = argDate.day, month = argDate.month;
  int year = argDate.year;
// Converts date to Rata Diem.
  int rd = epoch-1; // Rata Diem of December 31, 0
// Normalize date to one in standard clav.
  while (year < 1) {
    year += 400; rd -= clavLen;
  }
  while (year > 400) {
    year -= 400; rd += clavLen;
  }
  assert(1 <= year && year <= 400);
  rd += dayNum(day, month, year);
  year --;
  // year contains number of whole years, if any
  if (year == 0) return rd;
  assert(0 < year);
  rd += (year/100)*centLen; year %= 100;
  rd += (year/4)*quadLen;    year %= 4;
  rd += year*yLen;
  return rd;
}
int dayOfWeek(Date date) {
  uint y = 6+Date2RD(date)%7; // day 0 was a Saturday
  return y>7 ? y-7 : y;
}
int main() {
  printf("%d\n", dayOfWeek(mDate(14,7,1789)));
  return 0;
}

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s


%d bloggers like this: