Root finding using the Bisection Method

I was motivated to write this program after joining in an Open University discussion on root finding using the bisection method (also known as the binary search method or the interval halving method). Chris has nicely described this method and the the errors in determining the root in his blog.

For completeness I will describe the method here. Suppose we wish to determine the single root of a continuous function f(x) that we know lies somewhere between x=a and x=b. That is, we want to find an approximate (or exact) value for α where f(α)=0 in the interval a<x<b.

As f(x) is continuous, the graph of y=f(x) must change sign as the graph crosses the x-axis at x=α. If we consider the value of f at the midpoint c=(a+b)/2 of the interval a<x<b, then if f(c) is the same sign as f(a), then the graph of f must cross the x-axis in the interval c<x<b whereas if f(c) is the opposite sign to f(a) then the graph of f must cross the x-axis in the interval a<x<c. Clearly if f(c)=0 then c=α and the root has been found.

Thus, even if we haven't found the root, we have at least narrowed down its location to either the interval a<x<c or c<x<b. By repeating the process iteratively we can home in on the location of the root by confining it to smaller and smaller intervals. On the first iteration our estimate of the root α is c and we know that the error in determining this root is less than (b-a)/2 since the root lies somewhere between x=a and x=b. On each iteration the maximum error will halve and so we can stop the iterations when this error is less than some tolerance δ.

Conditional branching command =>

In the multiplicative inverses program and the approximate integration program use was made of the conditional branching command => . For example, in the integration program we had the statement

:B=1=>A+2(C-1)Y→X:

This is a useful short way of coding the statements

If B=1 then
   Store A+2(C-1)Y in X
EndIf

So if the condition that B=1 is true the statement after => is executed, i.e. A+2(C-1)Y is stored in X, otherwise it is skipped and the next statement after the separator command ':' is executed.

Multiplicative Inverses in Zn - Euclid's Algorithm

Previously we saw that the multiplicative inverse of 33 in Z₂₁₁ was 32 because 33x32=1 (mod 211). I also mentioned that this inverse could be calculated using Euclid's Algorithm. Here I will explain how this is done.

We begin by writing 211 in terms of multiples of 33 and a remainder, i.e. 211=6x33+13. We then split 33 into multiples of the remainder 13, i.e. 33=2x13+7. Repeating this process we continue until the remainder is 1. This repeated use of the division algorithm gives the following set of divisions and remainders:-

211=6x33+13
33=2x13+7
13=1x7+6
7=1x6+1

We now reverse the equations and express the remainder 1 in terms of the other remainders:-

1=7-1x6
6=13-1x7
7=33-2x13
13=211-6x33

We now eliminate the remainder 6 from the first equation using the second:-

1=7-1x(13-1x7)=-1x13+2x7

We now eliminate the remainder 7 from this equation using the third:-

1=-1x13+2x(33-2x13)=2x33-5x13

Finally we eliminate the remainder 13 from this equation using the fourth:-

1=2x33-5x(211-6x33)=-5x211+32x33

Rewriting this we get 33x32=5x211+1 which is what we obtained with the fx-50F.

It is not difficult to see that such a calculation is quite time consuming and it is easy to get lost in the process. If you had quite a few of these calculations to do, then the brute force method using a calculator seems much the better option!

Multiplicative Inverses in Zn - execution of the program

This program will give the multiplicative inverse of a number in Z_n. As discussed previously suppose we wish to find the inverse of 3 in Z₇. On running the program, the dispay shows A? This is the number whose inverse we require. Enter 3 and press [EXE]. The display then shows B? This is the number n. Enter 7 and press [EXE]. The display then shows AX-C+B-1=>X and the number 5. This indicates that the inverse 5 has been found and this matches what we expect.

We can look at bigger values of n. Consider n=211 which is a prime number. As 211 is prime then any number j>1 chosen from Z₂₁₁ will be coprime with 211 and so will have a multiplicative inverse. For example, running the program with j=33 and n=211 we find that its inverse is 32 in Z₂₁₁. This is because 33x32=1056=5x211+1.

What if j and n are not coprime. Running the program with j=8 and n=32 the calculator displays X<B=>Goto 0 and 0. This indicates that no multiplicative inverse has been found and this is to be expected as 8 and 32 have a common factors of 2, 4 and 8 and so are not coprime.

Multiplicative Inverses in Zn - the program

The logic for the multiplicative inverses program was discussed previously. Here is the logic converted into a program for the fx-50F:-

?→A:?→B:2→X:Lbl 0:B→C:While C<AX:C+B→C:WhileEnd:AX-C+B=1=>X▲X+1→X:X<B=>Goto 0:

This is pretty much how the program looks when keyed into the calculator in edit mode. Special Program Commands (? → => : Lbl While WhileEnd = < ▲ Goto) are input using the [SHIFT] [P-CMD] keys. The cursor key (marked Replay) is used to switch between various screens of commands. Memory variables A, B, C, and X are input using [ALPHA] [A], [B] etc keys. The ":" Separator Code command is used frequently and can be alternatively input using the [EXE] key.

This program is 57 steps in length.

Multiplicative Inverses in Zn - the logic

So how do we go about finding, say, the multiplicative inverse of 3 in Z₇? What we want is to find the number k such that 3xk=1 (mod 7). When you have a programmable calculator, the answer can be found in a systematic way. As Z₇ is the finite set {0,1,2,3,4,5,6}, we can assign each element of this set to k and then calculate the value of 3xk (mod 7). When we find that the result is 1, then the element we assigned to k is the multiplicative inverse of 3 in Z₇.

Following this method then, the results would be:-

3x0=0 (mod 7)
3x1=3 (mod 7)
3x2=6 (mod 7)
3x3=2 (mod 7) as 3x3=9=7+2
3x4=5 (mod 7) as 3x4=12=7+5
3x5=1 (mod 7) as 3x5=15=2x7+1

So we can stop the processing at this point as we have found the solution to the problem, i.e. 3x5=1 (mod 7). So 5 is the multiplicative inverse of 3 in Z₇.

Note that we can actually start the processing with k=2 since it is always true that jx0=0 (mod n) and jx1=j (mod n) (neither of which result in the value of 1) and we assume that we are not looking for the multiplicative inverse of either 0 or 1 (since 0 has no multplicative inverse and 1 is always its own inverse, i.e. 1x1=1 (mod n))

The logic for this processing is as follows:-

Input A (the number j whose inverse is sought)
Input B (the value of n)
Store 2 in X (the first value of k to try)
Label 0
Store B in C (C holds multiples of n)
While C is less than AX (AX is the product of jxk)
   Store C+B in C
WhileEnd
If AX-C+B=1 then
   Display X (this is the multiplicative inverse of j)
EndIf
Store X+1 in X (try the next value of k in the ordered set)
If X less than B then
   Goto Label 0
EndIf