Prime factorization with HP-71 BASIC

One of the very first projects I started  after  getting  my  new
HP71  was  a prime factorization program.   I wanted a relatively
short and simple program, preferably  using  some  fairly  common
algorithm,  so  that  comparison  with other machines and earlier
programs would be possible, and such that it  could  be  used  to
produce  all  factors  of  a  number,  either interactively or by
another program (e.g.   for computing Euler's  Phi-function),  as
well as for just primality testing.

After  some  thought,  I decided I wanted a subprogram that would
return the factors one at a time, along with their exponents, and
a main program that would call it  repeatedly  until  the  entire
number was factored.

That  left little enough choice on the subprogram parameters: the
number being factored (N), factor found (F) and its exponent (E).

One problem was that the sub had to be able to return  a  factor,
and later continue from where it was. Simple solution to this was
using F as input also, telling what potential factors had already
been tried, and dividing found factors out from N.  Note that the
DIM statement in line 10 clears F.

Now I knew how I wanted subprogram FACTOR to behave.  The problem
remaining was how to implement that.

The  algorithm chosen is a simple one: try 2, 3 and 5 first, then
all not divisible with them up to the square root of n.  This has
been implemented on any number of computers and programmable cal-
culators, notably including the HP-41.

The central loop of my first effort looked something like this:

 190 S=SQR(N) @ F=5
 200 F=F+2 @ IF NOT RMD(N,F) THEN 300
 210 F=F+4 @ IF NOT RMD(N,F) THEN 300
 ... (6 more similar lines)
 280 IF F<S THEN 200
 290 F=N

Some   experimenting   showed  that  constants  are  faster  than
variables in the 71, like in the 75.   (More surprising  was  the
discovery  that  longer  constants are actually faster than short
ones, except one-digit integers which have tokens of  their  own.
This  was  later  verified  by  IDS  3: the beast always reads 12
nibbles and then clears what isn't needed.)

Then it occurred to me that it isn't really necessary to save the
factor candidate in F each time. It only needs to be done once in
each loop - and then FOR-NEXT can be used. So:

 190 FOR F=7 TO SQR(N) STEP 30
 200 IF NOT RMD(N,F) THEN 300
 210 IF NOT RMD(N,F+4) THEN F=F+4 @ GOTO 300
 ...
 280 NEXT F @ F=N

At this point I naturally started looking for other  things  that
might not be necessary each time.  The result you can see in FAC-
TOR: instead of a separate IF NOT for each remainder, bunch  them
together with AND and use one IF NOT for the lot.   Unfortunately
all 8 won't fit on a single line, but four at a time isn't bad.

Testing four divisors simultaneously causes the problem that then
we'll only know that one of the four is a factor, so we must test
them again separately. This is done in lines 220-280 - which look
rather like the original search loop.

One more adjustment in the loop was needed  to  allow  continuing
from the previous factor. It is neither possible nor necessary to
start exactly from it: instead, we must find a number of the form
30k+7 not greater than the next potential factor.  The expression
(F DIV 30)*30+7 is smaller than F  always  except  when  F=30k+1;
then it gives 30k+7, and there are obviously no primes in between
(when k>0).

Factors 2, 3 and 5 are handled on line 170.  The odd-looking exp-
ression gives 2 for all F<2, 3 for 2<=F<3 and 5 for 3<=F.    Then
the  new  F is tried, using ON GOTO because the 71 does not allow
another IF here.   The functions SGN and RMD don't seem to bother
ON GOTO.   The program will fail, however, if OPTION ROUND POS or
NEG is in effect.

Finally, line 160 is a safeguard against non-positive inputs.

How fast it is, then? Discovering 9999999967 is a prime  takes  6
min  40  s;  999999999989 takes one hour 8 min 40 s, in a 71 with
634 KHz clock.   That makes it about 15% faster than  Jim  Horn's
program  in  PPC  CoJ  V1N2P22 (which runs unmodified in the 71),
even though Jim used a theoretically faster  algorithm  (skipping
multiples  of  7  too),  and some 10 times faster than best HP-41
programs I've seen.

In case you have been wondering why I didn't  use  FORTH  in  the
first place, it being an inherently faster language, only because
I hadn't got my FORTH/Assembler ROM then.  As soon as I got it, I
started ... - but that's another story. Maybe next issue.

                          Tapani Tarvainen [337]
                          Kiljanderinkatu 2 A 7
                          SF-40600 JyvÌskylÌ
                          Finland