Is it prime?

I wanted to create an algorithm that could be used to test whether a number is prime or not. I have previously written a program that finds the prime factors of a number and this can be used to determine whether a number is prime or not but it isn't very quick for this purpose. I wanted a faster algorithm that could give a yes/no answer for all positive integers up to about 1000 and could be extended to larger numbers if required.

One issue with the fx-50F calculator is that there is no quotient/remainder function. So if we divide 101 by 7 then we know that the quotient is 14 and the remainder is 3 (i.e 101=14x7+3). If the remainder was 0 then we would have known that 101 was divisible by 7.

One way round this is to use the sine function. We have that sin(πx)=0 if x is an integer and not otherwise (you need to use the fx-50F in radians mode - [SHIFT] [SETUP] [2]). Thus if we let x=101/7=14.428..., then sin(πx)=0.974...whereas if we let x=99/3=33, then sin(πx)=0. So we have a way of determining whether one positive integer is divisible by another.

This forms the basis of the processing. For a given number N we can trial divide N by each of the prime numbers P which are less than some prime M, say. If N is divisible by P then we can set a flag to indicate that N is not prime. Clearly, if N is one of the prime numbers P, then we need to check this and flag it accordingly. This method will be accurate as long as N is is less than M².

The list of programs so far

I thought I would, in a brief blog, list the five programs that I have written for the fx-50F so far:

• Root finding using the bisection method
• Multiplicative inverses in Zn
• Approximate integration
• Prime factors
• Temperature conversion

Root finding using the Bisection Method - exact solutions

As a final blog on this program, I want to cover the situation where an exact solution of f(x)=0 is found by the bisection method, i.e. where the program finds α where f(α)=0 and not just an approximation to α.

Consider the quadratic equation 160x²-428x-645=0. This equation has a root around x=-1 because f(-1.1)=19.4 and f(-0.9)=-130.2 (where f(x)=160x²-428x-645), so f changes sign between -1.1 and -0.9. If we now edit the program to have 160C²-428C-645 instead of C-COS(C) in the "subroutine" then on running the program with A being -1.1, B being -0.9 and D set to 1x10^-5 we obtain the root -1.075 in 3 bisections. This is the exact value for the root because the quadractic equation can be factorised to (x-3.75)(x+1.075)=0, so the roots are x=-1.075 and x=3.75.

This shows that the program can cope with the situation where it finds the root exactly. Of course, this doesn't always happen even if an exact value could be found. For example, rerunning the program with A set to -1.5, B set to -0.5 and D as before, we obtain C as -1.075004578 in 17 iterations.

Root finding using the Bisection Method - execution of the program

In a previous post I provided a program for finding the root of a function f(x) using the bisection method. To show how this works, consider the problem of calculating where the lines y=x and y=cos x intersect on a graph of y against x (see the above plot from Wolfram Alpha). There is a single point of intersection somewhere between x=0 and x=π/2 and to find it we need to obtain the solution of x=cos x, i.e. x-cos x=0. This trigonometric equation has no obvious solutions and so we can use the bisection method to find a numerical solution (where f(x)=x-cos x).

The program given is already set up for this. On running the program the display shows A? This is the lower limit of the interval that is being bisected. Enter [0] and press [EXE]. The display then shows B? This is the upper limit of the interval. Enter [SHIFT] [π] [/] [2] and press [EXE]. The display then shows D? This the error tolerance. Enter [1] [EXP] [-] [5] and press [EXE].

After a moment or two the display shows that the calculated root C is 0.739085126 and pressing [EXE] once more, the number of iterations X used in the calculations was 18. The maximum error in the root is  ±0.0000060, smaller than the tolerance of ±0.00001 requested, as expected (the maximum error can be calculated from C-A). So we can probably quote the root as being 0.7391 to four decimal places and be confident that this is correct (other maths software shows that the root is 0.739085133 to 9 decimal places).

Repeating the bisection with different intervals but the same level of tolerance in the error gives the following results:-

A=0.5, B=1, C=0.739082336, X=16
A=0.25, B=1.5, C=0.739091873, X=17
A=0.7, B=0.8, C=0.73908081, X=14

