Computational Projects

I realised math.floor(math.sqrt(n)) was based on floating-point arithmetic, so decided to restore the purity of my integer Python by writing a replacement. It finds the greatest square ≤ the input via a log(n)²-time binary representation–type method.

E.g. for 130, it tests 1², 2², 4², 8², 16² (too big), backtracks to 8², then tests (8+1)², (8+2)², (8+4)² (too big), then (8+2+1)², (8+2+2)² (too big), so the answer is 11.

My initial attempt (test 1², 2², 3² …) was too slow on the inputs I needed it for, so I considered how to traverse numbers faster. This was fast enough even for massive numbers where the floating-point method started giving slightly inaccurate answers.

def rf(n):          # returns floor(sqrt(n))
    i = 1; x = 0                # x stores our progress so far; i is tested, then incremented.
    while True:                 # (infinite loop)
        if (x + i) ** 2 <= n:   # if (x + i)^2 <= n,
            i <<= 1             # double i and try again.
        elif i != 1:            # else we've gone too far! as long as i is not 1,
            x += i >> 1; i = 1  # add the previous i to x and set it back to 1.
        else:                   # but if i is 1,
            return x            # then we are done.

(I added the comments in later to try to motivate some people on The Student Room to take a crack at using Python for CATAM :P.)

Update (Dec 2017): I realised while lying in bed last month that the binary searches I’d been implementing to practice for algorithms tests can solve this problem once an upper bound is determined by the first sequence of 2ⁿ tests. Both parts have complexity log(n), so the whole algorithm does, too. Upon looking it up, I learned that this is called an exponential search.

def euclid(a, b): return a if b == 0 else euclid(b, a%b)

A natural use of recursion.

The best part of my entire coursework coding-wise was the ker algorithm to find the kernel of a matrix over F₂ via Gaussian elimination. I had a bit of spare time to do it creatively, so took the opportunity to learn 2 new bits of Python: bitwise operations and generators.

Algebra in F₂ (integers modulo 2) resembles basic logic, so bitwise operations seemed like a natural fit for the Gaussian elimination part. So, I converted the input matrix into a list (i.e. column vector) of integers whose binary representations were the rows. Then:

• & (bitwise AND) is F₂ multiplication (per element along a row);
• ^ (bitwise OR) is F₂ addition;
• (B[i] >> j) & 1 equals B_ij.

i = 0; j = 0; p = []
# Convert matrix to row-echelon and find pivot variables p.
while i < m and j < n:
    for k in range(i,m):
        if (B[k] >> (n-1-j)) & 1 == 1:
            B[i],B[k] = B[k],B[i]
            for l in range(i+1,m):
                if (B[l] >> (n-1-j)) & 1 == 1:
                    B[l] ^= B[i]
            p.append(j)
            i += 1; j += 1
            break
    else:
        j+=1

My final factorisation algorithm needed to try vectors from the kernel as seeds until one worked. This motivated me to learn about generator functions.

There were f free variables to assign 0 or 1 in some order for each unique kernel vector (some bound variables were determined by back-subbing).

So, I counted with e from 0 to 2ⁿ – 1, read off the binary rep. of e as the pattern of free variables, found the bound variables and yielded the resulting kernel vector at each step.

This one-vector-at-a-time approach allowed me to fully harness the factorisation algorithm (i.e. test all the kernel vectors necessary) without always generating the full kernel. A more common (fast) approach was to settle with picking one kernel vector that seemed like it should work.

# Generate non-trivial kernel vectors.    
y = [0 for t in range(n)]
for e in range(1,min(1 << len(f),201)):
    # Assign 0/1 to free variables of y as per e in binary.
    for i in range(len(f)-1,-1,-1):
        y[f[i]] = e & 1
        e >>= 1
    # Determine pivot variables of y by back-substituting.
    for k in range(len(p)-1,-1,-1):
        l = p[k]
        x = 0; b = B[k]
        for s in range(n-1,l,-1):
            x ^= y[s] & b
            b >>= 1
        y[l] = x
    yield y
return

cff, mostly workmanlike, contains a suspicious trick for flattening a list of lists f:

[y for x in f for y in x]

No idea how it works but there’s a cyclical syntax going on there…

I had to test huge primes, with pauses in between sessions, so wrote it to yield primes from m to n. rf is integer square-root.

def primes(m,n):
    m = max(m,2)
    while True:
        for r in range(2,rf(m)+1):
            if m % r == 0:
                break
        else:
            yield m
        if m > n:
            break
        else:
            m += 1

This computes inverses modulo p (prime) by exponentiating by p–1. Numbers were growing very fast here, so I decided to multiply and take residue alternately. But exponentiation by squaring seemed like an unnecessary optimisation.

def modinv(a,p):
    a %= p
    if a == 0: print('0 not invertible.'); return
    x = 1
    for i in range(p-2):
        x *= a; x %= p
    return x

Clearly inspired by my bitshifty ker algorithm (above).

It’s possible to uniquely identify quotients and remainders in the space of polynomials over a field (here: F^p, the integers mod p), so we can compute aⁿ mod b, where a, b are polynomials.

mpa(r) is a lambda that multiplies polynomial r by a and reduces it mod b. The binary rep. of n is read right-to-left (n & 1, then n >>= 1); at each step, a is squared (becoming a², a⁴, a⁸ …) and, if the current digit of the binary rep. of n is 1, the current value of a is multiplied into the answer x (which is initially 1).

def exp(a,b,n,p):
    mpa = lambda r: div(mp(r,a,p),b,p)[1]
    x = [1]
    while n != 0:
        if n & 1: x = mpa(x)
        a = mpa(a)
        n >>= 1
    return x

Mathematically, it looks perfectly sensible to say you need to take the HCF of a degree 7 polynomial with a polynomial of the form X^pr – X, where p is a 2-digit prime and r ≤ 7. But you first consider if your algorithm can handle those inputs. You don’t chuck them in and then get confused when your computer starts generating a list of 80798284478113 zeroes that uses up all available RAM and almost requires a restart. First step of Euclid: replace X^pr – X with itself mod the other poly. So why not start by generating X^pr mod the other poly via exp? *taps nose*