An introduction to Axiom (7): Programming

(It’s been a while since I posted about Axiom – that’s what exam marking does!)

Axiom has a full and complete, and well-structured programming language. Rather than describe it in detail, I’ll just make a few general remarks, and provide a few examples.

A simple function

Here’s an example of a simple program, to solve an equation numerically by the bisection method. We implement the following algorithm:

INPUT: function f,
       values a,b for which f(a)f(b)< 0,
       a tolerance eps
WHILE f(a)>eps REPEAT
  c = (a+b)/2
  IF sign(f(a))=sign(f(c))
    THEN a = c
    ELSE b = c
RETURN a

For this algorithm at least, the program is very similar in Axiom:

bisect(f,a,b,eps) ==
  local c
  while abs(f(a))>eps repeat
    c:=numeric (a+b)/2
    if sign(f(a))=sign(f(c)) then
      a:=c
    else
      b:=c
  return(a)

Note that for ease of exposition I’ve done no sign checking on the values of f(a) and f(b). This little program can be saved in a file called bisect.input and read into Axiom with the command

)read bisect

And now it can be used:

Note the use of anonymous functions.

Axiom programs are made of blocks, which are groups of statements all with the same indentation, or alternatively delineated with parentheses; the statements being separated with semicolons. The above program could be entered as

bisect2(f,a,b,eps) == (while abs(f(a))>eps repeat (c:=numeric (a+b)/2;if sign(f(a))=sign(f(c)) then a:=c else b:=c),return a)

This is equally correct, but less readable.

Axiom contains the usual programming constructs: loops with for and while, branching, and recursion. The above program shows several of these.

An algebraic problem

If you factorize for a few values of , each factor seems to have non-zero coefficients only of 1 or -1. This might lead us to conjecture that this is the case for all . We can write some functions and programs to test this. First, the maximum coefficient of a polynomial:

maxCoeff(p) ==
  reduce(max,map(abs,coefficients(p)))

Next, the maximum coefficient taken over all factors of a polynomial:

maxFactCoeff(q) ==
  reduce(max,[maxCoeff(nthFactor(q,j)) for j in 1..numberOfFactors(q)])

Finally, the “driver” program which simply lists all values of for which the factorization of contains a factor one of whose coefficients has an absolute value greater than 1:

testCoeff() ==
  for k in 1.. repeat
    if maxFactCoeff(factor(x^(k::PI)-1))>1 then print(k)

Put all these into a single file, called say maxFactor.input, read it into Axiom with

)read maxFactor

and run it with

testCoeff()

You will obtain the following values:

   105
   165
   195
   210
   255
   273
   285
   315
   330
   345
   357

and so on, until either you hit CTRL-c to abort the program, or you hit your computer with a heavy brick. (I suggest the first method).

This problem can be equivalently stated in terms of the coefficients of cyclotomic polynomials, and in fact the sequence generated above is sequence A013590 in the Online Encyclopaedia of Integer Sequences: “Orders of cyclotomic polynomials containing a coefficient >= 2.”