Fundamental Algorithms

The Number of Combinations

The number of combinations nCr can be obtained by dividing the number of permutations nPr by r!, but if we calculate it as mentioned now, the process of calculating nPr happens to cause an overflow, while the true result nCr is representable in BASIC.
Using a formula nCr = nCr-1・(n-r+1)/r, we can avoid intermediate results unnecessarily increasing in their size.

10 INPUT n,r
20 LET C=1
30 FOR i=1 TO r
40    LET C=C*(n-i+1)/i
50 NEXT i
60 PRINT C
70 END

Binary Distribution

We want to calculate nCx px(1-p)n-x.
When n is a large number, nCx happens to exceeds the largest manipulable number, or px(1-p)n-x happens to fall below the smallest manipulable positive number.
To avoid those, we use logarithms.

10 INPUT p,n
20 INPUT x
30 LET u=x*LOG10(p)+(n-x)*LOG10(1-p)
40 FOR i=1 TO x
50    LET u=u+LOG10(n-i+1)-LOG10(i)
60 NEXT i
70 PRINT 10^u
80 END

Prime Factorization

We want to factor a natural number n into prime numbers.
We proceed to divide n by every prime number while n is divisible by it starting from the smallest prime number.
First, letting f = 2, we proceed to divide n by f while n is divisible by f (lines from 130 to 160). Next, we have to assign the next prime number to f, but there are no problem when we let f be a number added one. The repetition in lines from 130 to 160 is not performed if f is not a prime number, since n has been divided by all factors of f. This continues until f becomes a divisor of n.
This repetition ends when all the factors of n has been output and n becomes 1.

100 INPUT n
110 LET f=2
120 DO UNTIL n=1
130    DO WHILE MOD(n,f)=0
140       PRINT f;
150       LET n=n/f
160    LOOP
170    LET f=f+1
180 LOOP
190 END


When f becomes greater than √n,   f cannot be a factor. Therefore, if we revise the termination condition, we can ends the repetition of increasing f earlier.

100 INPUT n
110 LET f=2
120 DO UNTIL f>SQR(n)
130    DO WHILE MOD(n,f)=0
140       PRINT f;
150       LET n=n/f
160    LOOP
170    LET f=f+1
180 LOOP
190 IF n>1 THEN PRINT n
200 END

Note.
Use INTSQR(n) instead of SQR(n) on the rational operation mode.

As only even prime number is 2, if we particularly handle with 2, we can only consider the odd factors afterward.
Furthermore, as it is wasteful to calculate SQR(n) for every time, we calculate SQR(n) only when the value of n has changed.

100 INPUT n
110 DO WHILE MOD(n,2)=0
120    PRINT 2;
130    LET n=n/2
140 LOOP
150 LET f=3
160 LET s=SQR(n)
170 DO UNTIL f>s
180    DO WHILE MOD(n,f)=0
190       PRINT f;
200       LET n=n/f
210       LET s=SQR(n)
220    LOOP
230    LET f=f+2
240 LOOP
250 IF n>1 THEN PRINT n
260 END


Notice.
The above program cannot be rewritten as follows because the limit of a FOR statement is evaluated only when the control enters the FOR ~ NEXT loop.

90  REM wrong program
100 INPUT n
110 DO WHILE MOD(n,2)=0
120    PRINT 2;
130    LET n=n/2
140 LOOP
150 LET s=SQR(n)
160 FOR f=3 TO s STEP 2
170    DO WHILE MOD(n,f)=0
180       PRINT f;
190       LET n=n/f
200       LET s=SQR(n)
210    LOOP
220 NEXT f
230 IF n>1 THEN PRINT n
240 END


Decide whether a number is prime or not.

A program that decide whether a number greater than 2 is prime or not.

100 INPUT n
110 LET s$="prime"
120 FOR f=2 TO SQR(n)
130    IF MOD(n,f)=0 THEN
140       LET s$="composite"
150       EXIT FOR
160    END  IF
170 NEXT f
180 PRINT n;"is ";s$
190 END

