Clock Arithmetic


Oregon Curriculum Network
Discovering Math with Python

Chapter 4: CLOCK ARITHMETIC

Calling modulo arithmetic "clock arithmetic" is a way of calling attention to the finite cyclic nature of always taking remainders, factoring away the "modulus". What does that mean?

Clock Arithmetic is important to Number Theory, and Number Theory is important to crytography.

Cyptography, in the not so distant past, would not concern an average civilian, whereas today the HTTPS protocol in every mainstream web browser (Mozilla, Chrome, IE, Safari, Opera...) is capable of implementing what we called TLS, or Transport Layer Security.

The web browser client, and server, will handshake, meaning agree on what cryptographic strategy to use: first a public one, to exchange a symmetric key, and then a symmetric one, most likely AES at the time of this writing.

AES is the Advanced Encryption Standard, the winner of a global contest as adjudged by NIST. AES uses "clock arithemetic" internally in order to apply several layers of "mixing it up" by invertible methods. Enciphering is reverisble (the whole point) but only if the key is handy.

In Python, one of our primitive operators is %, for "modulo", which relates to the built-in function divmod( ). Lets see how % and divmod do their jobs:


In [1]:
(5 + 5) % 10 # no remainder


Out[1]:
0

In [2]:
divmod(10, 5) # 5 goes into 10 twice, no remainder


Out[2]:
(2, 0)

In [3]:
19 % 10 # divide by 10, give the remainder


Out[3]:
9

Try doing a number of these yourself. Guess the answer before hitting the Shift-Enter key.

You may find your first years of formal schooling did not emphasize cyclic arithmetic, except in the two most common forms, of telling time and scheduling appointments. We understand both the 24 hour day, and the 365 day year, comprise a modulus such that we're always talking about our degree of offset into each one. Adding hours or days, per a timedelta object, keeps us within the scope of some calendar.

The unit circle of 360 degrees, or 2 pi radians, is also treated in a clock-like fashion. Going around thousands of degrees from some present position, never takes us beyond 360. We're confined to a finite domain. Cyclic phenomena more generally are accommodated by our clock arithmetic, of complex numbers, of e to imaginary powers.

There's a kind of "closure" in this picture, which might be resonating with you by now, as one of those properties of a group. Will we find groups in this Chapter? You bet. Groups also enjoy the symmetry of a circle in that each member is paired with another "180 degrees opposite" although sometimes the inverse might be the same as the self.

Again, we call this "clock arithmetic" because even if we say "20 hours from now", or "30 hours ago" we'll still be somewhere on the circle marked out into 12 intervals, each representing one hour. The yearly calendar is likewise a kind of modulus. No matter how many days we add, we're still somewhere between 0 and 365 days into some year.

We also call it "modulo arithmetic" meaning the computations are always vis-a-vis some fixed modulus. This word "modulus" is our inspiration for the M-numbers coded below.

A Class for M-numbers

Lets build a class, the instances of which will multiply and add per some fixed modulus, meaning we're always factoring out the modulus and keeping the remainder.


In [4]:
class M:  # for "modulo"
    
    modulus = 10  # class level
    
    def __init__(self, val, modulus=None):
        if modulus:
            self.modulus = M.modulus = modulus # resettable
        else:
            self.modulus = M.modulus
        self.val = val % M.modulus
        
    def __add__(self, other):
        if self.modulus != other.modulus:
            raise ValueError
        return M((self.val + other.val) % self.modulus)
    
    def __mul__(self, other):
        if self.modulus != other.modulus:
            raise ValueError
        return M((self.val * other.val) % self.modulus)
    
    def __pow__(self, exp):
        raise NotImplemented
    
    def __repr__(self):
        return "(" + str(self.val) + " mod " + str(self.modulus)+ ")"
    
a = M(8)
b = M(7)
print(a, b)
print("(8 mod 10) * (7 mod 10) = ", a * b)
print("(8 mod 10) + (7 mod 10) = ", a + b)


(8 mod 10) (7 mod 10)
(8 mod 10) * (7 mod 10) =  (6 mod 10)
(8 mod 10) + (7 mod 10) =  (5 mod 10)

