# Category Archives: Software

## Meeting Julia

In my last post I mentioned the new language Julia. It deserves more than a single paragraph, so I thought I’d walk through a problem, and tackle it with the language.

The problem is a stock standard one: investigating the fitting of an SIR model of disease spread to an outbreak of influenza at an English boarding school. Of the 763 pupils, on the first day only one was infected, and the number of infected pupils changed over 14 days as:

1   2   3   4   5   6   7   8   9  10  11  12  13  14
3   6  25  73 222 294 258 237 191 125  69  27  11   4


We will try to find parameters $\beta$ and $\gamma$ so that when the equations

$\displaystyle{\frac{dS}{dt}=-\beta IS}$

$\displaystyle{\frac{dI}{dt}=\beta IS-\gamma I}$

$\displaystyle{\frac{dR}{dt}=\gamma I}$

are solved, the values of $I$ fit the data.

Note that the sum of all the right hand sides is zero; this indicates that $S+I+R$ is constant, which for this model is the population. (Other models allow for births, deaths and other population changes).

Finding the parameters can be done in several steps:

1. Set up a function for the system of differential equations.
2. From that function, create a second function which is the sum of squared differences between the data, and the computed values from the equations.
3. Minimize that sum.

The first function can be created in Julia (following the example here) as:

function SIR(t,x,p)
S = x[1]
I = x[2]
R = x[3]
beta = p[1]
gamma = p[2]
[-beta*S*I, beta*S*I-gamma*I, gamma*I]
end


The default behaviour of functions in Julia is to return the line before the “end” statement; thus in a function such as this one, there is no need for a “return” statement (although one exists).

To solve a system of equations we need access to the ODE package, which I’ll assumed you’ve installed with

Pkg.add()

This can be loaded into a running Julia session either with

using ODE

which adds all the ODE functions into the top namespace, or as

import ODE

which adds all the functions in the ODE namespace. So to call the Dormand-Prince method, you would use

ode45

if all methods were in the top namespace, or as

ODE.ode45

in the second case.

So, for example, we could attempt to solve the equations as:

p = [0.01, 0.1]
t, x = ODE.ode45((t,x)->SIR(t,x,p),[0:0.1:14],[762.,1.,0.]);


To plot the functions, we’ll use the plotting library Winston (Winston, Julia from 1984 – geddit?), which is fairly basic, but enough for our simple needs:

import Winston
wn = Winston