Note. If 1 is entered, the result is wrong. ( 1 is neither a prime number nor a composite number.)

As even numbers greater than 2 is all composite and the factors of an odd number are all odd, the following program decreases the time to decide.

100 INPUT n
110 IF n=2 then
120    LET s$="prime"
130 ELSEIF MOD(n,2)=0 THEN 
140    LET s$="composite"
150 ELSEIF n>1 THEN
160    LET s$="prime"
170    FOR f=3 TO SQR(n) STEP 2
180       IF MOD(n,f)=0 THEN
190          LET s$="composite"
200          EXIT FOR
210       END  IF
220    NEXT f
230 END IF
240 PRINT n;"is ";s$ 
250 END

Note. Use INTSQR(n) instead of SQR(n) in the rational operation mode.

Greatest Common Divisors

Euclidean algorithm, which find the greatest common divisor, can be written using recursion.

100 FUNCTION GCD(a,b)
110    IF b=0 THEN
120       LET GCD=a 
130    ELSE 
140       LET GCD=GCD(b, MOD(a,b))
150    END IF
160 END FUNCTION
170 INPUT a,b
180 PRINT GCD(a,b)
190 END


This recurrence can be rewritten without recursive call.
In the above program, line 140 is substantially equivalent to the following assignment.
(a, b) ← (b, MOD(a,b))
But this assignment cannot be written as follows, because the value of a has changed when the second LET statement is executed.
LET a=b
LET b=MOD(a,b)
As the same problem remains even if we reverse the order of LET statements, we prepare another variable t and substitute as follows.

LET t=b
LET b=MOD(a,b)
LET a=t

When these are repeated until b becomes 0, a becomes the greatest common divisor.

10 INPUT a,b
20 DO UNTIL b=0
30    LET t=b
40    LET b=MOD(a,b)
50    LET a=t
60 LOOP
70 PRINT a
80 END


This algorithm is often written putting the calculation of the remainder former, but this is equivalent to the above substantially.

10 INPUT a,b
20 DO UNTIL b=0
30    LET r=MOD(a,b)
40    LET a=b
50    LET b=r
60 LOOP
70 PRINT a
80 END

Fibonacci numbers

Fibonacci sequence {fn} is a sequence defined by the following recurrence.
f1 = 1, f2 = 1, fn = fn-1+fn-2


This algorithm can be written as follows using recursion.

10 FUNCTION f(n)
20    IF n=1 OR n=2 THEN LET f=1 ELSE LET f=f(n-1)+f(n-2)
30 END FUNCTION
40 INPUT n
50 PRINT f(n)
60 END

But this program is inefficient, since f(n-1) and f(n-2) are calculated separately.

Using an array.

DIM f(1000)
LET f(1)=1
LET f(2)=1
INPUT n
FOR i=3 TO n
   LET f(i)=f(i-1)+f(i-2)
NEXT i
PRINT f(n)
END

The above program saves waste of calculation, but wastes memories since the calculation uses only two former values.

Hence we rewrite as follows letting the two former values be d and e.

100 LET d=1
110 LET e=1
120 INPUT n
130 FOR i=3 TO n
140    LET f=e+d
150    LET d=e
160    LET e=f
170 NEXT i
180 PRINT f
190 END

Lines 150 and 160 mean renewal for the next repetition.

Binary Notation

Converting a natural number into binary notation
When a number is divided by 2, the remainder becomes the ones place digit. When the quotient is divided by 2, the remainder becomes the twos place digit. And when the quotient is divided by 2, the remainder becomes the 22s place digit.
That is to say, a repetition of dividing a number by 2, outputting the quotient and letting the quotient to be a new number converts a number to binary notation.
the following program outputs the digits from the ones place, that is, in the reverse order of usual notation.

10 INPUT n
20 DO UNTIL n=0
30    LET r=MOD(n,2)
40    PRINT r;
50    LET n=(n-r)/2
60 LOOP
70 END