OK, everything seems to be working, even though we haven't implemented powering yet. Eventually we'd like to go pow(M(3), -1) to get the inverse of M(3), such that M(3) times its inverse equals the multiplicative identity M(1).

Looking for Groups

But wait, does every M-number, with modulus set to 10, have an inverse? We can check that easily. First, lets make a list of all 10 M-numbers, (0) through (9):


In [5]:
elems = [M(n) for n in range(10)]
elems


Out[5]:
[(0 mod 10),
 (1 mod 10),
 (2 mod 10),
 (3 mod 10),
 (4 mod 10),
 (5 mod 10),
 (6 mod 10),
 (7 mod 10),
 (8 mod 10),
 (9 mod 10)]

Now we can do like a "times table" wherein we pick a single M-number, say M(5), and multiply it by every number in elems...


In [6]:
[M(5) * x for x in elems]


Out[6]:
[(0 mod 10),
 (5 mod 10),
 (0 mod 10),
 (5 mod 10),
 (0 mod 10),
 (5 mod 10),
 (0 mod 10),
 (5 mod 10),
 (0 mod 10),
 (5 mod 10)]

Interesting. In ordinary arithematic, the times table for 5 goes 0, 5, 10, 15, 20, 25, 30... and so on. Factoring out the 10s, leaving only remainders, we get M(0) or M(5) as our two M-numbers. We have no way to reach M(1) and so M(5) has no inverse. We don't have a group yet.

Lets try M(2):


In [7]:
[M(2) * x for x in elems]


Out[7]:
[(0 mod 10),
 (2 mod 10),
 (4 mod 10),
 (6 mod 10),
 (8 mod 10),
 (0 mod 10),
 (2 mod 10),
 (4 mod 10),
 (6 mod 10),
 (8 mod 10)]

Same thing. We cycle around and around, always stopping at the same stations, like a model train going in a circle. We never stop at station M(1). So M(2) has no inverse either. What about M(3)?


In [8]:
[M(3) * x for x in elems]


Out[8]:
[(0 mod 10),
 (3 mod 10),
 (6 mod 10),
 (9 mod 10),
 (2 mod 10),
 (5 mod 10),
 (8 mod 10),
 (1 mod 10),
 (4 mod 10),
 (7 mod 10)]

Aha! Now we're getting somewhere. M(3) * M(7) returns M(1), so these two are inverses of one another. One fact to notice immediately is neither has any factors in common with 10. In fact, both are prime numbers in the ordinary integer sense. M(9) is not prime, but again, 9 has no factors in common with 10. So does M(9) have an inverse? Lets check:


In [9]:
[M(9) * x for x in elems]


Out[9]:
[(0 mod 10),
 (9 mod 10),
 (8 mod 10),
 (7 mod 10),
 (6 mod 10),
 (5 mod 10),
 (4 mod 10),
 (3 mod 10),
 (2 mod 10),
 (1 mod 10)]

Indeed it does. M(9) * M(9) = M(1), meaning M(9) is its own inverse.

When positive integers have no factors in common, aside from 1, we say they're "relatively prime" or that they're "strangers" to one another. 3 and 10 are strangers, as are 9 and 10. Sometimes we write 7 | 10 = 1, meaning their greatest factor in common is 1. On the other hand, 6 | 10 = 2, as 2 divides into both without remainder.

What we're about to discover is the strangers to a modulus comprise a group, i.e. M(n) where n | N = 1 and M.modulus = N. A Cartesian product of such "minions" would demonstrate closure, inverse for everyone, naturally a neutral element, and associativity i.e. makes no difference if you go M(a) M(b) M(c) by starting with either M(a) M(b) or M(b) M(c). Commutativity is not a requirement for grouphood.

Finding Totatives

It'd sure be handy at this point, to have a function, gcd, that returns the greatest common divisor of two numbers. Then we could find all the strangers to 10, that are positive integers less than 10, pretty easily. They would have no common divisor with 10.

We call these strangers the "totatives" of 10, and the number of totatives, is called the "totient" of 10.

Note that we're using list comprehension syntax in this example.

