Computing totient(n) -- melding math and programming
Kirby Urner
So the idea here was to inject some computer programming
into a course with math content, that content being some
group and number theory.
We're killing two birds with one stone: learn some
programming in tandem with some math ("maths" -- UK).
Not a new idea.
What's a little bit new is that programming environments
have become easier to use, plus are freely downloadable
over the internet. We have several to choose from.
Some criteria should help us to pick an appropriate
environment. A lot of high schools seem to like Scheme.
PLT Scheme in particular is all about providing an
environment geared to students and teaching.
What content shall we cover? It will overlap quite a bit
with content already standard in middle and high school,
but then go further in some directions.
For example, the concepts of prime and composite number are
already deeply entrenched, as are even and odd. We can find
these ideas already well-formed in Euclid.
But not until Euler do we have the formal definition of
'totient': the number of positive integers < N with no factors
in common with N (i.e. relatively prime to N).
These positive integers < N and coprime to N are called
'the totatives of N'. The totient is just the number of
totatives.
Most contemporary curricula leave 'totient' aside until
some college course. But it's such a simple idea and so
perfect for learning a little programming.
We want to write something like:
totient(n): count(i) for 0<i<n where gcd(i,n)=1
How would we write this in Scheme? Here's one possible solution:
;; compute the number of positives < n that are
;; relatively prime to n
(define (totient n)
(let loop ((tot 0)(pos (- n 1)))
(if (> pos 0)
(if (= 1 (gcd n pos))
(loop (+ tot 1)(- pos 1))
(loop tot (- pos 1))
)
tot)
)
)
I'm not highly experienced in Scheme. A more streamlined
version of the totient function, even using this same algorithm,
is no doubt possible -- but even this one is pretty short.
Here's the new program in action (> means user enters):
> (totient 20)
8
> (totient 100)
40
> (totient 291929)
265380
There's a noticable delay on this last execution. The program
is basically counting down from n to 0, counting n's coprimes
along the way. As the numbers get larger, the execution time
takes proportionally longer. It's not an algorithm well suited
to large numbers.
Euler came up with a better one that requires knowing n's prime
factors -- and those may be hard to get as well, when n is trully
large.
Rather typically for Scheme, the proposed solution defines a
recursive loop, in this case using 'named let' syntax. On a
first iteration through the loop, tot is initialized to zero
and pos is set at one less than the argument n. This is the
meaning of the line:
(let loop ((tot 0)(pos (- n 1)))
...)
If pos and n are relatively prime, i.e. their gcd is one, then
the loop is called again, with tot, the number of totatives,
incremented. Otherwise, the loop is called without incrementing
tot:
(if (= 1 (gcd n pos))
(loop (+ tot 1)(- pos 1))
(loop tot (- pos 1))
)
In both cases, pos is decremented and so will inevitably reach 0.
When it does, tot is the returned value:
(if (> pos 0)
(...
)
tot)
Let's briefly look at a different language, likewise available
for free over the internet: Python. Python and Scheme both
support a "shell mode", meaning the student is able to interact
with her or his programs in one area, while defining programs
in another.
In interactive mode, these environments may be used much as
students use calculators, which may be useful for continuity.
More to the point, an interactive mode gives immediate feedback
regarding the correctness of syntax and programs.
In Python, gcd is not built in, as it is in Scheme; we need to
program it. Fortunately, Euclid's algorithm is quite compact:
def gcd(a,b):
"Euclid's Algorithm"
while b:
a, b = b, a % b
return a
Translation: as long as b is still non-zero, loop around
making b the next a, and the remainder after dividing b into
a the next b. When b becomes zero, a has divided it evenly,
so return a.
When explaining why this works, it helps to use planks. What
greatest length evenly divides two planks? "Evenly" here means
"with no leftover bit".
We assume a unit plank, and that both planks are integral
multiples of this unit.
If the shorter divides the longer evenly, we're done. The
shorter length is the gcd.
If not, then there's some remainder after dividing by the shorter,
and the gcd will have to divide this remainder, as well as the
shorter, in order to divide both the shorter and the longer (our
original goal).
So now we're interested in the gcd of the remainder and the
shorter. Notice how the remainder now plays the role of the
shorter, and the shorter plays the role of the longer. So we
divide the remainder into the shorter, to get a new remainder,
and so on. Eventually, we'll reach the gcd, which may be one.
Finally, the program for returning the totient is short:
def totient(n):
"short and sweet, but a bit slow"
return len([i for i in range(1,n) if gcd(i,n)==1])
Here's the program in action, in shell mode:
>>> totient(20)
8
>>> totient(100)
40
>>> totient(291929)
265380
Again, there's a noticable time lag as the last execution
completes.
At this point, it looks a lot like Scheme. However, under the
hood, these two languages use obviously different syntax.
We could try to imitate the Scheme solution in Python, also
using recursion:
def totient(n):
"An unPythonic implementation"
def loop(tot, pos):
while pos>0:
if gcd(pos,n)==1: return loop(tot+1,pos-1)
else: return loop(tot, pos-1)
return tot
return loop(0,n-1)
But unlike Scheme, Python does not implement tail-call elimination.
Each loop takes us deeper into the call stack, which has an upper
limit in Python. In Scheme, on the other hand, tail calls form
a genuine loop.
Here's a non-recursive approach in Python that doesn't involve
listing all the totatives first in order to count them, as did
our first Python attempt:
>>> def totient(n):
"""
Compute the number of positives < n that are
relatively prime to n -- good solution!
"""
tot, pos = 0, n-1
while pos>0:
if gcd(pos,n)==1: tot += 1
pos -= 1
return tot
>>> totient(10)
4
>>> totient(100)
40
>>> totient(291929)
265380
This is noticibly faster than the first Python solution, as it
doesn't involve keeping all the totatives in a list in order to
count them. It also looks more like the Scheme version, except
that it's not recursive.
Kirby
Scheme: TeachScheme!: http://www.teach-scheme.org/
Python: edu-sig page: http://www.python/org/sigs/edu-sig/
--
submissions: post to k12.ed.math or e-mail to k12...@k12groups.org
private e-mail to the k12.ed.math moderator: kem-mo...@k12groups.org
newsgroup website: http://www.thinkspot.net/k12math/
newsgroup charter: http://www.thinkspot.net/k12math/charter.html
Kirby Urner
>How would we write this in Scheme? Here's one possible solution:
>
> ;; compute the number of positives < n that are
> ;; relatively prime to n
>
> (define (totient n)
> (let loop ((tot 0)(pos (- n 1)))
> (if (> pos 0)
> (if (= 1 (gcd n pos))
> (loop (+ tot 1)(- pos 1))
> (loop tot (- pos 1))
> )
> tot)
> )
> )
>
>I'm not highly experienced in Scheme. A more streamlined
>version of the totient function, even using this same algorithm,
>is no doubt possible -- but even this one is pretty short.
>
Here's a slightly modified version that came to me after
I'd read the latest Harry Potter for a few more chapters:
;; totient : integer -> integer
;; compute the number of positives < n that are
;; relatively prime to n
(define (totient n)
(let loop ((tot 0)(pos (- n 1)))
(if (> pos 0)
(loop (+ tot (if (= 1 (gcd n pos)) 1 0))
(- pos 1))
tot
)
)
)
The idea is the same: call the named loop recursively
adding 1 or 0 to tot, depending out the outcome of the
gcd test. Quit when you've tested all positives down to 1
and return whatever tot has accumulated.
If you haven't played with DrScheme, why not give it a try?
It's free. The syntax checker will fight you tooth and nail
(its job) until you've started to get the hang of it --
then it'll neatly colorize your code.
Thanks to Harry Potter, I can't help thinking of learning
to write spells. Students might enjoy that analogy, even
if some oldsters become quite humorless when it comes to
mixing in such literary allusions.
Where we go with this totient concept is towards Euler's
Theorem and multiplicative groups modulo N, where the
elements being multiplied are N's totatives. You get a
group. And any of these totatives to the totient power,
modulo N, equals one. Showing what that even *means*
(by repeated examples, using one's own code) should
precede trying to prove it.
If N is a prime, then all positives < N are its totatives,
so the totient is simply p-1. We can check that using the
code we've already written (above), although it's also
intuitively obvious that primes can't have any proper
factors, so *of course* p has p-1 totatives:
> (totient 17)
16
> (totient 103)
102
Here I think that Scheme and Python approaches might diverge
quite a bit. In Python, a most promising approach is to
define "Modulo Numbers" with a built-in knowledge of how
to add, multiply etc., modulo something. We can reuse the
regular operators, * and +. So 2 + 2 might well not be 4
when the Modulo Numbers come to town.
I need to think some more about what a good Scheme approach
might be. No doubt this is already well explored in the
literature -- but I should do my own thinking as well.
And then there's always J, another language with an interactive
shell mode: http://www.inetarena.com/~pdx4d/ocn/Jlang.html
Such an embarrassment of riches! Maths mixed with programming
could turn out to be magical for at least some kids. More
students need to be at least exposed to the opportunity.
Kirby
Kirby Urner
And another thing: we need to check boundary conditions.
That means asking about totient(1) and totient(0).
And what will happen in our program if someone asks for
the totient of a non-integer or negative integer?
So it turns out there's this subtle bug in *both*
versions of the program (the Scheme version and the
Python version).
In initializing pos to n-1, where we're asking for
totient(n), I was making the assumption that we
don't want to bother with gcd(n,n) since that will
always be n, and we only care if the gcd is 1. But
what if n=1? gcd(1,1)=1 and therefore pos should
be incremented. So I should initialize pos to n,
in order to catch this case.
According to the MathWorld website, totient(0)=1 by
convention, and yet Mathematica itself (the product
around which MathWorld is built) allows EulerPhi[0]
to be zero. Note that our totient function was
symbolized by the greek letter phi by Euler, hence
this function name EulerPhi in Mathematica.
http://mathworld.wolfram.com/TotientFunction.html
In the interest of brevity, and because I don't see
why this convention is of great importance, I'll
not try to capture totient(0) as a special case
and return 1. It's OK with me for both programs
to return 0 at this point.
Finally, what if a negative or non-integer is passed
as an argument. Totient(n) is clearly not defined in
these cases. We're outside of totient's domain. It'd
be cool to trap these incorrect inputs and indicate
the error to the user somehow.
Here are enhanced versions of both programs with these
features:
Scheme:
;; totient : integer -> integer
;; compute the number of positives < n that are
;; relatively prime to n
(define (totient n)
(if (not (and (integer? n)(>= n 0))) "Incorrect input"
(begin
(let loop ((tot 0)(pos n))
(if (> pos 0)
(loop (+ tot (if (= 1 (gcd n pos)) 1 0))
(- pos 1))
tot
)
)
)
)
)
Python:
def totient(n):
"""
count totatives of n, assuming gcd already
defined
"""
if not (type(n)==type(1) and n>=0):
raise ValueError, 'Invalid input type'
tot,pos = 0, n
while pos>0:
if gcd(pos,n)==1: tot += 1
pos -= 1
return tot
Again, as per MathWorld and an earlier remark, the more
typical way to compute the totient of n is to break n
into prime factors and then "cross out" all multiples
of these primes. The totient will be the number of
integers surviving the cross-out. There's a way to
derive the number of remaining prime multiples per
each unique prime factor and to take the product of
these numbers as the total proportion of positives
=< n that will be removed. Therefore n times this
product will be the totient.
totient(n) = n * (1 - 1/p1)(1 - 1/p2)(1 - 1/p3)...(1 - 1/pm)
where p1...pm are the unique prime factors of n.
Gregor Lingl
you could fix this in the following way:
def totient(n):
"""
count totatives of n, assuming gcd already
defined
"""
if not (type(n)==type(1) and n>=0):
raise ValueError, 'Invalid input type'
tot,pos = 1, n-1
while pos>1:
if gcd(pos,n)==1: tot += 1
pos -= 1
return tot
Regards, Gregor