If you want to display as the most significant digit comes to the left, use string manipulation.

10 LET s$=""
20 INPUT n
30 DO UNTIL n=0
40    LET r=MOD(n,2)
50    LET s$=STR$(r)&s$
60    LET n=(n-r)/2
70 LOOP
80 PRINT s$
90 END



Converting a decimal fraction into binary notation

We want to convert a decimal fraction x with 0 ≤ x<1 to binary notation.
We proceed the reverse way of the above. That is, when the number is multiplied by 2, the integer part is the digit in the 0.1s place in binary notation. And the fractional part is multiplied by 2, then the integer part is the digit in the 0.01s place in binary notation.

10 OPTION ARITHMETIC DECIMAL
20 INPUT x
30 PRINT "0.";
40 DO UNTIL x=0
50    LET p=INT(2*x)
60    PRINT p;
70    LET x=2*x-p
80 LOOP
90 END

Full BASIC handles decimal fractions exactly. Thus the above is a correct program in Full BASIC.
Note that some decimal fraction cannot be represented in a finite binary fraction. For example, when 0.1 is entered in the above program, the program will not terminate for ever. (Select Break in the Run menu.)

Integer Power Algorithm

Power operations such as the indices are powers of 2 can be done in fewer multiplications as follows.
a2=a×a
a4=a2×a2
a8=a4×a4
a16=a8×a8

Furthermore, integer power operation can be done in few multiplications, even if the index is not a power of 2, for example, since 25 = 16 + 8+1,
a25=a16×a8×a

Since the operation that makes a natural number into a sum of powers of 2 is the same as that of converting to binary notation, an integer power operation an can be done rapidly.

100 INPUT a,n
110 LET p=1
120 LET b=a
130 DO UNTIL n=0
140    IF MOD(n,2)=1 THEN LET p=p*b
150    LET b=b*b
160    LET n=INT(n/2)
170 LOOP
180 PRINT p
190 END

Although the above program has no practically worthless because Full BASIC has power operation, this algorithm can be applied to power operation of matrices or power operation in modular arithmetic.


The following program calculates the nth power of a 2 by 2 matrix a.

100 DIM a(2,2),b(2,2),m(2,2)
110 MAT INPUT a
120 INPUT n
130 MAT m=IDN
140 MAT b=a
150 DO UNTIL n=0
160    IF MOD(n,2)=1 THEN MAT m=m*b
170    MAT b=b*b
180    LET n=INT(n/2)
190 LOOP
200 MAT PRINT m
210 END

The following program calculates an in arithmetic modulo k.

100 INPUT k
110 INPUT a,n
120 LET p=1
130 LET b=MOD(a,k)
140 DO UNTIL n=0
150    IF MOD(n,2)=1 THEN LET p=MOD(p*b,k)
160    LET b=MOD(b*b,k)
170    LET n=INT(n/2)
180 LOOP
190 PRINT p
200 END

Fermat Test

Using Fermat's little theorem
for an integer a with 0<a<k, ak-1 ≡ 1 (mod k) if k is a prime number,
we can conclude that k is not a prime number if the (k-1)-th power of a does not become 1 in arithmetic modulo k for some a with 1<a<k.
The following program runs the above test for a=2. Change the line 100 if you want to run for another a.

100 LET a=2
110 INPUT k
120 LET n=k-1
130 LET p=1
140 LET b=a
150 DO UNTIL n=0
160    IF MOD(n,2)=1 THEN LET p=MOD(p*b,k)
170    LET n=INT(n/2)
180    LET b=MOD(b*b,k)
190 LOOP
200 IF p<>1 AND k>a THEN PRINT k;"is not a prime number." 
210 END

An advantage of calculating in modular arithmetic is that numeric values only upto (k-1)2 have to be handled with accuracy.
The decimal mode of Decimal BASIC has a precision of over 16 digit in a numerical expression, the above program can decide numbers upto 8 digit. If the number is larger than that, use the rational operation mode.


Back