The for loop inside the square brackets makes x be every integer from 0 to modulus - 1, but then most of these get filtered out, because of factors in common with the modulus.


In [10]:
import math
modulus = 10
[x for x in range(modulus) if math.gcd(x, modulus) == 1]


Out[10]:
[1, 3, 7, 9]

In [11]:
modulus = 12
[x for x in range(modulus) if math.gcd(x, modulus) == 1]


Out[11]:
[1, 5, 7, 11]

That was pretty easy. Right away we're getting totatives. We could be using set comprehension syntax instead.

Using Python's "list comprehension" syntax, which allows an if clause for filtering, we were able to get the totatives of 10 and 12 respectively. Both 10 and 12 have a totient of four, meaning each has four totatives.

What's true is that the M-numbers that are totatives of some modulus (say 10), collectively form a group under multiplication.

If we confine ourselves to strangers to the modulus, we'll have an inverse for every element, closure, associativity, and M(1) will be included. That's a quick way to get a group.

We don't have closure for addition though. If we want to include addition, along with its neutral element zero, then we'll need to make our modulus a prime number. Then we will be guaranteed a field, a Galois Field to be more precise, though in German they say it differently.

Before we jump to fields though, lets build a function for computing totatives and use that to build some groups.

A "set comprehension" is just like a list comprehension in terms of syntax, but the curly braces make it for a set, instead of a list.

Sets are unordered and their elements are always unique (no duplicates). Since we know totatives are unique, we might as well practice using a set object.


In [12]:
def totatives(n):
    return {x for x in range(n) if math.gcd(x, n) == 1}

In [13]:
elems = {M(x, 12) for x in totatives(12)}
elems


Out[13]:
{(1 mod 12), (7 mod 12), (11 mod 12), (5 mod 12)}

What we just did there was compose a group of M-number objects, of twelves four strangers. If you're using this Notebook live, here might be a good time to insert some code cells of your own, testing out whether the Cayley Table, i.e. the everything by everything multiplication table, really shows closure, and 1 in every row (proof of inverse element).


In [14]:
M.modulus = 100
elems = {M(x, 100) for x in totatives(100)}  # set comprehension
elems  # print M-number totatives in no special order -- it's a set


Out[14]:
{(69 mod 100),
 (1 mod 100),
 (71 mod 100),
 (73 mod 100),
 (77 mod 100),
 (79 mod 100),
 (17 mod 100),
 (81 mod 100),
 (83 mod 100),
 (7 mod 100),
 (87 mod 100),
 (19 mod 100),
 (89 mod 100),
 (91 mod 100),
 (93 mod 100),
 (13 mod 100),
 (97 mod 100),
 (21 mod 100),
 (99 mod 100),
 (9 mod 100),
 (29 mod 100),
 (31 mod 100),
 (3 mod 100),
 (23 mod 100),
 (27 mod 100),
 (11 mod 100),
 (33 mod 100),
 (37 mod 100),
 (39 mod 100),
 (41 mod 100),
 (43 mod 100),
 (47 mod 100),
 (49 mod 100),
 (51 mod 100),
 (53 mod 100),
 (57 mod 100),
 (59 mod 100),
 (61 mod 100),
 (63 mod 100),
 (67 mod 100)}

Note that set objects cannot be ordered even if we might like them to be as by definition they're not sequences. Does Python have an OrderedSet corresponding to an OrderedDict?

Finding Inverses

If we confine ourselves to a set of strangers, we're safe in assuming there's always an inverse for any given element. That suggests a "brute force" way of finding any element's inverse: just multiply it by every element in the same set of strangers, until their product is 1 (the identity element). The function below accomplishes this:


In [15]:
def inverse(n : M, elems : set) -> M:
    for m in elems:
        if (n * m).val == 1:
            return m

In [16]:
inverse(M(79), elems)


Out[16]:
(19 mod 100)

In [17]:
M(19) * M(79)


Out[17]:
(1 mod 100)

Finding Totients

We'll need this concept in the next chapter on Public Key Cryptography. The totient of N is simply the number of totatives it has. We can use the len( ) to simply count them. That won't always be practical, when N gets huge, but for now this approach makes the concept easy to grasp:


