nxu/primes32.tal

( primes32.tal                                         )
(                                                      )
( Uses a simple trial-divion method to find primes.    )
(                                                      )
( To determine if x is prime we:                       )
(                                                      )
( 1. Check if x is 2 (prime)                           )
( 2. Check if x is even (not prime)                    )
( 3. Check if x is 3 (prime)                           )
( 4. Starting with i=5, we see if x%i is 0             )
(   a. We alternately increment i by 2 and 4           )
(   b. We stop when x < i*i or i=0xffff                )
( 5. If we didn't find an i, x is prime.               )
(                                                      )
( The reason we alternate our increment is because we  )
( know that x%6 must equal 1 or 5. if x%6 was 3 then x )
( would be divisible by 3.                             )
(                                                      )
( This method can be fast for some large composite     )
( numbers but is slower for large primes.              )
(                                                      )
( On my machine, checking 0x7fffffff took 0.5 seconds  )
( and checking 0xfffffffb took 0.9 seconds. Both are   )
( prime numbers.                                       )
(                                                      )
( Smaller primes also run fairly quickly: 0x17b5d was  )
( determined to be prime in 0.02 seconds.              )

%EMIT { #18 DEO }
%DIGIT { #00 SWP ;digits ADD2 LDA EMIT }
%SPACE { #20 EMIT }
%NEWLINE { #0a EMIT }
%EMIT-BYTE { DUP #04 SFT DIGIT #0f AND DIGIT }
%DEBUG { #ff #0e DEO }

|0100
    ( number to check comes first )
    #fffe #0001
    OVR2 OVR2 ;is-prime32 JSR2 ( test for primality )
    STH ;emit-long JSR2 SPACE STHr EMIT-BYTE NEWLINE ( output )
    #00 #00 DIV ( exit with /0 to make timing easier )
    BRK

( include 32-bit math library )
~math32.tal

( return 01 if x is a prime number, else 00 )
( works for x >= 2 )
@is-prime32 ( x** -> bool^ )
    ,&x1 STR2 ,&x0 STR2 ,&x0 LDR2 ,&x1 LDR2 ( store and reload x )
    DUP4 #0000 LIT2 &two 0002 ;ne32 JSR2 ( x is 2? )
    ,&not-two JCN POP4 #01 JMP2r ( 2 is prime )
    &not-two DUP #01 AND ( x x&1 )
    ,&not-even JCN POP4 #00 JMP2r ( x is even: not prime )
    &not-even DUP4 #0000 #0003 ;ne32 JSR2 ( x is 3? )
    ,&not-three JCN POP4 #01 JMP2r ( 3 is prime )
    &not-three #0000 ,&i0 STR2 #0005 ,&i1 STR2 ( x i<-5 )
    ,&two LDR2 ,&inc STR2
    ,&i0 LDR2 ,&i1 LDR2 ( load our candidate, i )
    ,&loop JMP ( jump over register data to loop label )
    [ &i0 0000 &i1 0000 &x0 0000 &x1 0000 &inc 0000 &mask 0006 ] ( registers )
    &loop ( x i )
    ,&x0 LDR2 ,&x1 LDR2 ,&i0 LDR2 ,&i1 LDR2 ( x i x i )
    DUP4 ;mul32 JSR2 ;lt32 JSR2 ( x i x<i*i )
    STH DUP2 #ffff EQU2 STHr ORA ( x i x<i*i||i=0xffff )
    ,&finished JCN ( x i )
    ,&x0 LDR2 ,&x1 LDR2 ,&i0 LDR2 ,&i1 LDR2 ( x i x i )
    ;mod32 JSR2 ;is-zero32 JSR2 ( x i x//i^ )
    STH #0000 ,&inc LDR2 ;add32 JSR2 ( x i+2 )
    ,&inc LDR2 ,&mask LDR2 EOR2 ,&inc STR2 ( inc<-inc^6 )
    ,&i1 STR2 ,&i0 STR2 ,&i0 LDR2 ,&i1 LDR2 ( write i+2 to register )
    STHr ( x i+2 x//i^ )
    ,&i-divides-x JCN ( x j )
    ,&loop JMP ( if x<j*j, loop )
    &i-divides-x POP4 POP4 #00 JMP2r ( since i divides x, not prime )
    &finished POP4 POP4 #01 JMP2r ( didn't find divisors, prime )

( print a long value as hex )
@emit-long ( x** -> )
    SWP2 SWP EMIT-BYTE EMIT-BYTE
         SWP EMIT-BYTE EMIT-BYTE JMP2r

( convenience for less branching when printing hex )
@digits
    30 31 32 33 34 35 36 37
    38 39 61 62 63 64 65 66