Recently, I came across the problem to compute expressions in the modulo arithmetic with integer numbers. Here is an example: We like to compute expressions modulo 31. Since 31 is prime the numbers 0 to 30 with addition and multiplication modulo 31 will be a field. Each number a between 1 and 31 will have a unique multiplicative inverse b, such that ab=1 modulo 31. E.g.

>mod(27*23,31)
     1

The inverse of 27 is 23 modulo 31.

How can one find the inverse element? The fastest way uses Euclid’s algorithm. If you compute the greatest common divider of 27 and 31 using Euclid’s algorithm you can derive x and y such that 27x+31y=1. Then x (taken modulo 31) will be the inverse of 27 modulo 31. Here is the computation for these numbers:

\(31=27+4, \quad 27=6 \cdot 4 + 3, \quad 4 = 3 + 1,\)

and by computing backwards

\(1 = 4 – 3 = 7 \cdot 4 – 27 = 7 \cdot 31 – 8 \cdot 27.\)

Note -8=23 modulo 31. We get our result 23.

Of course, we could program that in EMT or Maxima. Maxima contains already a function for the inverse modulo a number.

>&inv_mod(27,31)
                                   23

We also know the small Lemma of Fermat.

\(a^{p-1} = 1 \text{ mod } p.\)

Thus

\(a \cdot a^{p-2} = 1 \text{ mod } p.\)

This can also be used to compute the inverse. To get the result in EMT we need a trick which I am going to explain below. But in Maxima we can compute very large powers, albeit very slowly.

>&mod(27^29,31)
                                   23

For single computations or small numbers, symbolic computation is okay. But assume we need more speed. The computation above does not work in EMT, because the power is too large. Python, by the way, has infinite integers too. So the computation works there too.

>>> print 27**29
 323257909929174534292273980721360271853387
>>> print 27**29 % 31
 23

One of the reasons of this posting is to show how this can be done in EMT and other languages without an infinite integer arithmetic, and moreover much faster. The trick is the fact that we can as well take the modulo in each step of the multiplication. Here is the first simple code.

>mod(27^29,31) // WRONG RESULT!
   12
>function map powmod (a,n,p) ...
$  b=1;
$  loop 1 to n; b=mod(b*a,p); end;
$  return b
$  endfunction
>powmod(27,29,31) // CORRECT RESULT!
   23

However, this is still not the fastest way. For powers, we can use this trick.

\(x^{2y} = \left(x^y\right)^2, \quad x^{2y+1} = \left(x^y\right)^2 \cdot x.\)

This leads to the following code. By the way, the function matrixpower() of EMT uses the same technique. To understand the code, note that it is not efficient to use a mapping function recursively. Thus we call it from another function.

>function powmodrek (a,n,p) ...
$  if n==0 then return 1
$  elseif n==1 then return a
$  elseif mod(n,2)==0 then
$     b=powmodrek(a,n/2,p); return mod(b*b,p);
$  else
$     b=powmodrek(a,(n-1)/2,p); return mod(mod(b*b,p)*a,p);
$  endfunction
>function map powmod (a:number,n:nonnegative integer,p:index) ...
$  return powmodrek(a,n,p);
$  endfunction
>powmod(27,29,31)
   23
>powmod(27,1:30,31)
   [27,  16,  29,  8,  30,  4,  15,  2,  23,  1,  27,  16,  29,  8,  30,
   4,  15,  2,  23,  1,  27,  16,  29,  8,  30,  4,  15,  2,  23,  1]

The last computation is not effective, of course. If we need to compute all powers we should use a simple loop instead. But I included that example to show that powmod() can map to vector input.

As you see, the powers of 27 do not cycle through all the number 1 to 30. For other base numbers, they do.

>powmod(3,1:30,31)
 [3,  9,  27,  19,  26,  16,  17,  20,  29,  25,  13,  8,  24,  10,  30,
 28,  22,  4,  12,  5,  15,  14,  11,  2,  6,  18,  23,  7,  21,  1]

The second reason for this posting is to demonstrate the solution to another problem. We want to partition the numbers 1 to 30 into triplets with a+b+c=0 modulo 31. This is surprisingly easy. We have for all q different from 1

\(1+q+q^2 = \dfrac{q^3-1}{q-1}\)

Thus by the Lemma of Fermat

\(3^0 + 3^{10} + 3^{20} = \dfrac{3^{30}-1}{3^{10}-1} = 0 \text{ mod }31.\)

Our triples can be found easily from this. Just multiply the equation with powers of 3.

>for k=0 to 9;  ...
>  powmod(3,k,31)+"+"+powmod(3,10+k,31)+"+"+powmod(3,20+k,31)+"=0",  ...
>end;
 1+25+5=0
 3+13+15=0
 9+8+14=0
 27+24+11=0
 19+10+2=0
 26+30+6=0
 16+28+18=0
 17+22+23=0
 20+4+7=0
 29+12+21=0

It is also possible to find such a partition with any desired property using linear programming and LPSOLVE in EMT. This should be the topic of another posting, however.