and this confirms that quoting the root as 0.7391 to 4 d.p. is correct.

Root finding using the Bisection Method - the program

I discussed the logic for this program in a previous blog. Here is the logic converted into a program for the fx-50F:-

?→A:A→C:0→X:Goto 1:Lbl 2:Y→M:?→B:?→D:Lbl 4:X+1→X:.5(A+B)→C:.5(B-A)<D=>Goto 0:Goto 1:Lbl 3:Y=0=>Goto 0:If MY>0:Then C→A:Y→M:Else C→B:IfEnd:Goto 4:Lbl 0:C▲X▲Lbl 1:C-COS(C)→Y:X=0=>Goto 2:X≠0=>Goto 3:

This is pretty much how the program looks when keyed into the calculator in edit mode. Special Program Commands (? → => : Lbl = < > ≠ If Then Else IfEnd ▲ 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, D, X, Y, and M 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 132 steps in length.

Root finding using the Bisection Method - comments on the logic

In my previous blog entry I described one way in which to implement the bisection method. Most programmers would probably be aghast at my liberal use of the Goto command, but in this situation I think it is required.

There are necessarily two calculations of the function f(x) in different places of the logic; f(a) at the beginning of the program and f(c) during the iterations. As I only want to have the function in one place in the program (to make it easier to edit) I needed to have a subroutine, but the fx-50F doesn't have this functionality. To make matters worse, the control structures such as While/WhileEnd do not allow you to jump in and out of them using a Goto command and they also cannot be nested. So in the end I had to go back to basics and use Goto commands with a structure similar to a subroutine and a While/WhileEnd loop.

The 'subroutine' can be found at the end of the logic after the statement Display X (this is the actual end-point of the program). I have included another memory variable X which counts the number of iterations that were needed to obtain the root to the required accuracy. Conveniently, X has been used to identify where in the program the subroutine has been called from.

The function that appears in this subroutine is f(x)=x-cos(x) and I will say more about this in a later blog.

Root finding using the Bisection Method - the logic

We want to come up with an algorithm for finding the root of a function f(x) at x=α (where f(α)=0) using the bisetion method. We have 6 key pieces of information which we need to store in the calculator's memory as the iteration process described takes place. For the interval a<x<b being bisected, these are the value of a, b, δ, f(a), c=(a+b)/2, and f(c). Let the memory variables which store these values be A, B, D, M, C and Y, respectively. δ is the accuracy to which we wish to calculate α.

On the first iteration we key in the values a and b which span the root and store them in A and B. We also key in the value of δ and save it in D. We calculate c=(a+b)/2, f(a) and f(c) and store them in C, M and Y, respectively. Our estimate of the value of α in this first iteration is c and the error in the estimate is (b-a)/2. So if (b-a)/2 is less than our tolerance δ or if f(c)=0 then we are finished and we can display the contents of C as our estimate of α.

If, on the other hand, (b-a)/2 is greater than our tolerance δ and f(c) is not zero we can determine whether the root lies in the range a<x<c or c<x<b. If the product f(a)f(c)>0 then the root lies in the interval c<x<b otherwise it lies in the interval a<x<c. So if the product MY>0 we can store C in A and Y in M otherwise we store C in B. This gives us our new end points a and b (held in A and B) to process in the next iteration and, possibly, a new value of f(a) in M.

On the second and subsequent iterations we recalculate c=(a+b)/2 and f(c) using the values in A and B and store them in C and Y, respectively. We can then repeat the processing that we did in the first iteration.

The processing will look something like this:-

Input A
Store A in C
Store 0 in X
Goto Label 1
Label 2
Store Y in M
Input B
Input D
Label 4
Store X+1 in X
Store 0.5(A+B) in C
If 0.5(B-A) is less than D then Goto Label 0
Goto Label 1
Label 3
If Y=0 then Goto Label 0
If MY is greater than 0 then
   Store C in A
   Store Y in M
Else
   Store C in B
IfEnd
Goto Label 4
Label 0
Display C
Display X
Label 1
Store C-Cos(C) in Y
If X=0 then Goto Label 2
If X is not equal to 0 Goto Label 3