In [18]:
def totient(n):
    """how many totatives have we?"""
    return len(totatives(n))

totient(1_000_000) # using _ helps us parse the number, one million


Out[18]:
400000

The number one million has a totient of 400,000, meaning that many numbers from one to one million are co-prime to it.

Adding a Power Method

We now have the technology (tool set) we need to add the __pow__ method to our M class. A negative power triggers finding an inverse of a, the -1 part, then we switch back to positively powering. We're saying $M(n)^{-e}$ equals $(M(n)^{-1})^{e}$.

This is how we treat exponents normally, e.g. 2 ** -2 equals (2 ** -1)**2 or 1/4.


In [19]:
2 ** -2 == (2 ** -1) ** 2


Out[19]:
True

For example, pow(M(79), -2) means find the inverse of M(79) e.g. M(79)**-1, and then raise the result to the 2nd power.


In [20]:
class M:  # for "modulo"
    
    """
    Version 2, with inverting and powering added
    Version 3, not listed, is in groups.py and uses xgcd (next Chapter)
    """
    
    modulus = 10  # class level
    
    def __init__(self, val, modulus=None):
        if modulus:
            self.modulus = M.modulus = modulus # resettable
        else:
            self.modulus = M.modulus
        self.val = val % M.modulus
        
    def __add__(self, other):
        if self.modulus != other.modulus:
            raise ValueError
        return M(self.val + other.val, self.modulus)
    
    def __mul__(self, other):
        if self.modulus != other.modulus:
            raise ValueError
        return M(self.val * other.val, self.modulus)

    def __invert__(self):
        elems = {M(x) for x in totatives(self.modulus)}
        for m in elems:
            if (self * m).val == 1:
                return m
    
    def __pow__(self, exp):  # pow() and ** both trigger this method
        output = self
        if exp < 0:
            exp = abs(exp)  # make exp positive now
            output = ~self  # -1 means take the inverse
        elif exp == 0:
            output = M(1)   # return identity if exp is 0
        elif exp == 1:
            output = self   # return M if exp is 1
        if exp > 1:
            for _ in range(1, exp):  # use __mul__ (already defined)
                output *= self
        
        return output
    
    def __repr__(self):
        return "(" + str(self.val) + " mod " + str(M.modulus)+ ")"

In [21]:
M.modulus = 100
a = M(79)
a_inv = pow(a, -1)
a_inv


Out[21]:
(19 mod 100)

In [22]:
M(13) * M(13) * M(13)


Out[22]:
(97 mod 100)

In [23]:
M(13)**3       # confirm this is another way of saying it (we're testing)


Out[23]:
(97 mod 100)

In [24]:
M(13)**-1      # give me the inverse of M(13) please


Out[24]:
(77 mod 100)

In [25]:
M(13) * M(77)  # these must be inverses then


Out[25]:
(1 mod 100)

Now that we have our inverse function, we can ask questions such as:

  • What is the inverse of M(5, 7)?
  • What is the inverse of M(7, 5)?
  • What is the inverse of M(97, 100)?

In [26]:
print("What is the inverse of M(5, 7)?", ~M(5,7))
print("What is the inverse of M(7, 5)?", ~M(7,5))  # M(7,5) is same as M(2,5)
print("What is the inverse of M(97, 100)?", ~M(97, 100))


What is the inverse of M(5, 7)? (3 mod 7)
What is the inverse of M(7, 5)? (3 mod 5)
What is the inverse of M(97, 100)? (33 mod 100)

Primality Tests

Now that we have a powering method in place, lets take a look a Fermat's Little Theorem. The logic here is useful, in the sense that we sometimes get tripped up by IF / THEN statements.

IF a number p is prime, THEN M(A, p) ** p == A, where 0 < A < p. If some number n passes this test for a given A, where gcd(A, n) == 1, then that's evidence n might be prime, but not proof. If n in fact turns out to be composite, then a is called a Fermat Liar or Fermat Fooler.

Lets take a look. Is 511 prime? Lets find an A with no factors in common with 511.


In [27]:
from math import gcd
gcd(511, 25)


Out[27]:
1

