DELVING NUMERICALLY DEEPER (16K ZX81)
by Frank O'Hara,
from ZX Computing 2, August/September 1982, pages 28-31

Frank O'Hara from Surbiton in Surrey helped lan Logan decode the 8K ROM. Dr
O'Hara has continued his investigations into the operation of tne ZX81, and
here shares with us some of his discoveries, with notes on some programs on
elementary number theory for the ZX81.

Over the past year or so, Dr O'Hara has developed a few programs on elementary
number theory for the ZX81, having previously run a few such programs on a
programmable calculator, a Texas TI58, over about 2.5 years from mid-1978.

A couple of these programs are "one off", ie. ad hoc programs to solve a single
program. There is one of about 100 bytes of machine code which generates the
decimal representation of quite large powers of 2. He used this to obtain the
13395 digits of (2 to the 44497) minus 1, the largest known prime number,
discovered by Nelson and Slowinski using a Cray One computer in 1979. The
program took 2 hours 31 minutes to obtain this number on the ZX81.

Another even more exciting result was given by about 400 bytes of machine code,
including a multiple precision multiplication routine. This program actually
proved the primality of the first 15 Mersenne primes, up to and including (2 to
the 1279) minus 1, a number of 386 digits. It thus repeated some of the work of
"SWAC" in 1953, going far beyond what the desk calculators had done and
capturing some of the flavour of a historic moment (although 28 years later!).

Three other programs have a more general purpose flavour, I have called
them:

(a) SPRF: single precision prime factorization;
(b) MPRF: multiple precision prime factorization;
(c) FE24: Fermat's theorem used to test numbers up to 24 digits long for
          compositeness, ie. lack of primality.

The kernel of these 3 programs is the machine code multiple precision integer
division routine which finds the true integer quotient and true integer
remainder of an integer of arbitrary length with respect to another arbitrarily
long integer. Barden is quite mistaken when he describes this process as a
"cop-out" (How to Program Microcomputers, by William Barden, Jr. Sams,
Indianapolis, 1977, page 109). It is in fact very easy to program. It is a
simple extension of the standard restoring division of one or two bytes by one
byte. The shift is just a loop with its kernel as: LD A,(DE) : RLA : LD (DE),
A. The addition or subtraction is another loop centred on: LD A,(DE) : ADC A,
(HL) : LD (DE), A; and so on. In fact the only complexity arises when one has
to shorten the process in order to speed it up, as in the first of the 3
programs, SPRF. This program has a 4 byte dividend and 2 byte divisor, and
needs to use the exchange resisters H', L', D' and E' to gain speed.

The first program, SPRF, finds the smallest prime factor of any odd number from
5 to 4294967255 (2 to the 32 minus 1) in not more than 20 seconds. If the
number is prime, the program reports this. It goes about 1000 times as fast as
the Texas calculator did. This has been achieved by a series of improvements,
starting with a BASIC program that was only about 10 times as fast as the
calculator. The program contains about 400 bytes of machine code; 300 or so of
these are just a simple linear sieve designed to exclude multiples of 3, 5 and
7 as well as 2 and so gain a factor of 35/16 in speed. The speed has to be seen
to be believed. Numbers up to 7 digits long are dealt with instantaneously. The
largest 8, 9 and 10 digit primes in its range take 3, 10 and 20 seconds
respectively. The program can easily be adapted to print screenfuls of results
and so, for example, find the largest prime less than 2 to the 32 in one run.
By using random 8, 9 or 10 digit input one can use it to see primes probably
never seen before. (Only the first 10 or 11 million numbers are completely and
accurately listed as prime or composite.)

The second program, MPRF, is a general purpose prime factor finder and can deal
with numbers up to 77 digits long, as time permits. (The break key is active in
the machine code, so one can exit from it and look at the divisor.) Its speed
depends on the length of the number being examined. A 13 digit number is
analysed at about 10,000 divisors a minute; a 25 digit number is analysed at
about 3,000 divisors a minute; ie the divisor reaches the stated range when all
odd divisors are being tried. Dr. O'Hara used this program to check the factors
of the "repeated unit" numbers up to 30 ones (ie. one ninth of ((10 to the 30)
minus 1)) which had already been obtained, with much more labour, using a
calculator. Three of these results required the third program, too. Other uses
in conjunction with the third program are mentioned below. MPRF has 182 bytes
of machine code and quite a lot of BASIC to start and finish it.