wn.plot(t,x[:,1],"b",t,x[:,2],"g",t,x[:,3],"r",[1:14],data',"k*")


The green curve corresponds to $I$, the number of infected individuals at a given time, and the asterisks to the data. Clearly these parameters do not provide a very good fit.

The next step is to create a sum of squares function to minimize. To do this, we will solve the ODE system with the ode4 method which uses a fixed step size. This means we can ensure that there are computed values for $I$ at all the integer points:

function ss(b::Vector)
data = [3 6 25 73 222 294 258 237 191 125 69 27 11 4];
t,x = ode.ode4((t,x)->SIR(t,x,b),[0:0.05:14],[762.,1.,0.]);
sum((data-x[21:20:281,2]').^2)
end


Now this function can be optimized, using a method from the Optim package:

import Optim

Optim.optimize(ss,[0.001,1],method=:cg)


The “method” variable here defaults to the Nelder-Mead method, but for this particular function I found that the conjugate-gradient method gave the best results.

This produces quite a bit of information, of which the salient bits are:

 * Minimum: .0021806376138654117
.44517197444683443

* Value of Function at Minimum: 4507.078964


So, let’s try these parameters and see what happens:

p = [.0021806376138654117, .44517197444683443];
t, x = ODE.ode45((t,x)->SIR(t,x,p),[0:0.1:14],[762.,1.,0.]);
wn.plot(t,x[:,1],"b",t,x[:,2],"g",t,x[:,3],"r",[1:14],data',"k*")


And here’s the result:

As you can see, the fit is remarkably good.

As far as problems go, this is not particularly hard – conceptually speaking it’s probably at an undergraduate level. On the other hand, it’s not completely trivial, and does seem to me to give the computational environment a bit of a work-out. And Julia solves it with ease.

## The best Matlab alternative (3)

Over two years ago I wrote The best Matlab alternative with a follow-up a bit later, which seem to have engendered a far amount of discussion.

Well, things have moved on in the computational world, and a user is now spoiled for choice.

Octave has just had version 3.8 released; this version comes with a GUI (which is not enabled by default). However, the GUI is a very fine thing indeed:

As you see, it has all the bells and whistles you may want: command window and history, variable browser, file browser, editor, and a documentation reader. The GUI will become standard as of version 4.0.

I think that Octave has come along in leaps and bounds over the past few years, and as a drop-in Matlab replacement (if that’s what you want) I don’t think it can be bettered.

Scilab is at version 5.4.1, and is a mature product. I have come across a few comments in newgroups to the effect that Scilab is better than Octave for some advanced numerical work. However, I’m not sure if these are just opinions or are backed up with solid data and timings. I used to like Scilab a great deal, but now I don’t see that it offers anything over Octave other than Xcos (a dynamic systems modeller, similar to Simulink), and a few niche extras, such as solutions of boundary value problems.

Python was a language I hadn’t considered in my first post, but with the addition of NumPy and SciPy (plus numerous add-on packages) becomes a definite contender for serious number-crunching. The web is full of stories about people who have dumped Matlab for Python. And of course with Python you get a lovely programming language, and extensions such as Cython, which gives you C-type speed. To get an idea of the packages, see here. Python has a interactive shell IPython which adds enormously to Python’s ease of use.

Julia is the “new kid on the block”: a language which has been designed from the ground up for efficiency and speed, as well as convenience and ease of use. It is supposed to have the power of Matlab, the speed of C, and the elegance of Python. It is still in the early stages (the most recent version is 0.3.0), but it shows enormous promise, and has garnered an impressive list of users and developers. And there is already a small – but growing – list of add-on packages. The Julia web pages show some timings which seem to indicate that Julia is faster (sometimes by several orders of magnitude) than its contenders. (Well, they would, wouldn’t they? They would hardly show timings which were slower than other environments.)

My particular interest is image processing, in which Octave does very well (its image package is quite mature now), followed by Python (with its Mahotas and scikit-image packages). Scilab has its SIP toolbox, but I spent several days trying to install it and gave up. There is an Images package for Julia, but it doesn’t seem on a first glance to have the breadth of the others.

I ran a single test myself: to create a 1000 x 1000 random matrix, invert it, and find the trace of the product of the original and its inverse. In Octave:

tic; A = rand(1000,1000);Ai = inv(A);trace(A*Ai);toc

Scilab:

tic; A = rand(1000,1000);Ai = inv(A);trace(A*Ai);toc

Python (importing numpy as np):

%timeit A = np.matrix(np.random.rand(1000,1000));Ai = A.I;np.trace(A*Ai)

Julia:

tic();A=rand(1000,1000);Ai = inv(A);trace(A*Ai);toc()

Note that Octave, Scilab and Julia all have similar syntax; Python is slightly different because functions don’t all exist in the top namespace. This means, for example, you can use a function called “rand” from different packages, with different functionality, easily.

And here are the results (running the code three times each):

Octave: 1.08677; 1.04735; 1.10616 seconds

Scilab: 3.502; 3.425; 3.663 seconds

Python: 2.33; 2.28; 2.39 seconds

Julia: 0.47034; 0.35776; 0.36403 seconds

This is hardly definitive, but it does show how fast Julia can be, and in relation how slow Scilab can be.

These four environments are by no means the end. R is extensively used, and has over 5000 add-on packages at the CRAN repository. However, I’ve never used R myself. I believe it has branched out from its statistical foundations, and is now more of a general use environment, but still with the best statistical functionality in the business. Then there’s the Perl Data Language, apparently used by astrophysicists, and the GNU Data Language, about which I know nothing. See the Wikipedia page for a list of others.

(Note that I’ve not discussed Freemat this time; I don’t believe it’s a serious contender. I had a look through some of its source code, and most of its files seem to be cut-down Matlab/Octave files without any error checking. Anyway, there seems to be no need to use niche software when Octave is so mature and so well-supported by its developers and users.)

My recommendations? Use Octave, unless you have some niche requirements, in which case you may find Scilab more suitable. If you’re prepared to sacrifice maturity for speed, give Julia a go. If your wish is a mature programming language as your base, then Python.

## A first look a WeBWorK

In my last post (nearly three months ago) I commented on online assessment, and in particular on Pearson’s MyMathLab Global, which we have been using for our engineering mathematics. This is quite a good, full-featured and robust package; its only fault (if indeed a fault it be) is that it’s commercial. This works fine if you are structuring your course around a set text, but over the years I have moved further and further away from such texts. I have not found a text which could be used comfortably with my extraordinarily heterogeneous student cohort, who enter our courses with all different levels of mathematics competency, and all different levels of written and spoken English.

A good text must be carefully, simply, and clearly written, with plenty of useful examples, and above all well scaffolded exercises which move gently from the very simple to the complicated. Most (all?) textbooks in my experience suffer greatly from this. Our current preferred textbook: “Modern Engineering Mathematics” by Glyn James, has quite appalling exercises. Also, it’s bit of a heavy monster of a book, like most mathematics textbooks, and I don’t see why students should have to lug it around with them. But anyway, if you do buy the book, either in hardcopy, or in electronic form, you get a code which gives you one years access to MyMathLab and all ts trimmings.

So this last semester I have used my own notes and exercises for engineering mathematics, and sent students off to MyMathLab for their weekly online tests. There were problems:

1. Students who first started over a year ago (who, for example, might have had to take a break between their mathematics subjects) found that their access codes had expired. I spent the first few weeks of the semester writing almost daily emails to the publishers begging for some more access codes. To their credit, the local Pearson reps were in every way helpful, accessible, friendly, and understanding.
2. Right towards the end of the semester I kept finding students who had fallen through the cracks, so to speak, and hadn’t been able to access their tests at all. (And had got busy with other subjects and kept forgetting about them)
3. Some answers were marked wrong, when they weren’t. I had one student who sent me a screen shot to show that he had entered “12″ as the numerical answer to a problem, which was marked wrong: the correct answer, apparently, was “12/1″. There was also confusion with decimal places: a question would ask for say 4 correct decimal places, and no matter how carefully the students entered their answers they were marked wrong.

I have wanted to experiment with WeBWorK for a while now, but I couldn’t find a way to install it locally – until I realized that I could install it (and a web server) on my desktop at work, which runs Ubuntu 12.04. I did this (it took most of a day, with only a modicum of swearing, and messages to the WeBWorK forums), but now I have a working system, even if only running on a desktop, rather than on a dedicated server.

WeBWorK has many excellent features:

1. It’s free: open source, no less. You can pay for MAA to host your course if you don’t want to set up your own local service, but even that option is very inexpensive.
2. It has a huge problem library: the Open Problem Library has something around 25,000 problems from most areas of undergraduate mathematics. You might not find problems on, say, cohomology or modular forms, but these advanced topics would hardly be well served by an online homework system. We are assuming here we have a large class (in the many 100′s) of students in their first or second year of university study. And for this sort of basic mathematics, WeBWorK is terrific.
3. The authoring system for creating new problems is very straightforward: just a matter of some elementary editing. I tried authoring a problem in MyMathLab once, and although it’s quite possible, life’s too short.
4. WeBWorK seems to play very nicely with LaTeX and Sage, so that the mathematics is properly typeset, and you can have all the power of Sage available to you if you need it.

That’s as far as I’ve got so far. You can get a feel of WeBWorK without installing it by checking out some of MAA’s “Model Courses”, of which a list is here. There doesn’t seem to be any online hosted method for browsing the problem libraries to see what’s in them, as far as I know. I think this could be something worth making available so that teaching staff can decide whether they think WeBWorK would suit them and their courses.

In my current state of thinking, if I can get a better hosted system at my work (that is, not my desktop), I may well use WeBWorK for my next teaching semester.

## Online assessment in mathematics

This semester I’ve taken over a large first year subject, for which the previous subject convenor had organized to use MyMathLab Global for weekly testing.  The subject is based around Glyn James’ “Modern Engineering Mathematics”, a book which is OK for content, and pretty awful for exercises.  This means that as users of the text, we have access to the publisher’s (Pearson) online testing service.  For educators the idea is terrific: every week I simply pull from the exercise bank a set of 10 exercises corresponding to our current topic, and make them available for a week during which time students can have unlimited goes at them.  So in theory it’s an easy way for students to get some easy marks, and reduce the burden of marking weekly tests by the lecturer (me) and the subject tutors.

The subject I am teaching “Engineering Mathematics 2″ is a follow on subject from – wait, you;ve guessed it! – “Engineering Mathematics 1″.

What could be better?  Well, aside from the extraordinary ease of testing, I have begun to have doubts about the efficacy of the system.

1. It’s a commercial system, which means that students have to buy a “personal access code” to use it.  A code comes with the book if they purchase it.  However, a code lasts for only 12 months, and if students have bought the book (and its code) for Eng Maths 1 (as they all should), and if there’s been a break in their studies for any reason, then their code may have expired by the time they come to Eng Maths 2.  Then there are all the students who have got subject exemptions from Eng Maths 1 and never had a code in the first place.  The local publishers reps have been terrific and have provided me with lots of extra codes to hand out, but the onus is on me to get these dratted codes in the first place, and ensure the students get them.
2. The system only tests the final answer.  A typical question is:

Find the partial derivatives

$\displaystyle{\frac{\partial f}{\partial x}}$, $\displaystyle{\frac{\partial f}{\partial y}}$ and $\displaystyle{\frac{\partial f}{\partial y}}$ of $f(x,y,z)=\sin^{-1}(xyz)$.

and the students get three entry boxes for their results.  MyMathLab Global is very picky about the form of a solution, and if the form of a student solution differs from the “correct one” they are marked wrong.  For example,  we find that

$\displaystyle{\frac{\partial f}{\partial x}=\frac{4yz}{\sqrt{1-16x^2y^2z^2}}}$.

if a student leaves out the “1″ of the “16″ in error, the answers is just marked wrong.  If the student decides to write

$\displaystyle{\frac{4yz}{\sqrt{(1-4xyz)(1+4xyz)}}}$

it is marked wrong.  MyMathLab Global doesn’t seem to include a CAS to check if two answers are equivalent (aside from some very simple arithmetic and operator commutativity).  This has been a source of annoyance to my students, who may enter an answer which is mathematically correct, and yet marked wrong.  A very spirited attack on MyMathLab Global can be seen here.

3. The system checks only the final answer. We might well ask: what mathematical learning is actually being tested here? Surely we want our students to show that they have mastered a mathematical technique, process, or algorithm, and we would mark them on their working, not just on the final answer. At least, that is how I mark by hand. A question for which the student’s working is fundamentally correct but the final result is wrong will still get most marks from me. I regard the final result as just one part of a question. So suppose the student sets out to answer a question, and makes a small slip somewhere along the way which effects every succeeding step, and the final result. I would give almost compete marks – there’s only been one very small slip, and aside from that the logic and the working are exemplary. But MyMathLab Global can’t do this. An answer is either correct, or wrong.

I have been hoping to check out MAA WeBWorK as a comparison, but the only server I have access to runs a linux distribution (RHEL 5.9) which can’t run it. I think WeBWorK, from my understanding, is rather more nuanced in terms of mathematical equivalence of expressions, so probably overcomes item 2 above, and being open source completely overcomes item 1. I don’t think it overcomes item 3, though.

I will continue to use MyMathLab Global – its convenience is simply too great – in spite of my misgivings. But I think that online assessment in mathematics is still a very long way from testing the students’ ability to do mathematics, as opposed to their ability to obtain an answer. This latter is not mathematics: it is applied arithmetic; and noble and useful as it may be, it is only a tiny part of what we, as mathematics educators, should be testing.

## Easy Simpson’s rule

Simpson’s rule, as I’m sure all readers of this blog will know, is a simple rule for numerical integration (or quadrature). It can be obtained in two ways, either by determining the interpolating polynomial over three equally spaced points and integrating that polynomial, or by finding the values $w_i$ for which

$\displaystyle{\int_a^bf(x)\,dx = w_0f(a)+w_1f(a+h)+w_2f(b)}$

(where $h=(b-a)/2$) is accurate for $f(x)=1,x,x^2$.  Over the interval $[0,2]\,$ Simpson’s rule is

$\displaystyle{\int_0^2f(x)\,dx = \frac{1}{3}\Bigl(f(0)+4f(1)+f(2)\Bigr)}$

and over a general interval

$\displaystyle{\int_a^bf(x)\,dx = \frac{b-a}{6}\Bigl(f(a)+4f(a+h)+f(b)\Bigr)}.$

Simpson’s rule has the added advantage that it is also accurate for cubic polynomials.

On its own the rule is not hugely accurate, so one way of increasing the accuracy is to divide the interval into smaller, equally sized regions and apply Simpson’s rule to each of them. This provides the composite Simpson’s rule:

$\displaystyle{\int_a^bf(x)\,dx \approx \frac{b-a}{3n}\Bigl(f(a)+4f(a+h)+2f(a+2h)+4f(a+3h)+\cdots\Bigr.}$

$\displaystyle{\Bigl.+2f(a+(n-2)h)+4f(a+(n-1)h)+f(b)\Bigr)}$

where $n$ is even and $h=(b-a)/n$. Sometimes this rule is written

$\displaystyle{\int_a^bf(x)\,dx \approx \frac{b-a}{6n}\Bigl(f(a)+4f(a+h)+2f(a+2h)+4f(a+3h)+\cdots\Bigr.}$

$\displaystyle{\vspace{1cm}\Bigl.+2f(a+(2n-2)h)+4f(a+(2n-1)h)+f(b)\Bigr)}$

where there is no restriction on $n$ and $h=(b-a)/2n$. For a spirited criticism of this rule, both from practical and pedagogical perspectives, see here.

However, in a standard generic mathematics course, Simpson’s rule is as far as most students get with quadrature. The trouble is that the order of the weights (the coefficients of the function values):

$1,4,2,4,2,4,...,2,4,1$

is tricky to program. Either you have to run two loops, one for the 4′s and the other for the 2′s, or you have to do a little test at each stage to see whether you are at a 2 or 4 point, or you might invoke the modulus function. None of those ways are particularly straightforward, especially if, like many beginning students, you have little programming experience.

So here’s a super easy way. We start with end points $a$ and $b$ and an even integer $n$. We first create a list of all the values at which the function is to be computed:

$xs=[a+kh \mbox{\; for\;} k=0,1,2,\ldots n]$

where $h=(b-a)/n$, and then a list of the weights:

$w=[3-(-1)^k \mbox{\; for\;} k=0,1,2,\ldots n].$

This last list is $[2,4,2,4,\ldots,2,4,2]$.

We then compute

$\sum w_if(xs_i)$

However! note that our list of weights starts and finishes with 2′s instead of 1′s. We correct for this by simply subtracting $f(a)$ and $f(b)$ from the sum. Our “easy” Simpson’s rule is then:

$\displaystyle{\frac{b-a}{2n}\Bigl(\sum w_if(xs_i)-f(a)-f(b)\Bigr)}.$

This is almost trivial to compute. For example, on a TI-nspire CAS calculator, we might have:

$a:=0$

$b:=1$

$\mbox{Define\;}f(x)=e^{-x^2}$

$n:=16$

$h:=(b-a)/n$

$xs:=\mbox{seq}(a+k\cdot h,k,0,n)$

$w:=\mbox{seq}(3-(-1)^k,k,0,n)$

$(\mbox{sum}(w\cdot f(xs))-f(a)-f(b))\cdot (b-a)/n$

The TI-nspire has the nice attribute that a function applied to a list will produce a list of all function values, and the product of two lists is a list containing all the products of corresponding elements. The Casio Classpad has the same functionality with lists, and the above commands can be used.

This approach of using $3-(-1)^k$ to obtaining the weights for Simpson’s rule seems not widely used, but it seems to me to be absurdly simple. I’ve only come across it in one other place.

## Initial comparison of TI-nspire CAS and Casio Classpad 330 calculators

The TI-nspire CAS and the Classpad 330 are the top of the line CAS calculators, and are used in both secondary and tertiary education worldwide.  And each one has its passionate adherents and equally passionate detractors.  However, it’s hard to find a balanced comparison (would you expect to find the TI calculator favourably reviewed in a Casio forum?) TI have released a color version: the TI-nspire CX CAS; and Casio have announced a color version of the ClassPad, to be available in the second half of this year.

For the purposes of this comparison, I’ll use the current monochrome TI-nspire CAS, and the current ClassPad 330.  Here is what they look like:

Also this post will be entirely based on using the handheld calculators, not the computer software which is also available for each of them. And it will be by its very nature superficial; I hope to cover particular topics in more detail in later posts.

Interface

The TI-nspire has a screen size of 320×240 pixels; the ClassPad 160×240, exactly half the size of the TI-nspire. This means that the TI-nspire display is sharper. The TI-nspire has a menu driven system entirely driven by button presses; the ClassPad uses a touch screen controlled with a stylus. I find the TI-nspire much easier to use; I don’t particularly like jabbing the stylus at a screen with such poor resolution. It’s worth noting as a side issue that the stylus is easily lost; they seem to get looser in their housing with age. And many second hand ClassPads are advertised with “stylus missing”.

The ClassPad also defaults on readability. Here’s a couple of screenshots from Casio:

Note that at the top of the screen there is a menu bar, underneath of which is a “toolbar”. At the low resolution of the screen, I find the toolbar icons to be all but unreadable. Some of them are menu headers, which when selected produce a tools menu. This is a list of tiny icons, which to my eyes could mean anything. (Even with my reading glasses, and my eyesight, though not as sharp as once it was, is not that bad.) It may have been a cost issue: producing a touch screen cheaply, but I find the low resolution of the ClassPad a real hindrance to its use.

Commands on the TI-nspire are usually obtained with at most 4 keystrokes: a menu key, then at most 3 menu choices. Commands on the ClassPad, resolution notwithstanding, are also fairly easy to obtain, either with the menu bar, or from the onscreen keyboard. This image:

shows that Casio have in fact designed a very clever interface which makes good use of the limited space to obtain most commands and functions easily.

For comparison, here is a TI screenshot, showing a linear function and some values:

Note that the screen, though maybe cluttered for some, is still quite readable and sharp.

The TI-nspire uses a “Documents” system, and a document may have any number of different “pages”: calculation, geometry, lists & spreadsheet, data & statistics and so on. Variables may be transferred between pages, so that a variable, say a list created and named in a spreadsheet, can be opened in a calculator page. The ClassPad has “eActivities”; here is a description from a Casio page:

An eActivity is a new, exciting innovation. It is small packet of learning that can be written by anyone, even yourself! It is loaded into the calculator and a student is then free to work through it. It is normally a scripted learning experience that has all the functions of the calculator at its disposal. eActivities can be used to present students with assignments, learning opportunities or in fact to build your own simple software…

Mathematical strength

I think that the calculators are much of a muchness here. True, the ClassPad offers some goodies such as the Laplace and Fourier transforms, which might give it a slight edge. However, there are also a few anomalies. The TI-nspire provides LU factorization of a matrix $M$ in the form of the matrices $P,\,L,\,U$, where $P$ is a permutation matrix for which $PM=LU$; the ClassPad just returns $L$ and $U$. This means that given a matrix $M$ the ClassPad will return $L$ and $U$ matrices for which possibly $LU\ne M$. An example is the matrix

$\displaystyle{M=\begin{bmatrix}0&4&3\\ 1&2&-1\\ 3&1&2\end{bmatrix}}$

for which

$\displaystyle{P=\begin{bmatrix}0&0&1\\ 1&0&0\\ 0&1&0\end{bmatrix},\,L=\begin{bmatrix}1&0&0\\ 0&1&0\\ \frac{1}{3}&\frac{5}{12}&1\end{bmatrix},\,U=\begin{bmatrix}3&1&2\\ 0&4&3\\ 0&0&-\frac{35}{12}\end{bmatrix}}$

This may be a small worry, but surely it’s a simple thing and worth getting right?

I haven’t applied (as yet) a CAS test-suite to either calculator. You have to remember that in terms of raw computing power, these calculators are lightweights, and can only be expected to solve relatively small problems. But here’s one test I did yesterday. On the TI-nspire, I found the factorization of Cole’s number:

$2^{67}-1=193707721 \cdot 761838257287$

performed almost immediately. Impressive! Just to check this wasn’t a fluke, I obtained the following factorization:

$24!-1=625793187653 \cdot 991459181683$

in a little under three minutes. I tried the first factorization on the ClassPad and got nothing in a reasonable time (several minutes). So either the TI-nspire uses a better factorizing algorithm, or it has a more powerful processor.

This is an issue I’ll try exploring in more detail in later posts.

Variables and programming

The TI-nspire has three ways of assigning a value to a variable:

• Define x=3
• 3 $\rightarrow$ x (the arrow is obtained by using the “sto” key)
• x:=3 (The “:=” is obtained with two key presses: crtl, :=)

• 3 $\Rightarrow$ x
• x:=3

without, however, any single key to obtain “:=”. And only the arrow method can be used in programs.

Both calculators use a version of Basic on the handhelds. Apparently the scripting language lua can be used with both, but I suspect only on the associated computer software, and not directly on the handhelds. And they both provide the full gamut of control statements (loops, tests etc). Both calculators distinguish between functions and programs, although in slightly different ways.

The ClassPad offers very nice sequence handling: you can set up a number of interlinked sequences, and set ‘em off. For example, you can enter

$\displaystyle{\begin{array}{rl}a_{n+1}&=(4-b_n+c_n)/4\\ a_0&=1\end{array}}$

$\displaystyle{\begin{array}{rl}b_{n+1}&=(2-2\times a_n-c_n)/5\\ b_0&=1\end{array}}$

$\displaystyle{\begin{array}{rl}c_{n+1}&=(11-a_n)/3\\ c_0&=1\end{array}}$

(which are the sort of equations which may be encountered in Jacobi’s iteration method for linear equations), in the Sequence page, and just run them. The TI-nspire has a “Lists & Spreadsheet” page, in which you can perform similar calculations, but not as elegantly as this.

As a mini-comparison, here are for-loops on each:

 TI-nspire ClassPad For i,1,10x:=x+1EndFor For 1$\Rightarrow$i to 10x+1$\Rightarrow$xNext

Graphics

Both calculators support 2D graphing (they would hardly have a market if they didn’t). The ClassPad also supports 3D graphing, which is supposed to be a huge feather in its cap, but in fact I find this to be of limited use. Just as an experiment, I attempted to plot the monkey saddle surface $z=x^3-3xy^2$ and couldn’t make any sense of the result. Even with zooming, moving the graph around with the stylus, I couldn’t obtain any intuitive feel of the shape of the graph. Partly this is because there doesn’t seem to be any hidden line algorithm available. However, the drag and drop capability of the ClassPad, with the stylus, is one area which it shines. The TI-nspire’s touchpad, although it works, is clumsy in comparison.

Which is better?

I much prefer (for my purposes) the TI-nspire’s sharper screen and interface. However, when the new colour model fx-CP400 is released I will certainly reconsider. I have no loyalty to either brand. In the end, it comes down to what you and everybody else is using. My daughter is using a ClassPad 330 for her year 11 mathematics (its use is mandated by her school) and is very satisfied with it. I would imagine that even with the models discussed in this post (the monochrome TI-nspire already superseded, the ClassPad 330 about to be), either one would happily manage pretty much all school, and most college, mathematics.

## Maxima on Android

Such a long long time since I last posted!  Very poor blogging etiquette.  Anyway, to break the drought, here’s news of the best thing to come my way in the last six months – fully featured Maxima running on Android (this screenshot from its website):

At a conference I attended late last year, sponsored by calculator companies, the Maple reps (both of whom I knew), claimed that calculators were an old-fashioned outmoded technology, especially now that with smartphones, ipods etc, you can have access to the Internet in your pocket (Wolfram Alpha, anybody)?  I still think there’s a lot to be said though for a device which you can use independently of a network, and in fact I believe most calculator users use very little of their functionality.  How many CAS calculator users, I wonder, program on or with their machines?

But now here’s a fully featured, mature CAS running on a portable device – the first such one, as far as I know.  There are in fact other CAS’s available for portable systems: PocketCAS, Mathomatic, Yacas (which I’ve written about before), no doubt others.  None of them, as far as I can tell, work particularly well, or are in any way complete.  It doesn’t take long before you run up against a wall of something they are unable to do.

So I see Maxima as being the first complete portable CAS – and I’m delighted to see it on Android.  It works fine on a smartphone, and is quite wonderful on a tablet.

## The best Matlab alternative (2)

A year or so ago, I published a blog post called The best Matlab Alternative.  In that post, I discussed the merits of Scilab, Octave and Freemat, and decided that for me, Octave suited me best.  However, Scilab is probably the more full featured and powerful.

I’ve since discovered that there’s another contender, in some ways better than all of the above: Python.  My requirement for a “Matlab alternative” was that:

1. It must install under Windows. This is absolutely necessary for me, at least, since most of my students use Windows. A very few use OSX, and hardly any use Linux.
2. It must, as far as possible, be identical in use with Matlab. The commands should be the same, and simple scripts and functions should transfer between the two with no effort.
3. It should have a reasonably nice interface, as similar as possible to Matlab’s own. A basic command line only interface running in a terminal won’t cut it with my students.

Suppose we drop condition 2, and require software which can do pretty much most of what Matlab can do, in a nice visual environment, we find Python coming out looking very good indeed.  Now, I have been a lover of Matlab, and its open source alternatives, for many years, and I do like the rapid prototyping possible with Matlab, as well as its ease for dealing with large data types.  And vectorization (where a single command is automatically applied to every element of a multi-dimensional array) provides some wonderfully elegant and simple code.  (It can also of course provide for dense obtuseness.)

However, I have also become aware of some criticisms of Matlab.  A major one is vendor lock-in: the code, the language, the IDE, is all owned by The Mathworks.  If they decide to radically change something, you have no alternative but keep paying the license fees to keep up.  Second is that the language has no object-oriented capabilities, at least in a natural way;(see comment below) that the language has no standard as such – it’s just whatever The Mathworks decide to provide.  Thirdly is the namespace: all functions exist as top-level objects, and all functions (“m-files”) in the Matlab path (which usually will include multiple directories) are available when you start up.  Fourth is that an m-file can only contain one function (and maybe some sub-functions), whose name must be the same name as the file.  This means that if you create some new functions – and most users do – you need to be sure that it doesn’t conflict with any other functions.  This problem is exacerbated if you include lots of toolboxes.  The Image Processing Toolbox alone adds 270 new functions to the Matlab namespace  and there are 38 other toolboxes available.  That’s a lot of functions.

One list of criticisms is available here; another (called “Abandon Matlab”) is here.  Note that the second site recommends (as well as Python), R, which the first site also pans.  There’s no accounting for tastes.

So, back to Python.

Python is a deservedly popular language, supported by a huge user base, tons of books, articles, websites, conferences, seminars, and so on.  It is fully open source, and it is managed by a non-profit organization called the Python Software Foundation.

Python has lots of libraries for scientific computing, starting with numpy (“numeric python”) and scipy (“scientific python”).   And there are others – such as scikits-image which extends scipy’s image processing capabilities (the Python Imaging Library does not seem to be currently in active development).  Many more are listed here.

What’s more, Python has lots of IDEs which mean that you can use it pretty much as you would use Matlab, Octave or Scilab.  For a nice visual IDE with all the trimmings, there are Eric IDE, Spyder, or PyDev, which is an Eclipse plugin, plus lots more, both open source and commercial.

For those of us who are satisfied with less, there’s iPython, an enhanced shell (which is what I like to use).  And there are again lots more, including editors and all sorts of other goodies for development, here.

Let’s finish with a tiny example: to create a 1500 x 1500 random matrix, invert it, and check the trace of the product of the matrix and its inverse.  I’ll do all this in the iPython shell, which I invoke from my linux command line with ipython.

Once in the iPython shell, we import the pylab modules; this includes numpy, scipy, and the plotting library matplotlib:

In[1]: from pylab import *


(In fact this can be done from the command line with iPython, but this command should work in any Python environment.)

Now to define the matrix, compute its inverse, and check the trace of their product:

In [2]: m = mat(rand(1500,1500))
In [3]: mi = linalg.inv(m)
In [4]: np.trace(m*mi)
Out[4]: 1500.0


And of course all this is nice and fast (on my linux laptop):

In [5]: %time m = mat(rand(1500,1500)); mi = linalg.inv(m); np.trace(m*mi)
CPU times: user 6.52 s, sys: 0.21 s, total: 6.73 s
Wall time: 3.78 s


Note: there is now a follow up to this post called (now this is original and clever): The best Matlab alternative (3).

## Fast modular exponentiation

I’m a “program first, think later” kind of bloke, and I often write programs which are very slow, until I see some way later of vastly improving their speed. This was the case recently, when I was experimenting with the PRNG

$x_0=1,x_1=1, x_i=7^{x_{i-1}}+7^{x_{i-2}}\pmod{p}$

with $p=2^{31}-1$.

Modular exponentiation is computationally expensive. If you use the standard square and multiply algorithm, then you need a squaring for every bit, and an additional multiplication for every 1. In a 32 bit number (in which we would expect on average 16 ones), this is 48 square/multiply operations for each exponentiation. This meant that putting my PRNG though the TestU01 “Big Crush” suite took over 188 hours. (I just had it grind away in the background.)

However, since our exponents $n$ are always 32 bits, we can write

$n=aT+b$

where $T=2^{16}$ and $a,b< T$. Given $n$, the values $a,b$ can be found with simple bit operations:

a = n>>16; b = n&0xFFFF;

Then

$7^n=7^{aT+b}=(7^T)^a(7)^b$.

Suppose we have two arrays H[] and L[] each if size T, and containing all the powers of $7^T$ modulo $p$ and $7$ modulo $p$. Since

$7^{2^{16}}=118315438\pmod{p}$

we can precompute H and T using the following pseudo-code

H[0] = 1
L[0] = 1
for i from 1 to T-1:
H[i] = 118315438 * H[i-1] (mod p)
L[i] = 7 * l[i-1] (mod p)


Then

$7^n = H_aL_b$

with a and b computed as above. This reduces the exponentiation from 48 multiplications/squarings to just 1. When I ran Big Crush again on my speedier implementation, it took just 11 hours.

For those of you who like such things, here is a C program to test it. On my laptop it can produce 10^8 random values in just over 9 seconds.


/*
The Fibonacci exponential generator:

x[0]=1, x[1]=1, x[i] = 7^x[i-1]+7^x[i-2] (mod p) for n, with p=2^31-1.

Compile with

gcc testexpfib.c -o testexpfib -lm -std=gnu99 -g
*/

#include <stdlib.h>
#include <stdio.h>
#define p 0x7FFFFFFF
#define T 65536
#define ulong unsigned long

ulong n,a,b,i,xH,xL,yH,yL,t1,t2,t3,temp,x,y;
ulong high[T],low[T];

ulong modp(ulong z) {return (((z)&p)+((z)>>31));}

ulong mod2p(int q, ulong x) // returns (2^q)x mod p
{
return ((x>>(31-q))+((x<<q)&p));
}

ulong mul_modp(ulong x,ulong y) // returns xy mod p
{
xH = x>>16;
xL = x&0xFFFF;
yH = y>>16;
yL = y&0xFFFF;
t1 = mod2p(1,xH*yH);
t2 = modp(mod2p(16,xH*yL)+mod2p(16,xL*yH));
t3 = xL*yL;
return (modp(modp(t1+t2)+t3));
}

ulong expmodp(ulong x) // returns 7^x mod p
{
a = x>>16;
b = x & 0xFFFF;
return (mul_modp(high[a],low[b]));
}

ulong expfibp(void)
{
temp = modp(expmodp(x)+expmodp(y));
x = y;
y = temp;
return (y);
}

int main(int argc, char *argv[])
{
high[0]=1;
low[0]=1;
for (i = 1; i < T; i++)
{
high[i] = mul_modp(118315438,high[i-1]);
low[i] = mul_modp(7,low[i-1]);
}

if (argc<2)
n = 1000000;
else
n = strtol(argv[1],NULL,10);

x = 1;
y = 1;

for (i = 0; i < n+1; i++)
{
ulong s = expfibp();
if (i%(n/10)==0)
printf("%10lu, %lu\n",i,s);
}

return 0;
}


## A conceptually simple PRNG

I like pseudo-random number generators, and I’ve been experimenting with them for the last few weeks. Now, Stephan Mertens has developed a PRNG which he calls YARN (Yet Another Random Number), defined in general by

$x_i=ax_{i-1}\pmod p$

$z_i=g^{x_{i-1}}\pmod p$

$u_i=z_i/p$

where $gcd(a,p)=1$ and $g$ is a primitive root modulo $p$. You can read his material here. The final values $u_i$ will be uniformly randomly distributed over the interval $(0,1)$. The exponential definition of $z_i$ removes the linearity of the $x_i$ values. It is clear, however, that the period can only be at most $p$: after at most $p-1$ iterations, the $x_i$ values will start repeating, and hence so will the $z_i$ values.

A better (but slower) PRNG is defined as

$x_0=7, x_1=7, x_i=7^{x_{i-1}}+7^{x_{i-2}}\pmod p$ for $i\ge 2$

with $p=2^{31}-1$. This simple PRNG (which I’ve called a “Fibonacci exponential generator”), passes all the Crush tests of the TestU01 suite of randomness tests, currently the most stringent known. Here’s a little C program called fexp.c, which uses the FLINT libraries, to compute and display some of these values:

/*
A program to test the Fibonacci exponential  PRNG

x[i] = b^x[i-1] + b^x[i-2] mod p,

with b = 7, p = 2^31-1.

Compile with:
gcc fexp.c -o fexp -lflint -lmpir -ltestu01 -l probdist -lmylib -l -std=gnu99 -g
*/

#include <stdlib.h>
#include <stdio.h>
#include <flint.h>
#include <fpmz.h>
#define ulong unsigned long
#define pp 0x7FFFFFFF

fmpz_t x,y,new,r1,r2,b,p;
ulong i,n,t;

ulong fexp(void)
{
fmpz_powm(r1,b,x,p); // r1 = b^x (mod p)
fmpz_powm(r2,b,y,p); // r2 = b^y (mod p)
fmpz_add(new,r1,r2); // new = r1 + 22
fmpz_mod(new,new,p); // (mod p)
fmpz_set(x,y);       // x = y
fmpz_set(y,new);     // y = new
return (fmpz_get_ui(y));
}

int main (int argc, char *argv[])
{
if (argc < 2)
n = 1000000;
else
n = strtol(argv[1],NULL,10);

fmpz_init_set_ui(b,7);  // Initialize values for use in Flint functions
fmpz_init_set_ui(p,pp);
fmpz_init_set_ui(x,7);
fmpz_init_set_ui(y,7);

for (i = 0;i < n+1;i++)
{
t = fexp();
if (i%(n/10)==0)
printf("%10lu,%10lu\n",i,t);
}

return 0;
}


It is known that if

$x^r+x^s+1$

is a primitive polynomial, (with $r>s>0$), then the “lagged Fibonacci generator”

$x_i = x_{i-r}+x_{i-s}\pmod m$

where $m=2^M$ will have a maximum period of $(2^r-1)2^{M-1}$. See this page of Richard Brent and the wikipedia page for more details. So, if for example we set $(r,s) = (6972593,3037958)$ and $M=32$, then the resulting generator will have a maximum period of

$(2^{6972593}-1)2^{31}\approx 2^{6972624}$

which is far far greater than the period of the Mersenne twister, one of the most popular PRNGs in current use. At the cost of keeping a state vector of 6972593 4-byte values in memory (roughly 28 megabytes), then a single PRNG defined by that recurrence would exceed all the world’s current and projected uses for random numbers. However, lagged Fibonacci generators pass some tests of randomness, but not all. And partly this is due to the strictly linear nature of the recurrence.

To remove the linearity and apply the exponential generator, suppose we have a finite field $GF(2^{32})$. If $g$ is a primitive element of that field, then we may define, say

$x_0=1, x_1=x, x_i=g^{x_{i-r}}+g^{x_{i-s}}$ for $i\ge 2$

Such a PRNG is both extremely robust (passes very stringent tests of randomness), and assuming that the period is obtained from the base recurrence, has a very long period.