Factorizing a number can be done by several methods. The first is a naive way, which iterates over every natural number between 2 and the square root of n, then dividing by it as long as the number is divisible by it.

def factor(n):
    factors = []
    i = 2
    while i*i <= n:   # Iterate up to the square root of n.
        exp = 0       # Start with exponent 0
        while n % i is 0:  # Increment exponent and divide while divisible.
            exp += 1
            n /= i
        if exp > 0:
            # Note: All prime factors of n that are less than i are
            # already gone, therefore i itself is definitely prime if we're here.
            factors.append((i, exp)) # Append the prime factor.
            i += 1
    # If n is not yet 1, then it is now prime.
    if n > 1:
        factors.append((n, 1))
return factors

This requires us to iterate over n^0.5 numbers in the worst case (which is when n is prime). An upper bound for factorizing all numbers 1<x<n is O(n^1.5) (though admittedly that is a loose upper bound since 1<n<10^⁴ contains 12.3% prime numbers, 1<n<10^⁵ has 9.6% and 1<n<10^⁶ has only 7.8%. Still, the ratio changes slowly enough to be considered constant.

A slight performance improvement can be had by pre-calculating the set of prime numbers between i and n (which can be done in O(n*log(n*log(n)*loglog(n)) time with the Sieve of Erastothenes) and then reusing it for every number to be factored. After incurring an initial linear-logarithmic cost, we then do the following:

def factor(n):
    global primes
    factors = []
    for p in primes:
        exp = 0
        while n % p is 0:
            exp += 1
            n /= p
        if exp > 0:
            factors.append((p, exp))
        if p*p > n:
            break
    if n > 1:
        factors.append((n, 1))
    return factors

It is apparent, however, that this produces merely the practically constant improvement of only checking prime numbers instead of all numbers. It's down by a factor 10 if n is small (up to maybe 20 when n is 10^8), but it's still in the O(n^1.5) upper bound.

Can we improve on this?

Yes.

The Sieve of Erastothenes can be slightly modified without changing its O(n*log(n)*loglog(n)) complexity, and the modified sieve will let us factor primes in constant and any number in logarithmic time.

First, the modified sieve:

def sieve(n):
    # We need to fill the first 2 spots with 0's to cleanly ignore the case i < 2.
    table = [0]*2 + [1] * (n-2) # A list filled with 2 0's and n-1 1's.
    for i in xrange(n):
        if table[i-7-] is 1:
            for j in xrange(2*i, n, i): # Iterate over all multiples of i
                table[j] = i # This is where we'd normally mark j as not prime.
    return table

The output of this algorithm is a list that assigns each non-prime number its largest prime factor, and each prime number the number 1.

For example, the output for 0<=n<=10 is:

[0, 0, 1, 1, 2, 1, 3, 1, 2, 3, 5]

When factorizing now, we just need to look into this table to find the next prime factor, instead of finding it by iteration:

def factor(n):
    global sieve
    if sieve[n] > 1:
        # Divide n by the prime factor and recurse.
        return [sieve[n]] + factor(n / sieve[n])
    # n is prime, stop recursion.
    return [n]

(Note that the output of this algorithm is not identical to the previous. Instead of [(2, 5), (5, 5)] for 10^5, it will produce [5, 5, 5, 5, 5, 2, 2, 2, 2, 2].)

What was our worst case earlier is now the best: If a number is prime, the algorithm will immediately halt. The best case from earlier (a small prime to a high power, like 2^20) is now the worst case - but it will still be done in logarithmic time. 2^20 needs 20 recursions and so on. (If the recursion depth bothers you, the algorithm can also be formulated iteratively.)

Since our worst case for a single number is now logarithmic, the whole problem has gone from O(n^1.5) to O(n*log(n)), which is about the complexity we spent building the sieve in the first place. From linear-root to linear-logarithmic is quite an improvement for big numbers.

In our Logic class, we just started on Goedelization of logical formulae - that is, the encoding of the language of first-order (arithmetic) logic into numbers, such that every formula gets a unique natural number. (For an analogy, take the English language and consider every word as a number in base 27, where "a" is 1 and "z" is 26. Something similar works for logic.)

Needless to say, these numbers are unbelievably large for even very short formulae.

For example, the first four natural numbers themselves are (of course) terms in arithmetical logic. However, the language does not actually allow writing these numbers out in their literal form: Rather, a term must be constructed only of the basic operators ("+" and "*") and the basic constants ("0" and "1"), and every term must be enclosed in parentheses.

Our encoding specifies 0xd for 0, 0xe for 1, 0x7 and 0x8 for parentheses, and 0xc for +.

So the hexadecimal code for the term equivalent to 3 would read:

0x77ece8ce8

The decimal form of this number is 32191187944. That's our Goedel number of the natural number 3.

(You don't want to know what the Goedel number for 32191187944 is. I did, but my inefficient recursive function crapped out.)

Here's an algorithm for generating the representation and then converting it to decimal:

def basic_num(n):
    if n < 0:
        raise ValueError()
    elif n <= 1:
        return str(n)
    else:
        return "(" + basic_num(n-1) + "+1)"
 
def goedel(s):
    if type(s) is int:
        s = basic_num(s)
    code = {
        # Logical operators
        "~":0x2, # NOT
        "&":0x3, # AND
        "|":0x4, # OR
        "E":0x5, # EXISTS
        "A":0x6, # ALL
        "(":0x7,
        ")":0x8,
        "=":0x9,
 
        # Arithmetic operators
        "<":0xa, # <=
        "+":0xb, # ADD
        "*":0xc, # MULT
        "0":0xd, # 0
        "1":0xe, # 1
    }
 
    g = 0
    for char in str(s):
        g = g * 16 + code[char]
    return g

Today, I was saved by the knowledge of yet another neat little algorithm. Well, not 'saved' as such; I could have done the job with something less efficient and elegant, but it certainly didn't hurt.

For a programming exercise, I needed to find the GCD of two numbers.

The GCD is the largest factor that both numbers are divisible by. If it is equal to 1, the numbers are called coprime. The amount of numbers in {1 .. n} that are coprime with n is called ?(n) and is important in cryptography.

So I was writing this little application that calculated ?(n), and needed a GGT(a, function.

Now, my first impulse was to do this:

do for i from min(a, to 1 if (a mod i == 0 and b mod i == 0) return i end do

You take the lower of the two numbers, and begin counting down to 1. As soon as you reach a number which both are divisible by (their modulo is 0) you've got your GCD. This works, let me stress that. However, you're doing min(a, modulo calculations there. When comparing 101 with 100, you are literally calculating the remainder of 101 divided by all numbers from 100 to 1 before you finally reach the result.

And a nagging feeling was in my head that I'd seen a more elegant solution for this before.

It took a while and a bit of paper to get it right from my bad memory, but it is this:

if (a < swap a with b do while b > 0 c = a mod b a = b b = c end do return a

First, you make sure that a is the larger of the two numbers (swapping a and b if necessary; after all, GGT(a, is commutatively equivalent to GGT(b,a)).

Then, you calculate the remainder from a divided by b. You shift the numbers left: a becomes b, b becomes the remainder you just calculated. As soon as b becomes 0, you've found your greatest common divisor.

Example: 35 mod 15 is 5. 15 mod 5 is 0. The GGT(35,15) is 5. In the example above (101,100), you will see that it only takes a single modulo calculation - 101 mod 100 - to realize that the result is 1. For this input, it is 100 times as efficient.

Try it. Think of any two numbers, and go through the above algorithm. It's not obvious at first sight, but you'll quickly see that it always works.

Oh, and since I tagged this Java, here the implementation:

Of course, it is trivially easy in PHP, too:

This is a little algorithm snippet that has come in handy for me three times in the past month since I came up with it, so I thought I should share it.

A lot of algorithms deal with matrices, or two-dimensional arrays. Looping through each cell in an array, either to find something or update them all, has been done for ages (and the commonly used loop indices are 'i' and 'j', folks, not 'x' and 'y'. Please.)

Now, what isn't as common, but still occurs quite frequently, is the task of checking all the eight neighbours of a certain cell. This is done in labyrinth-traversing algorithms, searching for word groups in a matrix of letters, and (famously) in Conway's Game of Life, whose rules hinge on the status of each cell's eight neighbours. In the following example, I'll use the latter task, because it's very straightforward: In a boolean[][] array, find the number of neighbours that are 'true'.

Having written (and also corrected) a number of exercises dealing with such algorithms, I've seen it all: The record WTF is held by someone who made eight functions: checkLowerLeft(x,y), checkUpperRight(x,y), etc.

But this is common as well:

To be fair, it's not really easy to normalize. When we're taught in Coding Style 102 ("Loops"), we hear that we should use loops whenever we encounter a sequence of identical commands with only a single changing index. Or a nested loops if we need to traverse a rectangle. But what if we need to traverse a ring - more importantly, a ring that's only 3x3 cells large?

The nested loop is perfectly acceptable of course:

Still eight lines, however. Not that much of an edge over the eight repetitions of the single command.

Enter modulo mathematics and integer division. The fact that we have a fixed, small number of positions to check allows us to do something you should never, ever do when traversing a whole matrix: Loop over a single index and calculate the positions from this.

Because the height increment (-1..1) and the width increment (-1..1) are functionally dependent on a single index that goes from 0 to 8:

And in this way, using only two tiny math expressions, we've unified the two loop indexes into a single one. I'm reasonably sure this also shortens compile and run time.

In other words, this:

So I implemented Conway's Game of Life in Java. On a torus, a body whose edges are mapped onto each other to give it an infinite area easily navigated using modulo arithmetics.

It took 140 lines, but that's not spectacular - being the solution to an excercise I wrote on paper earlier today, it's written more for elegance and readability than number of lines, and has plenty of comments.

The source files are here, and attached to this post in a tarball. I usually make a point of not distributing compiled binaries, and given this is Java, you wouldn't be able to use it without a VM (and the compiler that usually comes with it). So you'll have to compile these yourself. But it's pretty straightforward, only two classes.

Have fun!

The exams are over. I got 2.7 on Data Communication, 2.7 on Statistics and 2.0 on Computer Architecture (which would be C+ and B, respectively). Now I have some time to myself. To celebrate that, I fiddled a bit with PHP and implemented a text search algorithm I'd learned. It's called the Boyer Moore Sunday algorithm, and works as follows:

• Begin with the offset 0.
• Test the substring of haystack that begins on offset:
• Begin with needle-offset 0.
• Compare the corresponding characters of needle and the substring of haystack.
• Continue with 4. until the string is found or there is a mismatch.
• If there is a mismatch, examine the character in haystack immediately after the substring.
• Take the last position of this character in needle, or -1 if it does not occur.
• Shift the offset in haystack to the right by the length of needle, and to the left by the position you just calculated.
• Continue with 2.

The implementation is here. Note that PHP differs from compiled languages in some syntax points - for example, a String is simply an array of characters and can be treated as such. The return value is typeless, it can be a boolean or an integer.