This should get more interesting than the previous two parts. We want a usr function that tests if a number between 0 and 65535, inclusive, is a prime number. We expect an a result of this usr function either a -1 (true) or a 0 (false). I really should split this whole thing into four parts to keep it somewhat manageable (for those who follow along):

  1. thinking about the problem
  2. formulating a rough algorithm in words
  3. formulating an assembly language routine
  4. testing an confirming if the code is correct

[ part 1 | part 2 | part 4 | part 5 | bonus ]

As for the assembly language routine, I will limit myself in this part to a generator of possible prime numbers. I can use this generator to test primality of a number that was input through the usr function in Basic. But since I wanted to finish this article in a day, the completed function will have to wait until part 4.

Coding correctly takes time, after all.

Some logical reasoning about a prime number sieve algorithm

Using the article on Geeks for Geeks on prime numbers, I did some simple logical reasoning to determine an algorithm to find out whether or not a number is prime. We can easily generate prime numbers between 2 and 43, inclusive. This means we can quickly test a number’s primality up to the value 1850. After that, we need to use the formula 6 * n ± 1 for values of n larger than 7 thereafter, though not all of those are primes.

  • a prime number is a whole number larger than one that is only divisible by one and itself.
  • the number 2 is the only even prime number
  • the number 3 is the lowest odd prime number
  • the number 5 is the one but lowest odd prime number
  • the number 7 is the two but lowest odd prime number, and also the lowest consecutive odd number larger than 1 that is prime
  • any odd number not divisible by either 2, nor 3 is a candidate prime number, meaning any candidate prime number larger than 3 is of the formula:
    … 6 * n - 1, or
    … 6 * n + 1,
    … with n equal to a whole positive number
    … after having tested for the values 2, 3, 5, and 7, n must be a whole number larger than one
  • a non-prime is the product of two primes, and the largest possible product of two primes that is equal to a non-prime is the square of a prime, in other words, if a number isn’t divisible by prime numbers that are less than or equal to the square root of that number, that number is prime, in yet other words:
    … when testing a number for divisibility by primes, then highest possible prime number to test for is:
    … p <= sqrt(n), where
    … p is a prime number and n is the number to be tested for primality
  • since 7 is the lowest consecutive odd prime, we can determine the highest value of n that gives a number that is always prime, never a non-prime
    … (6 * n - 1 < 7 * 7) or (6 * n + 1 < 7 * 7),
    … (n < (7 * 7 + 1) / 6) or (n < (7 * 7 - 1) / 6)
    … n < 8
  • summarizing (some of) the above, primes between 2 and 43 can be determined by being equal to:
    … 2, 3, 5, 7
    … or calculated as:
    … 6 * n - 1, 6 * n + 1 for 1 < n < 8
    … the lowest non-prime that can’t be easily determined by using these primes is an odd number larger than the square of 43, which is 1851.

An algorithm for testing primality

So, up to a value of 1850 we can generate primes to test if a number is prime. For higher numbers we can only generate possible primes, though that isn’t a problem, since if a number is divisible by a product of primes, it isn’t prime either. It would be nice, though, to have a method to eliminate some of the non-primes to test with. The article on Geeks for Geeks mentions some additional methods to do just that. However, let us first use what is somewhat easy to understand, without a degree in mathematics.

As a first approach, using the 6 * n ± 1 formula to generate primes to test with is a solid method. We need a “next prime” function, where the input value k determines the k-th prime, with k >= 0. For primes up to 43 we can simply test if a number is equal to the output of this function, since every output number is a prime. However, we only have to test for primes which square values are smaller than the number that is being tested.

  • clear primality flag (assume it’s not a prime)
  • endless loop
    • fetch next prime
    • if number < prime * prime then exit loop
    • if number == prime then set primality flag, exit loop
  • return primality flag

If a prime number from the prime number generator function becomes larger than 43, we must stop testing for equality (since the value that is output isn’t necessarily prime), and start testing for divisibility instead. We want it to be tested as follows:

  • set primality flag (assume it’s a prime)
  • endless loop
    • fetch next prime
    • if number < prime * prime then exit loop
    • if number mod prime == 0 then it’s no prime, exit loop
  • return primality flag

To distinguish between the two:

  • if number <= 43 then compare with primes, else test with “primes”

Putting it all together

Time to fetch the code editor and start mashing keys.

Let’s first test our prime number generator. It has to output a whole positive number based on a whole non-negative number input. However, remember that we only have a function in ROM to output a signed integer, so we have to convert our unsigned integer to a signed one, and convert it back to an unsigned integer in Basic, by subtracting and adding 32768, respectively.

*= $c000

; usrfunction 0.3
; generate a candidate prime number, 
;     based on the 16-bit unsigned integer input, 
;     representing the order of the prime
; output the prime number as a signed 16-bit integer
;     to work around the limitation of the Basic ROM
;     (there's no function to convert a unsigned integer 
;      into a real number)
;     in basic, add 32768 to the output of the usr function
; for example (ignoring that the output is signed)
;     prime(0) = 2
;     prime(1) = 3
;     prime(2) = 5
; etc.

getadr = $b7f7
givayf = $b391
fac1int = $64
fac2int = $66
        jsr getadr      ;convert real into unsigned int in .A and .Y (hi/lo)
        jsr prime       ;generate the nth prime
        sec
        sbc #%10000000  ;subtract $8000 (32678) to make it a signed integer
        jmp givayf      ;make it real and return to Basic

