Legendre prime counting function: Difference between revisions

Content added Content deleted
(Added C)
(→‎{{header|Vlang}}: correction to increase precision to extend range...)
Line 3,099: Line 3,099:


[inline]
[inline]
fn half(n int) int { return (n - 1) >> 1 } // convenience function!
fn half(n u64) i64 { return i64((n - 1) >>> 1) } // convenience function!
[inline] // floating point divide faster than integer divide!
[inline] // floating point divide faster than integer divide!
fn divide(nm u64, d u64) int { return int(f64(nm) / f64(d)) }
fn divide(nm u64, d u64) u64 { return u64(f64(nm) / f64(d)) }


[direct_array_access]
[direct_array_access]
fn count_primes_to(n u64) i64 {
fn count_primes_to(n u64) i64 {
if n < u64(3) { return if n < i64(2) { i64(0) } else { i64(1) } }
if n < u64(3) { return if n < i64(2) { i64(0) } else { i64(1) } }
rtlmt := int(sqrt(f64(n)))
rtlmt := u64(sqrt(f64(n)))
mxndx := (rtlmt - 1) / 2
mxndx := i64((rtlmt - 1) / 2)
arrlen := mxndx + 1
arrlen := int(mxndx + 1)
mut smalls := []u32{ len: arrlen, cap: arrlen, init: u32(it) }
mut smalls := []u32{ len: arrlen, cap: arrlen, init: u32(it) }
mut roughs := []u32{ len: arrlen, cap: arrlen, init: u32(it + it + 1) }
mut roughs := []u32{ len: arrlen, cap: arrlen, init: u32(it + it + 1) }
mut larges := []i64{ len: arrlen, cap: arrlen,
mut larges := []i64{ len: arrlen, cap: arrlen,
init: i64(((n / u64(it + it + 1) - 1) / 2)) }
init: i64(((n / u64(it + it + 1) - 1) / 2)) }
cullbuflen := (mxndx + 8) / 8
cullbuflen := int((mxndx + 8) / 8)
mut cullbuf := []u8{len: cullbuflen, cap: cullbuflen, init: u8(0)}
mut cullbuf := []u8{len: cullbuflen, cap: cullbuflen, init: u8(0)}


mut nbps := 0 // do partial sieving for each odd base prime to quad root...
mut nbps := i64(0) // do partial sieving for each odd base prime to quad root...
mut rilmt := arrlen // number of roughs start with full size!
mut rilmt := arrlen // number of roughs start with full size!
for i := 1; ; i++ { // start with 3; breaks when partial sieving done...
for i := i64(1); ; i++ { // start with 3; breaks when partial sieving done...
sqri := (i + i) * (i + 1)
sqri := (i + i) * (i + 1)
if sqri > mxndx { break } // breaks here!
if sqri > mxndx { break } // breaks here!
if (cullbuf[i >> 3] & masks[i & 7]) != u8(0) { continue } // not prime!
if (cullbuf[i >>> 3] & masks[i & 7]) != u8(0) { continue } // not prime!
cullbuf[i >> 3] |= masks[i & 7] // cull bp rep!
cullbuf[i >>> 3] |= masks[i & 7] // cull bp rep!
bp := i + i + 1 // partial sieve by bp...
bp := u64(i + i + 1) // partial sieve by bp...
for c := sqri; c < arrlen; c += bp { cullbuf[c >> 3] |= masks[c & 7] }
for c := sqri; c < arrlen; c += i64(bp) { cullbuf[c >>> 3] |= masks[c & 7] }


mut nri := 0 // transfer from `ori` to `nri` indexes:
mut nri := 0 // transfer from `ori` to `nri` indexes:
for ori in 0 .. rilmt { // update `roughs` and `larges` for partial sieve...
for ori in 0 .. rilmt { // update `roughs` and `larges` for partial sieve...
q := int(roughs[ori])
r := roughs[ori]
pi := q >> 1 // skip recently culled in last partial sieve...
rci := i64(r >>> 1) // skip recently culled in last partial sieve...
if (cullbuf[pi >> 3] & masks[pi & 7]) != u8(0) { continue }
if (cullbuf[rci >>> 3] & masks[rci & 7]) != u8(0) { continue }
d := q * bp
d := u64(r) * u64(bp)
larges[nri] = larges[ori] -
larges[nri] = larges[ori] -
(if d <= rtlmt { larges[int(smalls[d >> 1]) - nbps] }
(if d <= rtlmt { larges[i64(smalls[d >>> 1]) - nbps] }
else { i64(smalls[half(divide(n, u64(d)))]) })
else { i64(smalls[half(divide(n, d))]) })
+ i64(nbps) // adjust for over subtraction of base primes!
+ nbps // adjust for over subtraction of base primes!
roughs[nri] = u32(q) // compress for culled `larges` and `roughs`!
roughs[nri] = r // compress for culled `larges` and `roughs`!
nri++ // advance output `roughs` index!
nri++ // advance output `roughs` index!
}
}
Line 3,142: Line 3,142:
mut si := mxndx // update `smalls` for partial sieve...
mut si := mxndx // update `smalls` for partial sieve...
for pm := ((rtlmt / bp) - 1) | 1; pm >= bp; pm -= 2 {
for pm := ((rtlmt / bp) - 1) | 1; pm >= bp; pm -= 2 {
c := smalls[pm >> 1]
c := smalls[pm >>> 1]
e := (pm * bp) >> 1
e := (pm * bp) >>> 1
for ; si >= e; si-- { smalls[si] -= (c - u32(nbps)) }
for ; si >= e; si-- { smalls[si] -= (c - u32(nbps)) }
}
}
Line 3,159: Line 3,159:
p := u64(roughs[ri])
p := u64(roughs[ri])
m := n / p
m := n / p
ei := int(smalls[half(int(m / p))]) - nbps
ei := int(smalls[half(u64(m / p))]) - nbps
if ei <= ri { break } // break out when only multiple equal to `p`
if ei <= ri { break } // break out when only multiple equal to `p`
ans -= i64((ei - ri) * (nbps + ri - 1)) // adjust for over addition below!
ans -= i64((ei - ri) * (nbps + ri - 1)) // adjust for over addition below!
Line 3,189: Line 3,189:
This took 2 milliseconds.</pre>
This took 2 milliseconds.</pre>


The above code is about as fast as the same algorithm in the fastest of languages such as C/C++/Nim due to the greatly reduced computational complexity of O(n**(3/4)/((log n)**2)); as shown it is run on an Intel i5-6500 (3.6 GHz single-thread boosted). The above code can find the number of primes to a trillion (1e12) in about a hundred milliseconds, but has a segmentation fault for values much larger than that for unknown reasons, as the code has been translated reasonably faithfully from Nim and Crystal, perhaps due to limitations in the size of objects that the default garbage collection can handle.
The above code is about as fast as the same algorithm in the fastest of languages such as C/C++/Nim due to the greatly reduced computational complexity of O(n**(3/4)/((log n)**2)); it can compute the number of primes to 1e11 in about 22 milliseconds on the same machine as shown where it is run on an Intel i5-6500 (3.6 GHz single-thread boosted). The above code can compute the number of primes to 1e16 in about a minute, although there is a precision/accuracy problem for counting ranges above about 9e15 due to the speed-up optimization of using floating point division operations (precise only to 53-bits, which is about this limit).


=={{header|Wren}}==
=={{header|Wren}}==
Retrieved from "https://rosettacode.org/wiki/Legendre_prime_counting_function"