Rebol3 Code Examplex
Legendre prime counting function
Count primes up to a limit using Legendre’s function.
Rebol [
title: "Rosetta code: Legendre prime counting function"
file: %Legendre_prime_counting_function.r3
url: https://rosettacode.org/wiki/Legendre_prime_counting_function
]
legendre-pi-count: function/with [
"Count primes up to n using Legendre's phi formula"
n [number!]
][
if n == 0 [return 0]
if n == 2 [return 1]
primes: primes-up-to square-root n
prime-count: length? primes
result: (phi n prime-count) + prime-count - 1
][
primes: sieve: _ ;; shared state across calls
cache: make map! 100000 ;; phi memoization table
prime-sieve: function [
;; Returns a bitset where bit i is true if i is prime
sieve-size [integer!]
][
sieve: make bitset! sieve-size
repeat i sieve-size - 1 [sieve/:i: true]
sieve/1: false ;; 1 is not prime
repeat s-pos sieve-size [
if sieve/:s-pos [
p: s-pos * s-pos ;; start crossing off at p^2
while [p <= sieve-size] [
sieve/:p: false
p: p + s-pos
]
]
]
sieve
]
sieve: prime-sieve to integer! sqrt 1e9 ;; one-time sieve up to sqrt(10^9)
primes-up-to: function [
;; Extract list of primes from sieve up to n
n [number!]
][
out: clear []
repeat i n [ if sieve/:i [append out i] ]
out
]
phi: function [
;; Legendre's phi: count integers <= x not divisible by first a primes
x a
][
key: as-pair x a ;; pack (x, a) into a pair! for cache lookup
unless result: cache/:key [
sum: 0
while [a > 1] [
either x <= pa: primes/:a [
done?: true ;; x has no more factors to sieve
break
][
-- a
sum: sum + phi x // pa a
]
]
result: either done? [1] [x - sum - (x // 2)]
cache/:key: result ;; memoize before returning
]
result
]
]
start: stats/timer
for exp 0 9 1 [
print ajoin ["π(10**" exp ") = " legendre-pi-count 10 ** exp]
]
print ["This took" stats/timer - start]