prime:
        ;generate the nth prime number
        ;input  .A and .Y as hi/lo byte of 16-bit unsigned integer
        ;       indicating the order of the prime
        ;output .A and .Y as hi/lo byte of 16-bit unsigned integer
        ;       containing the nth prime number on the number line
        ;e.g. prime(0) -> 2, prime(1) -> 3, prime(2) -> 5, etc.
        ;only numbers generated up to 43 are guaranteed prime,
        ;    larger numbers may be prime, 
        ;    or co-prime (product of 2 or more primes)
        cmp #0
        bne prime2      ;order is too big, skip 0 and 1
        cpy #0          ;order == 0?
        bne prime1      ;no, check for 1
        ldy #2          ;prime(0) = 2
        rts
prime1:
        cpy #1          ;order == 1?
        bne prime2      ;no, check for 2 and up
        ldy #3          ;prime(1) = 3
        rts
prime2:
        ;calculate primes for (order > 1)
        ;in the formula (6 * n ± 1), where n is equal to (order / 2),
        lsr             ;n = (order / 2)
        sta fac1int
        tya
        ror             
        sta fac1int+1
        asl fac1int+1   ;k = n * 2 
        rol fac1int
        lda fac1int     ;m = n * 2
        ldx fac1int+1
        sta fac2int
        stx fac2int+1
        asl fac2int+1   ;m = n * 4
        rol fac2int
        clc             ;add m to k, to get (6 * n)
        lda fac1int+1
        adc fac2int+1
        sta fac1int+1
        lda fac1int
        adc fac2int
        sta fac1int
        tya             ;retrieve original lo byte of order
        and #1          ;is it odd?
        cmp #1
        beq prime3      ;then add 1 to (6 * n)
        sec             ;else subtract 1 from (6 * n)
        lda fac1int+1
        sbc #1
        tay             ;lo byte of ((6 * n) - 1) in .Y
        lda fac1int
        sbc #0          ;hi byte of ((6 * n) - 1) in .A
        rts
prime3:
        clc             ;add 1 to (6 * n)
        lda fac1int+1
        adc #1
        tay             ;lo byte of ((6 * n) + 1) in .Y
        lda fac1int
        adc #0          ;hi byte of ((6 * n) + 1) in .A
        rts

Give our code a spin, testing its validity

I pasted the source code into Virtual 6502 Assembler, copied the hexdump, and converted with a calculator app into decimal values of data statements in the Basic program “testusrf 0.3” below.

0 rem testusrf 0.3
10 poke 785,0:poke 786,192:rem set usr address to $c000
20 forn=49152to49246:readb:poken,b:next:rem read in machine code
30 data 032,247,183,032,012,192,056,233
31 data 128,076,145,179,201,000,208,014
32 data 192,000,208,003,160,002,096,192
33 data 001,208,003,160,003,096,074,133
34 data 100,152,106,133,101,006,101,038
35 data 100,165,100,166,101,133,102,134
36 data 103,006,103,038,102,024,165,101
37 data 101,103,133,101,165,100,101,102
38 data 133,100,152,041,001,201,001,240
39 data 011,056,165,101,233,001,168,165
40 data 100,233,000,096,024,165,101,105
41 data 001,168,165,100,105,000,096
50 for x = 0 to 9:gosub 100:next
99 end
100 print "x =";x;", usr(x) =";usr(x)+32768:return

Running this program for the first 100 “primes” gave this output:

list 50

50 for x=00000 to 00099:gosub 100:next

ready.
run
 2  3  5  7  11  13  17  19  23  25  29
 31  35  37  41  43  47  49  53  55  59
 61  65  67  71  73  77  79  83  85  89
 91  95  97  101  103  107  109  113  11
5  119  121  125  127  131  133  137  13
9  143  145  149  151  155  157  161  16
3  167  169  173  175  179  181  185  18
7  191  193  197  199  203  205  209  21
1  215  217  221  223  227  229  233  23
5  239  241  245  247  251  253  257  25
9  263  265  269  271  275  277  281  28
3  287  289  293  295
ready.

I noticed this:

print usr(21845)+32768,usr(21846)+32768
 65533     1

ready.

This is because the prime number generator doesn’t check for any value overflow. The value 1 is actually the value 65537, but the 17th bit of the latter number doesn’t exist, so the value wraps around to 1.

print (65533 - 1)/6
 10922

ready.

So the prime function gives valid answers up to the 10922th prime, which is more than enough to test if the largest 16-bit unsigned value (65535) is prime. For that we need at most a prime of the square root of that value, which is slightly less than 256. These “primes” are already visible in the first 100 “primes”.

Did I do it right?

To check, you can compare the output of the prime number generator with the List of Prime Numbers on Wikipedia. We are only concerned with primes less than 256, since it’s the square root of the value 65536. 65536 is already higher than the maximal value of a 16-bit unsigned integer number (65535). To speed things up, we could’ve just as well hard-coded the list of 54 primes that are less than 256. The routine would need more bytes, though.


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

The prime numbers from the Wikipedia page are all included in the output of “testusrf 0.3” in the listing above. The prime generator does output some extra (non-prime) numbers, that don’t do any harm when testing for primality.