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.              )

%DUP4 { OVR2 OVR2 }
%POP4 { POP2 POP2 }

|0100
  #ffff #fffb                 ( n** ; number to check )
  DUP4 is-prime32             ( n** is-prime^ ; test for primality )
  STH emit/long #2018 DEO     ( [is-prime^] ; emit n )
  STHr emit/byte #0a18 DEO    ( ; emit boolean )
  #ff0f DEO BRK               ( ; exit )

( 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^ )
  DUP4 ,&x1 STR2 ,&x0 STR2                    ( x** ; store x )
  DUP4 #0000 #0002 u32-ne ?{ POP4 #01 JMP2r }   ( x** ; return if 2 )
  DUP #01 AND ?{ POP4 #00 JMP2r }             ( x** ; return if even )
  DUP4 #0000 #0003 u32-ne ?{ POP4 #01 JMP2r }   ( x** ; return if 3 )
  #0002 ,&inc STR2                            ( x** ; inc<-2 )
  #0000 #0005 DUP4 ,&i1 STR2 ,&i0 STR2        ( x** i** ; i<-5 )
  &loop                                       ( x** i** )
    LIT2 [ &x0 $2 ] LIT2 [ &x1 $2 ]           ( x** i** x** )
    LIT2 [ &i0 $2 ] LIT2 [ &i1 $2 ]           ( x** i** x** i** )
    DUP4 u32-mul u32-lt                           ( x** i** x<ii^ )
    STH DUP2 #ffff EQU2 STHr ORA              ( x** i** x<ii|i=0xffff^ )
    ?{                                        ( x** i** )
      ,&x0 LDR2 ,&x1 LDR2 ,&i0 LDR2 ,&i1 LDR2 ( x** i** x** i** )
      u32-mod u32-is-zero                         ( x** i** x//i^ )
      STH #0000 ,&inc LDR2 u32-add              ( x** i+2** [x//i^] )
      LIT2 [ &inc $2 ] #0006 EOR2 ,&inc STR2  ( x** i+2** [x//i^] ; inc<-inc^6 )
      DUP4 ,&i1 STR2 ,&i0 STR2 STHr           ( x** i+2** x//i^ ; i<-i+2 )
      ?{ !&loop }                             ( x** i+2** ; if x<j*j, loop )
      POP4 POP4 #00 JMP2r                     ( 0^ ; since i divides x, not prime )
    } POP4 POP4 #01 JMP2r                     ( 1^ ; didn't find divisors, prime )

@emit
  &long SWP2 emit/short
  &short SWP emit/byte
  &byte DUP #04 SFT emit/char
  &char #0f AND DUP #09 GTH #27 MUL ADD #30 ADD #18 DEO JMP2r