The third program, FE24, tests a number N by raising 2 to the power N-1,
continually reducing modulo N. If the result is not equal to 1, then N is
composite. Otherwise N is called a "near-prime" or a "pseudo-prime". In fact N
is nearly always prime in this case, and the proof of its primality can be
completed provided N-1 can be completely factorized. This involves MPRF again.
Dr. O'Hara used FE24 in proving the primality of 19 ones, 23 ones and a 15
digit factor of 27 ones. He also used it to obtain results which he could not
get with the calculator. He found the largest N digit prime for N up to 13. He
has now extended that to 20, and is still looking at the next four, FE24 is
mostly in BASIC. To raise 2 to the power M, where M is usually N-1, it first
gets the binary decomposition of M. Then it repeatedly squares and, as
necessary, doubles the residue, starting from 1, and reducing each time mod N.
Only the reduction is in machine code, about 75 bytes. The multiple precision
multiplication is in BASIC. The program is fast enough for its purpose. It
tests a 24 digit number in under 20 minutes. Because of its powerful indirect
method it is thus more than a million times faster in achieving its aim than
MPRF on its own would be.

It would be nice to find a better technique for factorization than MPRF allows.
MPRF can factorize an arbitrary 13 digit number in not more than 5 hours. To
factorize larger numbers would be impractical. So far Dr. O'Hara has been lucky
in that the largest run needed was about 40 minutes for a factor of about
300,000. But it would be nice if a powerful technique like the use of quadratic
sieves, as described by Donald Knuth in "The Art of Computer Programming"
(Addison-Wesley 1969; vol. 2, pages 345-347) could be implemented on the ZX81.
Dr. O'Hara has studied this, but does not think it is feasible. The array
facilities available in Z80 machine code seem to be insufficient. In any case,
it may be that a 16 or 32 bit microprocessor would be needed.

A final note on what is perhaps his most spectacular result to date. The
beautiful and justly famous factorization of 17 ones into the product of the
two primes 2071723 and 5363222357 took under 20 minutes with MPRF, compared
with 20 hours on the calculator. Of course, there is a gain of a factor of 17
in time here, since one only needs to try every 34th divisor. So one reaches 2
million in the time it would normally take to reach 120,000.


The Assembler For SPRF

1. In order to allow plenty of room for BASIC, Dr. O'Hara started by putting
RAMTOP at 96,0 ie. at 24576d (24K; with 16K RAM attached). He then used
addresses 25471 to 25913 to hold the machine code (it has many subroutine calls
so is non-relocatable) and these addresses for other purposes:

26496-7: to save the contents of D'E' (not necessary, I later found).
26498-9: to save the contents of H'L' (essential to save H'L', but it could
         have been pushed on to the stack).
26510-1: to hold the divisor.
26512-5: to hold the number being tested.
26516-7: to hold the square root of the number being tested.
26518:   to hold a flag, 1 for a prime number, 0 for a composite number.

Once the number is entered (as a string) the BASIC stores it, sets the divisor
to 1 and enters the square-root ("bug-proof", because of the defective ROM). It
then calls the machine code and stays there until it is ready to announce
primality or print some factors. The latter will involve picking up the divisor
(factor) and the quotient (co-factor) and printing up to ten digits of the
quotient. So the repeated parts of testing for primality are all in machine
code for speed.

[List 600 of the program for instructions and examples values to try]

2. The Assembler. This contains 443 bytes. The first 9 just save H'L' and D'E'.
The next 351 are just a simple linear sieve which excludes multiples of 3, 5
and 7 as well as 2 and hence gains a factor of about 35/16 in speed. This part
is very repetitive, and only the beginning and end are shown below. Instead of
just adding 2 to the divisor (DE) it adds this sequence of numbers in a
perpetual loop (after trying the values 3, 5, 7 and 11):

2, 4, 2, 4, 6, 2, 6, 4, 2, 4, 6, 6, 2, 6, 4, 2, 6, 4, 6, 8, 4, 2, 4, 2, 4, 8,
6, 4, 6, 2, 4, 6, 2, 6, 6, 4, 2, 4, 6, 2, 6, 4, 2, 4, 2, 10, 2, 10.