In [28]:
n = M(25, 511) # n will power modulo 511
n ** 511


Out[28]:
(298 mod 511)

Here we say 25 is a "witness" to the fact that 511 is not prime. If it were, it would pass this Fermat test. How about 513?


In [29]:
n = M(25, 513)
n ** 513


Out[29]:
(1 mod 513)

Nope, that's not passing either. Remember we expect A back out again, of 25 in this case. 25 is a witness that 513 is not prime. How about 523 then?


In [30]:
n = M(19, 523)
n ** 523


Out[30]:
(19 mod 523)

In [31]:
n = M(25, 523)
n ** 523


Out[31]:
(25 mod 523)

In [32]:
n = M(31, 523)
n ** 523


Out[32]:
(31 mod 523)

OK, at last it appears we've struck gold. All of our As are testifying that 523 is a prime. Could all of these be liars?

Lets introduce the set of Carmichael Numbers, precisely those composites which sneak through the Fermat Test for any suitable base 1 < A < n.

561 is the lowest Carmichael number.


In [33]:
liars = [M(n, 561) for n in (19, 25, 31) if gcd(n, 561) == 1]
liars


Out[33]:
[(19 mod 561), (25 mod 561), (31 mod 561)]

In [34]:
[n**561 for n in liars]  # all pass the Fermat primality test


Out[34]:
[(19 mod 561), (25 mod 561), (31 mod 561)]

And yet 561 is a composite number, with prime factors (3, 11, 17).

Another test for primality, called the AKS Test, is more fool-proof. We'll need the coefficients of the polynomial (x - 1) ** p, where p is our candidate prime. Those coefficients may be read of Pascal's Triangle. If all coefficients but the first and last (which are 1) are divisible by p, then p is prime.

Lets try that out, first by writing a Python generator for successive rows of Pascal's Triangle:


In [35]:
def pascal():
    row = [1]
    while True:
        yield row
        row = [(i+j) for i,j in zip(row + [0], [0] + row)]

gen = pascal()
for i in range(10):
    print(next(gen))


[1]
[1, 1]
[1, 2, 1]
[1, 3, 3, 1]
[1, 4, 6, 4, 1]
[1, 5, 10, 10, 5, 1]
[1, 6, 15, 20, 15, 6, 1]
[1, 7, 21, 35, 35, 21, 7, 1]
[1, 8, 28, 56, 70, 56, 28, 8, 1]
[1, 9, 36, 84, 126, 126, 84, 36, 9, 1]

That seems to be working. So Pascal's Triangle could be at the heart of a new iterator for producing successive prime numbers.


In [36]:
def primes():
    yield 2
    p = 1
    gen = pascal()
    next(gen) # [1]
    next(gen) # [1, 1]
    while True:
        p += 2 # 3, 5... check odd rows only, once 2 yielded
        next(gen)     # skip even row
        r = set(next(gen)[1:-1]) # drop 1s on both ends and dedupe
        if sum([coeff%p for coeff in r]) == 0:
            yield p

ps = primes()
print(str([next(ps) for _ in range(100)]))


[2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71, 73, 79, 83, 89, 97, 101, 103, 107, 109, 113, 127, 131, 137, 139, 149, 151, 157, 163, 167, 173, 179, 181, 191, 193, 197, 199, 211, 223, 227, 229, 233, 239, 241, 251, 257, 263, 269, 271, 277, 281, 283, 293, 307, 311, 313, 317, 331, 337, 347, 349, 353, 359, 367, 373, 379, 383, 389, 397, 401, 409, 419, 421, 431, 433, 439, 443, 449, 457, 461, 463, 467, 479, 487, 491, 499, 503, 509, 521, 523, 541]

Homework:

Watch this Youtube on the Chinese Remainder Theorem (CRT). Be prepared to explain its utility. What Youtube do you find most elucidating, regarding the CRT?

So now it looks like we have another fun tool for exploring Group Theory, while learning Python at the same time!

Now lets take what we've learned in this chapter and apply it to an important topic: Public Key Cryptography.

Back to Chapter 3: A First Class
Chapter 5: Public Key Cryptography
Introduction