(Comp.sys.handhelds)
Item: 1892 by jurjen at cwi.nl
Author: [Jurjen NE Bos]
  Subj: HP48: The ultimate factoring program (up to 18 digits)
  Date: Fri Feb 01 1991 16:22 

This is a program that can factor large integers (up to 63 bits) on the
HP48.  It is inspired by the postings of Klaus Kalb of last week.  In fact, I
used two of Klaus' routines (see below).
The program has the following features:
- very fast (average 4 seconds for numbers of 12 digits,
  1 minute for 18 digits.)
- uses two different algorithms for factoring (trial division and Pollard rho)
  and proves all factors prime (using Selfridge's theorem).
- combines RPL (the language in the manual) with forth (what's in the ROM) and
  machine code.  This gives speed and servicability at the same time.
- The recursive calls are nicely shown during the computation.
- A separate prime testing routine.

There's also a bug:
- It doesn't quite work with wordsizes that are a multiple of four.  This is
  due to the carries.  If you set the wordsize to 64 bits, it will be
  modified to 63 to make it work.

[Note: Jurgen fixed this bug; see his posting at the end.  I put his bugfix
into his program; this debugged version is what's in FACTOR on the disk. -jkh-]

How to use:
Download the object at the end of this posting.  Use ASC-> to convert it into
a directory.  Save it under the name FACTOR.  Enter the directory.  Press
DEMO.  See what happens.
Remark: factoring numbers with large factors (more than 7 digits) can take
up to ten minutes.
Mind the bug above: if you set the wordsize to a multiple of four, factoring
might take forever.  [No, it won't; see bugfix note above.  -jkh-]

Some interesting numbers to try:
3^37+1 (to get this, type 63 STWS #3 #37 #1 NEG POWMOD 1 +).  The program
  thinks it found a large prime factor, but quickly finds out it isn't.
100264053529.  This number looks even more prime.

(Intermezzo:)
Some puzzles for experienced HP48 users:
        (For information, send Email.  Please realize these are PUZZLES, so
        they aren't trivial, and explanation is complicated.)
- Compute 'SIN(180_\^o)' .  You do not get 0.  Why?
- The time to compute '253!' is much longer than the time to compute '253.1!'.
- To compute the factorial for large negative x, use
'\pi*x/(SIN(\pi*x)*(-x)!)'
  instead of 'x!' for faster results.  (Don't forget RAD.)

Happy Hacking,
        Jurjen Bos

-------------------Listing of the FACTOR directory (do not download)------
For the nosy people among you, here the listing of the routines:
[FYI; download FACTOR or FACTOR.ASC from disk.  -jkh-]

DEMO    @ Factor a random integer and time the computation
\<< DEC TICKS RAND 1000000 * R\->B RAND 1.E12 * + RAND 1.E18 * + FACTOR
  TICKS ROT - B\->R 1.220703125E-4_s * SWAP
\>>

PRIME?  @Prime test     [NOTE: See bugfix posting below.  -jkh-]
\<< # 0d +
  CASE DUP # 2d <
    THEN DROP 0
    END DUP # 210d BGCD # 1d ==
    THEN RCWS 63 MIN STWS 0 'L' STO LPRT 'L' PURGE
    END DUP # 121d <
    THEN B\->R { 2 3 5 7 } SWAP POS 0 \=/
    END DROP 0
  END
\>>

FACTOR  @The factoring routine  [NOTE: See bugfix posting below.  -jkh-]
\<< RCWS 63 MIN STWS # 0d + 0 'L' STO LFAC 'L' PURGE
\>>

MULMOD  @Modulo multiplication of binaries
(In machine code.
a b n -> a*b mod n
)

POWMOD  @Modulo powering of binaries
(In machine code.
a e n -> a^e mod n
)

BGCD    @Greatest commond divisor of binaries
(In machine code.  Written by Klaus Kalb.
)

TRIAL   @Trial division routine
(In machine code and forth.  Machine code part written by Klaus Kalb.
)

MILRAB  @Miller-Rabin test
(In forth.
n base -> 1 if n is pseudoprime 
)

LFAC    @Internal factoring routine
\<< DUP 'L' INCR DISP "Trial division" 'L' INCR DISP # 2000d TRIAL
  IF DUP # 1d ==
  THEN DROP { }
  ELSE 1 \->LIST        @List of composite factors
  END
  WHILE DUP SIZE
  REPEAT SWAP OVER 1 GET
    CASE DUP # 4012009d <
      THEN B\->R +
      END DUP L 1 - DISP "Primetest" L DISP DUP LPRT
      THEN
        IF DUP # 1000000000000d <
        THEN B\->R
        END +
      END "Pollard \Gr" L DISP
      WHILE P\Gr DUP # 1d ==
      REPEAT DROP
      END 2 \->LIST ROT SWAP + SWAP
    END SWAP 2 OVER SIZE SUB
  END DROP 'L' "
" OVER DECR DISP 1 STO-
\>>

LPRT    @Internal primtesting routine
\<<
  CASE DUP # 3d MILRAB NOT
    THEN 0
    END DUP # 25000000000d >
    THEN SRT DUP
    END DUP # 2d MILRAB NOT
    THEN 0
    END DUP # 1373653d <
    THEN 1
    END DUP # 5d MILRAB NOT
    THEN 0
    END DUP # 25326001d <
    THEN 1
    END DUP # 7d MILRAB NOT
    THEN 0
    END DUP # 3215031751d \=/
  END SWAP DROP
\>>

SRT     @Test primality of large numbers with Selfridge test.
\<< DUP 1 - LFAC # 3d \-> n l a
  \<< "Selfridge" L DISP 9 CF 1 1 l SIZE
    FOR k
      IF l k GET SWAP OVER \=/ 8 FC?C OR
      THEN
        CASE a n 3 PICK / n POWMOD # 1d DUP2 \=/
          THEN SWAP 3 PICK # 0d + n POWMOD \=/ 8 + SF
          END + 'a' STO+ RAND .1 \>=
          THEN
          END n a MILRAB
          THEN
          END 9 SF
        END
      END 8 FS? 9 FS? -
    STEP DROP 9 FC?C
  \>>
\>>

P\Gr    @Pollard rho factoring algorithm.
\<< \-> n
  \<< RAND n B\->R * R\->B DUP NEWOB
    WHILE n Code n BGCD DUP # 1d ==
    @ Code is a machine code routine that does 50 rounds
    REPEAT DROP
    END ROT ROT DROP2 n
  \>> OVER /
\>>
--------------------------------- end of FACTOR directory --------------------

(Comp.sys.handhelds)
Item: 1893 by jurjen at cwi.nl
Author: [Jurjen NE Bos]
  Subj: HP48: factoring program, patch
  Date: Fri Feb 01 1991 16:22 

The bug mentioned in the previous posting can be succesfully solved by the
following minor patches.  This is typed by hand.  The easiest way to patch is
just using VISIT instead of trying to download it.

Happy hacking,
        Jurjen Bos

[NOTE: The following bugfix has already been incorporated into the version of
FACTOR that's on this disk.  It is listed here FYI only.  -jkh-]

PRIME?
\<< # 0d + # FFFFFFFFFFFFFFFFh SR AND
  CASE DUP # 2d <
    THEN DROP 0
    END DUP # 210d BGCD # 1d ==
    THEN 0 'L' STO LPRT 'L' PURGE
    END DUP # 121d <
    THEN B\->R { 2 3 5 7 } SWAP POS  0 \=/
    END DROP 0
  END
\>>

FACTOR
\<< # 0d + # FFFFFFFFFFFFFFFFh SR AND 0 'L' STO LFAC 'L' PURGE
\>>