\ prime.fs -- WIP
\ Miller-Rabin Primality Test
\ Version: 2026-06-19 Copyright Gan Uesli Starling
\ License: None. Free for all uses.

\ http://www.prime-numbers.org/   Download all prime number list up to 400 digits for $19.99

\ Dependent on words defined in files. These to be INCLUDED prior, outside of this file.
\ edit_me.fs
\ CWD.fs
\ defs.fs
\ dec_to_bin.fs
\ math.fs
\ bin_to_bcd.fs
\ bbs.fs

0 [IF] ( How to Proceed )

  1. Generate a big integer. If it's even, add 1. Call it n.
  2. Check if it's prime.
      a. Optional: find the GCD between n and the product of a moderate number of small primes > 2. If the GCD isn't 1, go to #1.
      b. Do Miller-Rabin primality test (https://en.wikipedia.org/wiki/Miller%E2%80%93Rabin_primality_test), if composite, go to #1.
      c. Do Baillie-PSW primality test. (https://en.wikipedia.org/wiki/Baillie%E2%80%93PSW_primality_test), if composite, go to #1.
  3. n is almost certainly prime. 

[THEN]

0 VALUE pri_0 \ Allocated array for large prime

\ Allocate or erase
: pri_0.anew? pri_0 0= IF a.new TO pri_0 ELSE pri_0 a.erase THEN ; 

\ Free if allocated
: pri_0.free? pri_0 0<> IF pri_0 a.free 0 TO pri_0 THEN ;

VARIABLE pri_bits
VARIABLE pri_cells
 
\ Align requested bits to stack width
: pri.bits.align ( u -- u ) \ Assume u = 511 on 64-bit Forth
  1 CELLS 8 *            ( 511 64 )
  2DUP MOD               ( 511 64 63 )
  DUP 0= 
  IF 2DROP ELSE - + THEN ( 512 )
;

\ Enforce limit: 2 cells <= bits <= half availble
: pri.bits.limit ( u -- u )
  DUP >R
  2 CELLS 8 * MAX
  A_LGTH_BITS 2 / MIN
  DUP R> <> IF
    CR ." Oops! Limit enforced! Bits = " DUP . 
    ." Min = 2 cells, max = A_LGTH_BITS / 2 " CR
  THEN
;

: rand.seed.churn
  CR ." To init rand_seed, enter a passphrase. "
  0 kybd.chances.inf
  rand_seed @
  pass_buff COUNT NIP ( cnt )
  0 DO
    pass_buff 1+ I + C@   \ Get one char
    I 8 CELLS MOD LSHIFT  \ Offset it to left
    XOR                   \ XOR with prior
  LOOP
  rand_seed !
;

: ask.prime.bits ( -- u )      \ Expect un-related 'u' on stack as input to "find.primes"
  2 CELLS 8 *                  \ Two cells is min bits
  A_LGTH 4 * >us               \ Half of available is max bits
  CR CR ." How many bits preferred (" 
  DUP . ." to "
  us@ . 8 EMIT ." )? "         \ Show max allowed bits
  CR ." > "                    \ Prompt on new line
  PAD 5 BLANK PAD 4 ACCEPT >R  \ Get user input
  PAD R> string.to.integer     \ ASCII to integer
  DUP us> >                    \ Above max...
  ROT 2 PICK > OR              \ ...or below min?
  IF DROP RECURSE THEN         \ If so, ask again
  pri.bits.align               \ Align to stack width
  pri.bits.limit               \ Keep within bounds, post-alignment
  DUP pri_bits !               \ Keep 
  8 CELLS / pri_cells !        \ Store cell-width of desired primes
;

\ Random with non-zero HSB
: hsb.nonzero ( -- u )
  BEGIN 
    random
    DUP MASK_HSB AND 0=
  WHILE
    DROP
  REPEAT
;
 
\ Fill array with a candidate for prime
: pri.candidate ( c-addr -- )
  DUP a.erase >R
  pri_cells @
  DUP >R
  0 DO               \ Stack up N cells to load
    I 0= IF
      hsb.nonzero    \ Insure HSB non-zero
    ELSE
      random
    THEN 
  LOOP               \ N random cells on stack
  1 OR               \ LSB must be odd
  R> R@ a.u!         \ Store into array
  R> DUP a.pos!      \ Ensure is positive
;

a.new CONSTANT a_tmp

\ SMALL PRIMES 

\ Primes < HEX 4FF

2 3 5 7 11 13 17 19 23 29 31 37 41 43 47 53 59 61 67 71 73 79 83 89 97 
101 103 107 109 113 127 131 137 139 149 151 157 163 167 173 179 181 191 
193 197 199 211 223 227 229 233 239 241 251 257 263 269 271 277 281 283 
293 307 311 313 317 331 337 347 349 353 359 367 373 379 383 389 397 401 
409 419 421 431 433 439 443 449 457 461 463 467 479 487 491 499 503 509
521 523 541 547 557 563 569 571 577 587 593 599 601 607 613 617 619 631 
641 643 647 653 659 661 673 677 683 691 701 709 719 727 733 739 743 751 
757 761 769 773 787 797 809 811 821 823 827 829 839 853 857 859 863 877 
881 883 887 907 911 919 929 937 941 947 953 967 971 977 983 991 997 
1009 1013 1019 1021 1031 1033 1039 1049 1051 1061 1063 1069 1087 1091 
1093 1097 1103 1109 1117 1123 1129 1151 1153 1163 1171 1181 1187 1193 
1201 1213 1217 1223 1229 1231 1237 1249 1259 1277 1279

CREATE small_primes 207 ALLOT ALIGN \ An array for small primes

\ Stopwatch timed 207 primes at 47 seconds, Miller-Rabin at 43 seconds.
100 VALUE small_primes_lmt \ Do all, or only some (207 max)?

\ Load 207 single-byte primes into array
\ Highest bit truncated for values above 255, replaced ex-post-facto
: small.primes.load 207 0 DO small_primes 206 + I - C! LOOP ;

small.primes.load \ Do that now.

: small.primes ( c-array max min -- c-array 0|1 ) \ zero-based: max:1..207, min: 0..max-1
  ."  Small Primes Division: " \ Inform of doing small primes
  1 -ROT
  DO
    a_tmp a.erase
    DUP 0= IF LEAVE THEN      \ Composite, so done
    small_primes I + C@       \ Get a small prime's lower byte
    I 171 >
    IF $400 OR                \ Know range by order in list
    ELSE    
      I 134 >                   \ Range 1031 and up
      IF $300 OR                \ Range 768-1021
      ELSE
        I 96 >
        IF $200 OR              \ Range 512-767
        ELSE
          I 53 >
          IF $100 OR THEN       \ Range 256-511
        THEN
      THEN
    THEN
    DUP . DUP >R
    a_tmp a.lsc e!            \ Load into a_tmp
    OVER a_tmp num_? a.a.mod  \ Divide
    num_? a.0<> AND           \ Any remainder?
    R> BS.u
  LOOP 
  S"  Small Primes Division: " NIP BS \ Wipe from screen
;

\ MILLER-RABIN

a.new CONSTANT exp_of_2 \ Array for current exponent of 2
a.new CONSTANT p-1
a.new CONSTANT pri_d
a.new CONSTANT pri_a

: pri.decompose ( n -- n' exp_of_2' ) \ c-addrs
  exp_of_2 a.erase
  BEGIN
     DUP num_2 num_? a.a.mod  ( n )
     num_? a.0=               ( n ) \ Remainder of zero?
  WHILE
    DUP a.rshift_1            ( n ) \ n = n/2
    1 exp_of_2 0 a.u+         ( n ) \ increment
  REPEAT
  exp_of_2                    ( n exp ) 
;

\ PRIMALITY TESTING

3 CONSTANT MIL_RAB  \ How many rounds of Miller-Rabin

: pri.d ( c-addr -- c-addr c-addr ) \ p (odd!), d
  DUP p-1 a.copy    ( p )
  p-1 a.even!       ( p )       \ p-1 set for later
  pri_d             ( p d )
  p-1 OVER a.copy   ( p d )     \ d = p-1 (even)
  BEGIN             ( p d )
    DUP a.even?     ( p d 0|1 )
  WHILE             ( p d )
    DUP a.rshift_1  ( p d )     \ d /= 2
  REPEAT
;

\ Setup for Miller-Rabin, Part 1
\ Add to stack a random: pri_a = 2 + [2...p-2] % (p-4)
\ Ref: Pick a random number 'a' in range [2, n-2]
: mil.rab.a ( p d -- p d a ) 
  a.new 2 PICK OVER a.copy    ( p d new )      \ new = p
  DUP num_2 a.a-              ( p d p-2 )      \ new = p-2
  num_2 a.new a.random.range  ( p d p-2 num_2 rand )
  -ROT 2DUP a.a- DROP         ( p d rand p-4 )
  2DUP pri_a a.a.mod          ( p d rand p-4 ) \ pri_a = rand MOD(p-4)
  a.free DROP                 ( p d )
  2 pri_a 0 a.u+ pri_a        ( p d a )
;

\ Setup for Miller-Rabin, Part 2
\ Ref 1: Compute: x = pow(a, d) % n
\ Ref 2: If x == 1 or x == n-1, return true.
: mil.rab.b ( p d a -- p d a 0|1 )  \ x = a^d % p
  2 PICK a_mod a.copy ( p d a )     \ a_mod = p
  2DUP SWAP           ( p d a a d ) 
  a.new DUP >R a.a^   ( p d a )     \ x in R          
  R@ a.1=             ( p d a 0|1 ) 
  OVER p-1 a.= OR     ( p d a 0|1 ) \ p-1 set in pri.d
  R> a.free
;

\ Single round of Miller-Rabin
: mil.rab.c ( p d i -- p 0|1 )   \ i is arbitrary indicator: 1, 2, etc.
  ."  Miller-Rabin: "            \ Inform of doing Miller-Rabin...
  DUP >R DUP . ." of 3 "         \ ...which round: 1, 2, etc.
  1 = IF                         \ First round of mil.rab ?
    mil.rab.a     ( p d a )      \ New 'a' on 1st round
  ELSE
    pri_a DUP a.rshift_1         \ a = a/2 on 2nd and 3rd (faster)
  THEN
  mil.rab.b       ( p d a 0|1 )
  NIP NIP         ( p 0|1 )
  R> BS.u
  S"  Miller-Rabin: of 3 " NIP BS
;

\ Miller-Rabin calls for plural rounds.
: mil.rab ( p -- p 0|1 )
  1 >us             ( p )        \ Init TRUE to AND with later
   MIL_RAB 1+ 1 DO  ( p )        \ N rounds of Miller-Rabin
    I 1 = IF
      pri.d         ( p d )      \ Make new on 1st round.
    ELSE
      pri_d         ( p d )      \ Reuse on 2nd and 3rd
    THEN
    I mil.rab.c     ( p 0|1 )    \ Again if passed
    us> AND DUP >us ( p 0|1 )    \ Must pass consecutively
    0= IF                        \ Exit on first failure
      LEAVE
    THEN
  LOOP              ( p d )
  us>              ( p 0|1 )    \ Consecutive results, AND'd together
;

\ Is candidate array prime?
: prime? ( c-addr -- c-addr 0|1 )    \ candidate, bool
  small_primes_lmt 
  0 small.primes      ( p 0|1 )      \ Test against from 1 to 207 primes stored as single bytes
  IF                  ( p )          \ Do Miller-Rabin test
    mil.rab           ( p 0|1 )
  ELSE
    FALSE             ( p 0 )        \ Not a prime
  THEN   
;

\ Find one prime
: find.prime ( u -- )
  CR ." Prime Candidate #"
  0
  pri_0 pri.candidate \ Load an odd random
  BEGIN
    SWAP 1+ DUP . ( show count ) SWAP
    2 OVER 0 a.u+ \ Increment by 2
    prime?
    2 PICK BS.u   \ Wipe count from screen
  UNTIL
  S" Prime Candidate #" NIP BS \ Wipe from screen
  a.show
  ." Seed: " rand_seed @ HEX U. DECIMAL
  ." Candidates: " .
;


\ Display as pictured string: 98.76%
: ##.##% ( u -- )
  0 \ Single now double
  <# [CHAR] % HOLD # # [CHAR] . HOLD #S #>
;

\ Show calculated accuracy of prime numbers
: mil.rab.pct ( - str )
  10000 10000
  MIL_RAB 0 DO 4 / LOOP
  - ##.##% TYPE
;

\ Init Blum Blum Shug seed
: bbs.seed
  pri_cells @ 2 *    \ Double the bits of primes target
  DUP >R
  0 DO
    random           \ Stack up N cells to load
  LOOP
  R> bbs_s a.u!      \ Store into BBS seed array
;
 
\ Find N primes
: find.primes ( u -- )
  rand.seed.churn \ Offset PRNG 
  ask.prime.bits  \ How big the primes 
  bbs.seed        \ Init seed for Blum Blum Shug
  pri_bits @ 1024 * 710 ( avg dist between 1024-bit primes ) / 4 / \ Avg candidates between primes
  CR
  CR ." Each odd-numbered candidate is initially vetted against the first " small_primes_lmt . ." primes. "
  CR MIL_RAB . ." rounds of Miller-Rabin follow. Accuracy thereby calculates to " mil.rab.pct ." . "
  CR ." An average spacing of " . ." odd-number candidates exists between " pri_bits @ . 8 EMIT ." -bit primes. "
  CR ." A very rare few are adjacent odd numbers. Extreme gaps separate certain others. "
  CR ." After finding a prime, the next hunt's jump-in point will be non-sequential. "
  CR
  pri_0.anew? 



  0 DO
    I find.prime
  LOOP
  pri_0.free?
;


  100 find.primes


0 [IF] ( TESTING )

: append-to-file ( addr u -- )
  file-name r/w open-file throw >r ( fileid )
  dup file-size drop rot reposition-file throw ( reposition to end )
  r@ write-file throw
  r> close-file throw ;

  
  : pri.setup
    S" 30000000091" rsa_p a.load \ P
    S" 40000000003" rsa_q a.load \ Q
    S" 6212296298"  bbs_s a.load \ BBS rand seed
    bbs.calc.n                   \ Loads rsa_n
  ;
  
  : test.pri.candidate
    CR CR ." Testing pri.candidate from file prime.fs " CR
    CR ." ENTER how many candidate primes to display: "
    PAD 3 ACCEPT PAD SWAP string.to.integer CR
    0 DO
      pri_0 pri.candidate
      CR pri_0 a.show
    LOOP
  ;
  
  : test.small.primes
    CR CR ." Testing small.primes from file prime.fs "
    CR CR ." Firstly, a definite prime: " 
    CR rsa_p a.show.bcd
    rsa_p small_primes_lmt 0 small.primes
    IF ."  has no prime factor below 1280! " ELSE ."  Composite prime found! " THEN
    DROP
    CR CR ." Of 50 random values, the following no factor below 1280:"
    50 0
    DO
      pri_0 pri.candidate
      pri_0 small_primes_lmt 0 small.primes
      IF 
        CR pri_0 a.show.bcd
      THEN
      DROP
    LOOP
    CR ." Done. " CR
  ;
  
  : test.pri.d
    CR CR ." Testing pri.d from file prime.fs " CR
    pri_0 pri.candidate
    CR ." P = " pri_0 a.show
    pri_0 pri.d
    CR ." D = " DUP a.show
    a.free 
    DROP
  ;
  
  : test.mil.rab.a
    CR CR ." Testing mil.rab.a from file prime.fs " CR
    pri_0 pri.candidate
    CR ." P = " pri_0 a.show
    pri_0 pri.d
    CR ." D = " DUP a.show
    BEGIN 
      4 0 DO
        mil.rab.a
        CR ." A = " a.show
      LOOP
      CR ." Press SPACE to continue, ENTER to quit. " CR
      KEY 13 =
    UNTIL
    2DROP  
  ;
  
  : test.mil.rab.a.x
    CR CR ." Testing mil.rab.x from file prime.fs " CR
    pri_0 pri.candidate
    CR ." P = " pri_0 a.show
    pri_0 pri.d
    CR ." D = " DUP a.show
    mil.rab.a
    CR ." A = " DUP a.show
    mil.rab.x
    CR ." PRIME? "
    IF ." YES!" ELSE ." NO!" THEN
    a.free a.free
    DROP  
  ;
  
  : test.prime?
    CR CR ." Testing prime? from file prime.fs " CR
    pri_0 pri.candidate
    CR ." Candidate = " pri_0 a.show
    pri_0 prime?
    CR ." PRIME? "
    IF ." YES!" ELSE ." NO!" THEN
  ;

\ ADAPT AND MOVE BELOW TO ABOVE TEST
  
  \ pri.setup
  \ test.small.primes
  \ test.pri.candidate
  \ test.pri.d
  \ test.mil.rab.ab
  \ test.mil.rab.a.x []
  \ test.prime?
  \ find.prime

\ https://www.mobilefish.com/services/big_number/big_number.php
\ https://www.numberempire.com/primenumbers.php


\ Verified Primes
\ CF1B3429809A08867
\ 33606C77923AF4DE57B
\ 56207F488B7823E2F6C7
\ ADF82C30F636409F34A0B
\ 1280677A5D089561175377
\ 24E2BCE2C6CA13AF228843
\ 24E2BCE2C6CA13AF2288B3
\ A7AFE553C68C21C94750F12C4412167
\ 95E784D74C7F43F0D2C7D2FAF490202AB
\ 106FDDD4F18835695B0043BADAA4092A4F
\ 118EF59E7C11B1C2B0BE96DBFF66147236FB132F
\ 13CB0355C31B5F89D1D250829EDF4DC8C7B7DE21D
\ 1607B2715FAC35BE4BAADF3966DA849521890DEFA323
\ 4DED71E114B26A648CD383BC34ACFF46C816F8C4994D81B
\ 7B5F938DA73671ED1A615911DA984E5E950A7A5B8E63349CE7F74125CA395A85431
\ 358316F279CCA7E4C7047171A9C3D23F16FC9D5222483C0DC07F14B6C8D94EE3B310364E6D7

  \ pri_0 a.free
  
[THEN]