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

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

Multiplicative Inverses in Zn

In Cryptography it is often necessary to find the multiplicative inverse of a number in Z_n. Z_n is the set of integers {0,1,2,3,...,n-1} so, for example, Z₇ is the set {0,1,2,3,4,5,6}. If we are multiplying two elements in Z₇, say 3 and 5, then we can obtain a result in Z₇ if we use modular arithmetic i.e. 3x5=15=2x7+1, so 3x5=1 (mod 7).

If j and k are elements of Z_n, then k is the multiplicative inverse of j (and vice versa) if jxk=1 (mod n). So in the above example, 5 is the multiplicative inverse of 3 in Z₇. Note that j has a multiplicative inverse in Z_n if and only if j and n are coprime (i.e. j and n have no common factor apart from 1).

One way in which to calculate the multiplicative inverse of a number in Z_n is to use Euclid's Algorithm. Armed with a programmable calculator, however, the multiplicative inverse can be found by brute force.

Approximate Integration - execution of the program

This program will produce the definite integral of the function f(x)=x² between the limits x=A and x=B (B>A). On running the program, the display shows A? (the lower limit of integration). To begin with, enter 0 and then press [EXE]. The display then shows B? (the upper limit of integration). Enter 3 and press [EXE]. The display then shows D?. This is the number of pairs of strips that the interval between x=A and x=B is to be divided into. Enter 2 and press [EXE]. The display then shows M and 9.

This value is correct because the integral of f(x)=x² between the limits x=0 and x=a is (1/3)a³ and with a=3, this is (1/3)3³=9. This is a useful way to check that the code has been implemented correctly.

Let's try something harder. It is very difficult to obtain an expression for the integral of the function f(x)=exp(-x³) with respect to x but an estimate of the definite integral can be obtained numerically. Changing the function in the program from x² to exp(-x³) can be done by following the notes in this previous post. For the new function key in [SHIFT] [e^x] [(-)] [ALPHA] [X] [x³] [)].

Integrating this new function between the limits of x=0 and x=1 we obtain the following results as the number of pairs of strips is increased from 1 to 9:-

D=1 M=0.816311175
D=2 M=0.807843586
D=3 M=0.807572494
D=4 M=0.807530216
D=5 M=0.807518914
D=6 M=0.807514894
D=7 M=0.807513181
D=8 M=0.807512351
D=9 M=0.807511911

As can be seen, the results show that an estimate of this integral is 0.808 to 3 significant figures (as there is no change in the fifth significant figure for D>5) and this value agrees with other estimates.

Approximate Integration - further comments on the program

One of the other things to note about this program is that I have had to reuse several of the memory variables. The fx-50F comes with six memory variables (A, B, C, D, X and Y) and one independent memory (M) but this was not sufficient for the purposes of the program. For example, I calculate the value of the x-ordinate and store it in memory variable X, but then immediately calculate the y-ordinate, x², and store that in memory X too. This can be done because the x-ordinate is not needed after the y-ordinate has been calculated.

Approximate Integration - comments on the program

An important part of the program I presented previously is defining the function that is being integrated. I chose the simple function f(x)=x² and the statement in the code :X²→X: achieves this. In order to integrate a different funtion the program must be edited. For example, if I wanted to have f(x)=cos(x), then I would replace :X²→X: with :cos(X)→X:.

To edit a program it is best to ensure that your calculator is in insert mode by keying [SHIFT] [INS]. This will insert text before the vertical cursor and not overwrite existing text. To edit the program key [MODE] [6] [1] and then your program number. Use the replay button to move the blinking vertical cursor | to after X², so that it looks like :X²|→X:. Then press [DEL] twice and replace the function by keying [COS] [AlPHA] [X] [)], say.

It is a bit tedious, but I haven't found a way round this on the fx-50F.

Approximate Integration - the program

The logic for the approximate integration program was presented previously. Here is the logic converted into a program for the fx-50F:-

?→A:?→B:?→D:0→M:(B-A)÷(2D)→Y:1→C:While C≤D:1→B:Lbl 0:
B=1=>A+2(C-1)Y→X:B=2=>A+2CY→X:B=3=>A+(2C-1)Y→X:X²→X:
YX÷3→X:B=3=>4X→X:XM+:B+1→B:B≤3=>Goto 0:C+1→C:WhileEnd:
M▲