Each time the divisor is set, the main division subroutine is called. This does
an ordinary restoring division, as explained in Zaks and Barden, with the
further refinement shown in the ROM division routine, of allowing a full 32
bits (16 here) in the divisor by saving any bit of the quotient which drops
into the carry. The remainder is tested for zero. If it is non-zero, the
divisor is tested against the square root of the number. If it is greater, the
program returns to BASIC to report a prime.

Step Label   Opcode           Comments
1            EXX
2            LD (6782),HL     Save H'L'.
3            LD (6780),DE     Save D'E'.
4            EXX
5            LD DE,0003       Set divisor to 3.
6            CALL 64E7,DIVN   Call main division subroutine.
7            LD DE,0005
8            CALL 64E7,DIVN
9            LD DE,0007
10           CALL 64E7,DIVN
11           LD DE,000B
12           CALL 64E7,DIVN
13     SIEVE INC DE           Add 2 to divisor by incrementing
14           INC DE           DE twice.
15           CALL 64E7,DIVN
16           INC DE           Add 4 to divisor by incrementing
17           INC DE           DE 4 times.
18           INC DE
19           INC DE
20           CALL 64E7,DIVN
21           INC DE
22           INC DE
23           CALL 64E7
24           INC DE
25           INC DE
26           INC DE
27           INC DE
28           CALL 64E7
29           LD HL,0006       Add 6 to divisor by using HL and
30           ADD HL,DE        the EX DE,HL instruction.
31           EX DE,HL         Eight and ten will be added in the
32           CALL 64E7        same way.

... sieve continues to step 205, ending with:

201          LD HL,000A       Add 10 to divisor by using HL and
202          ADD HL,DE        the EX DE,HL instruction.
203          EX DE,HL
204          CALL 64E7
205          JP 63AO,SIEVE    Loop for ever in the sieve.
206     DIVN EXX
207          LD HL,(6790)     N is loaded into D'E'H'L'. It would
208          LD DE,(6792)     be more efficient to do this at steps
209          EXX              5-8. (One day 1 must draw a flowchart, just
for fun!)
210          LD HL,0000       The remainder is set to zero in HL.
211          LD BC,20         Initialize count to 32 decimal.
212          AND A            Clear the carry Hag.
213     DIVL EXX              Enter the division loop.
214          ADC HL,HL        Shift the remainder dividend-
215          EX DE,HL         quotient left in HLD'E'H'L'.
216          ADC HL,HL
217          EX DE,HL
218          EXX
219          ADC HL,HL
220          JR C,650A,SAVE   If a bit drops into the carry, go and retrieve it for the
                              quotient.
221          SBC HL,DE        Trial subtract the divisor.
222          JRNC,650D,NRST   Go, if no carry, to no restore.
223          ADD HL,DE        Add back the divisor if there was carry.
224          AND A            Clear the carry and go with no bit
225          JR 650E,CONT     for the quotient.
226     SAVE AND A            Force no restore and one for the
227          SBC HL,DE        quotient here.
228     NRST SCF              Set the carry flag: one for the quotient.
229     CONT DJNZ 64F6,DIVL   Loop back for each bit of dividend (32 times).
230          EXX
231          ADC HL,HL        Move last bit into quotient.
232          EX DE,HL
233          ADC HL,HL
234          EX DE,HL
235          EXX
236          LD A,H           Now test the remainder.
237          OR L
238          JR Z,6527,FACT   Go if it is zero.
239          AND A            Clear the carry.
240          LD HL,(6794)     Put square root of N into HL,
241          SBC HL,DE        Subtract divisor.
242          RET NC           Return to sieve if more to do.
243          LD A,01          Otherwise, set flag for a prime
244          JR 652C,EXIT     and go to EXIT.
245     FACT XOR A            Reset flag for a factor.
246          LD (678E),DE     Save factor for BASIC.
247     EXIT LD (6796),A      Save flag for BASIC.
248          EXX
249          LD HL,(6782)     Restore H'L'.
250          LD DE,(6780)     Restore D'E'.
251          EXX
252          POP HL           Discard sieve return address.
253          RET              Return to BASIC.