_top      = 100000000   // 100m is largest magnitude this code can handle
_max      = 11000       // Number of primes in list
_prime    = 0
_notPrime = 1
uint16    primes(_max)
uint8     bits(_top / 7)
ptr       gBase : gBase = @bits( 0 )


void local fn init
print @"\n   A brilliant number is the product of two primes of the same magnitude."
print @"\n                         First 100 brilliant numbers:"

ptr     p = @primes(0)
uint64 n1 = 2, n2 = 2

while n1 < _max               // Create list of primes via sieve of Eratosthenes
if primes(n1) == _prime
% p, n1                    // add n1 to list of primes
p += 2
n2 = n1 << 1
while n2 < _max
primes(n2) = _notPrime   // mark all multiples of n1 as not prime
n2 += n1
wend
end if
n1 ++
wend

end fn


void local fn brilliants
ptr pn1 = @primes(0), pn2 = pn1
uint64 mag = 10, n, limit = _top / 3

do

while {pn2} < mag
n = {pn1} * {pn2}  // mark product of 2 primes as brilliant
if n >= limit then exit fn
bits( n >> 3 ) =  bits( n >> 3 ) or bit( n and 7 )
pn2 += 2
wend

pn1 += 2 : pn2 = pn1
if {pn1} > mag then mag *= 10 : if mag == _top then mag = limit

until 0
stop
end fn


void local fn show
uint64 count = 0, mag = 1, b, n, v
ptr    p     = gBase

while count < 100
v = peek int( p )
b = 0

while v
if v and (1 << b)
v -= (1 << b)
n = ((p - gBase) << 3) + b
count++
print fn StringWithFormat(@"%7d", n);
if count mod 10 == 0 then print

if n >= mag
mda_add 2 = count : mda_add = n
mag *= 10
end if
end if
b++
wend

p += 4
wend

//============================

while mag < _top           //search for next magnitude

while p < @bits(mag >> 3) // count up to byte before magnitude
v = [p]

while v
count++
v = v and (v-1)       // clear lowest 1 bit in v
wend

p += 8                     // move to next 64 bits
wend
// least-of-magnitude will be next number found

while peek int( p ) == 0 : p += 4 : wend

//v = peek int( p )
b = 0

while peek int( p )
if peek int( p )  and (1 << b)
count++
poke int p, peek int( p ) - (1 << b) // Bit has been counted so remove it

// get value of brilliant number
mda_add   = @(((p - gBase) << 3) + b)  //add it to array
mda_add 2 = count
mag      *= 10
end if
b++
wend

p += 4
wend

//============================

print @"\n                First brilliant number of each magnitude:"
count = 0

while mda(count)
printf @"%26d is brilliant number #%d", mda_integer(count), mda_integer 2(count)
count++
wend

end fn


window 1, @"Brilliant Numbers"
CFTimeInterval t : t = fn CACurrentMediaTime
fn init
fn brilliants
fn show
printf @"\n   Time = %.2f milliseconds.", 1000 * (fn CACurrentMediaTime - t)

handleevents