This is pretty much how the program looks when keyed into the calculator in edit mode except it all appears on one continuous line. 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, D, M, X and Y 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 136 steps in length.

Approximate Integration - the logic

I previously gave a formula (Simpson's Rule) for finding the approximate integral of a function f(x) between x=a and x=b (b>a). This consisted of a sum of n bracketed terms corresponding to the approximate area under the curve of n pairs of strips, with each strip being of width h.

This sum provides the logic for calculating the area under the curve and hence a value for the definite integral as follows:-

Input A (the lower limit of the integration)
Input B (the upper limit of the integration)
Input D (the number of pairs of strips)
Store 0 in M (this will hold the sum of the area)
Store (B-A)/2D in Y (the value of h, the width of a strip)
Store 1 in C (C indicates which pair of strips is being processed)
While C is less than or equal to D
   Store 1 in B (B now indicates one of three ordinates)
   Label 0
   If B=1 then
      Store A+2(C-1)Y in X (left ordinate)
   EndIf
   If B=2 then
      Store A+2CY in X (right ordinate)
   EndIf
   If B=3 then
      Store A+(2C-1)Y in X (middle ordinate)
   EndIf
   Store f(X) in X
   Store YX/3 in X
   If B=3 then
      Store 4X in X
   EndIf
   Store M+X in M (the summation)
   Store B+1 in B
   If B is less than or equal to 3 then
      Goto Label 0
   EndIf
   Store C+1 in C
WhileEnd
Display M

Approximate Integration

Sometimes it is necessary to obtain a numerical estimate of the definite integral of a function f(x) between x=a and x=b (where b>a). This integral is represented by the area bounded by the curve y=f(x), the x-axis and the ordinates at x=a and x=b. One reasonably accurate method of estimating this area, and hence this integral, is to use Simpson's Rule. Here the area is approximated by

(h/3)(y₀ + 4y₁ + y₂) + (h/3)(y₂ + 4y₃ + y₄) + ... + (h/3)(y_2n-2 + 4y_2n-1 + y_2n)

The area is divided into 2n strips of width h=(b-a)/2n and the ordinates are given by y₀=f(a), y₁=f(a+h), y₂=f(a+2h) etc. The method uses a quadratic of the form AX²+BX+C to approximate the function between three adjacent points on the curve (i.e. between two adjacent strips).

Prime Factors - increasing efficiency

Previously I mentioned that for simplicity I have just incremented the prime factor variable X in my program by one each time, so when X is 3, the next prime factor to be tried is 4. As I said then, clearly 4 is not a prime number and so it is ineffecient to be using this number in the search for prime factors (the same is also true of 6, 8, 9, 10...etc.).

So, how could this be improved? Well, the answer is to add some logic like this after incrementing X:-

If X=4 then
   X=X+1
If-End

So this will mean that 4 is skipped and the next prime number that is tried is 5. A similar approach can be used when X=6 (the next prime being 7), but at X=8 we need to add 3 as the next prime number is 11.

I will leave it to the reader to make these adjustments to their code.

Prime Factors - execution of the program

A program for finding the prime factors of a number on the fx-50F was described previously. Once it has been entered into the calculator you can press [Prog] [N] to run the program (N is the number of the program area where the code is stored). A? is displayed on the screen, prompting you to enter a number to be factorised.

Try the number 63 which we looked at before. So enter 63 and press [EXE]. The calculator displays Then X and 3. So, 3 is the first prime factor. Press [EXE] again and the calculator displays Then X and 3 again. So 3 is the second prime factor. Finally, press [EXE] once more and the display shows Then A and 7. The display of Then A rather than Then X indicates that this is the last prime factor to be found. So 63=3x3x7 which is correct. Press [AC] to exit the program.

If you try the number 61, the calculator displays Then A and 61 indicating that 61 is only divisible by itself and 1, so 61 is a prime number and has no other factors.

It is inadvisable to use this program for numbers much bigger than 500 as the processing may take a while to complete.

Prime Factors - the program

The logic for the prime factors program was presented previously. Here is the logic converted into a program for the fx-50F:-

?→A:2→X:Lbl 1:X→B:While B<A:B+X→B:WhileEnd:If B=A:Then X▲
A÷X→A:Goto 1:IfEnd:X+1→X:If X²>A:Then A▲IfEnd:Goto 1

This is pretty much how the program looks when keyed into the calculator in edit mode except it all appears on one continuous line. Special Program Commands (? → : Lbl While WhileEnd If Then 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 and X are input using [ALPHA] [A], [B] or [X] keys. The ":" Separator Code command is used frequently and can be alternatively input using the [EXE] key.

This program is 66 steps in length.

Prime Factors - Comment on the logic

There are a few things that need to be mentioned about the logic of my prime factors program.

The first is that for simplicity I have just incremented the prime factor variable X by one each time, so when X is 3, the next prime factor to be tried is 4. Clearly, 4 is not a prime number and so it is ineffecient to be using this number (the same is also true of 6, 8, 9, 10...etc.). However, it does not cause a problem with the logic, it just means that the program takes longer to run. I will address this issue again later.

Secondly, I have stopped the program when the square of the prime factor variable X is larger than the number A that we are testing. Why is this? Well, suppose our number is 63 and we have found the prime factors, 3, 3 and 7. We might try and see if the next prime number 11 is a factor of 63, but 63 divided by 11 is about 5.7. So we would already have found the factor 5 (if it was a factor) when looking for factors between 3 and 7. So in our test to see where to stop searching for prime factors we note that 7 squared is 49 (which is less than 63) and 11 squared is 121 (which is greater than 63), so it is clear that we can stop our prime factor search at 7.

Finally, we note that it would be much more efficient to use division than multiples of prime factors. For example, 63 divided by 2 gives 31 remainder 1. 63 divided by 3 gives 21 remainder zero. So, if we had knew the remainder, we could easily determine whether 2, 3, ... etc were factors or not. Unfortunately, the fx-50F is missing this functionality and so we have had to restort to using multiples of prime factors.

Prime Factors - the logic

How could we go about determining the prime factors of a composite number like 63 using the fx-50F? One way to do this is as follows. We begin with the smallest possible prime factor 2 and see if 63 is a multiple of 2. Multiples of 2 are 4, 6, 8, 10, . . . etc, which is just 2+2, 4+2, 6+2, 8+2, . . . etc. For each multiple, we can check to see if it is less than or equal to 63. If (a) it is less than 63, then we need to check the next multiple in the sequence. If (b) it is equal to 63, then 63 is divisible by 2. If (c) it is greater than 63, then 63 is not divisible by 2.

If (b) 2 is or (c) 2 is not a prime factor of 63 we need to continue looking for prime factors. If 2 is a prime factor of 63 we can divide 63 by 2 and obtain a new number with which to start the process all over again. If 2 is not a prime factor of 63 we can move to the next biggest prime factor 3 and start comparing 63 to multiples of 3. This process is repeated until we have found all the prime factors of 63.

The logic for a program to do this would be:-

Input A (the number whose prime factors are sought)
Store 2 in X (the first prime factor to try is 2)
Label 1
Store X in B (B holds the multiple being tested)
While B less than A
   Store B+X in B
While-End
If B=A then
   Display X (the prime factor found)
   Store A/X in A (the next number to be tested)
   Go to Label 1
If-End
Store X+1 in X (try the next prime factor)
If X squared is greater than A then
   Display A (the last prime factor to be found)
If-End
Go to Label 1

Prime Factors

Let's now consider something more interesting. In mathematics we often want to know whether a number like 63 is divisible by anything other 1 and 63. We know that in this case 63 is divisible by 3, 7, 9, and 21. But what about a number like 61? Well, no whole number will divide into 61 other than 1 and 61, so 61 is what is known as a prime number.

So consider the number 63 again. Using our known factors of 3, 7, 9 and 21 we find that 63=3x21 or 63=7x9 (there is no other combination of two factors). But 9=3x3 and 21=3x7, so in fact 63=3x3x7. This last representation is irreducible, as 3 and 7 are prime numbers themselves and so 3x3x7 represents the prime factors of 63. In fact, composite numbers like 63 are all constructed from a unique set of prime factors.