Legendre prime counting function: Difference between revisions
→{{header|Wren}}: Add non-memoizing partial-sieving V language version above...
GordonBGood (talk | contribs) (→{{header|Haskell}}: make memoized code compile; add non-memoized versions...) |
GordonBGood (talk | contribs) (→{{header|Wren}}: Add non-memoizing partial-sieving V language version above...) |
||
Line 2,022:
π(10^8) = 5761455
π(10^9) = 50847534</pre>
=={{header|Vlang}}==
There doesn't seem to be much point to contributing memoized (Hash table) versions as they seem to be the result of the task author not realizing that there is a simple way of keeping the computation complexity from having exponential growth with counting range by limiting the "splitting" of the "phi" calculation when the result must be just one. Currently, anonymous functions and closure recursion doesn't work very well (which is why the V language is still in beta) so those slower algorithms don't work properly. However, the "partial sieving" Nim version has been [translated from the non-memoizing Nim version](https://rosettacode.org/wiki/Legendre_prime_counting_function#Non-Memoized_Versions), which explanations can be referenced for an understanding of how it works:
'''The Non-Memoizing Non-Recursive "Partial-Sieving" Algorithm'''
<syntaxhighlight lang="vlang">// compile with: v -cflags -march=native -cflags -O3 -prod MagicCount.v
import time
import math { sqrt }
const masks = [ u8(1), 2, 4, 8, 16, 32, 64, 128 ] // faster than bit twiddling!
[inline]
fn half(n int) int { return (n - 1) >> 1 } // convenience function!
[inline] // floating point divide faster than integer divide!
fn divide(nm u64, d u64) int { return int(f64(nm) / f64(d)) }
[direct_array_access]
fn count_primes_to(n u64) i64 {
if n < u64(3) { return if n < i64(2) { i64(0) } else { i64(1) } }
rtlmt := int(sqrt(f64(n)))
mxndx := (rtlmt - 1) / 2
arrlen := mxndx + 1
mut smalls := []u32{ len: arrlen, cap: arrlen, init: u32(it) }
mut roughs := []u32{ len: arrlen, cap: arrlen, init: u32(it + it + 1) }
mut larges := []i64{ len: arrlen, cap: arrlen,
init: i64(((n / u64(it + it + 1) - 1) / 2)) }
cullbuflen := (mxndx + 8) / 8
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 rilmt := arrlen // number of roughs start with full size!
for i := 1; ; i++ { // start with 3; breaks when partial sieving done...
sqri := (i + i) * (i + 1)
if sqri > mxndx { break } // breaks here!
if (cullbuf[i >> 3] & masks[i & 7]) != u8(0) { continue } // not prime!
cullbuf[i >> 3] |= masks[i & 7] // cull bp rep!
bp := i + i + 1 // partial sieve by bp...
for c := sqri; c < arrlen; c += bp { cullbuf[c >> 3] |= masks[c & 7] }
mut nri := 0 // transfer from `ori` to `nri` indexes:
for ori in 0 .. rilmt { // update `roughs` and `larges` for partial sieve...
q := int(roughs[ori])
pi := q >> 1 // skip recently culled in last partial sieve...
if (cullbuf[pi >> 3] & masks[pi & 7]) != u8(0) { continue }
d := q * bp
larges[nri] = larges[ori] -
(if d <= rtlmt { larges[int(smalls[d >> 1]) - nbps] }
else { i64(smalls[half(divide(n, u64(d)))]) })
+ i64(nbps) // adjust for over subtraction of base primes!
roughs[nri] = u32(q) // compress for culled `larges` and `roughs`!
nri++ // advance output `roughs` index!
}
mut si := mxndx // update `smalls` for partial sieve...
for pm := ((rtlmt / bp) - 1) | 1; pm >= bp; pm -= 2 {
c := smalls[pm >> 1]
e := (pm * bp) >> 1
for ; si >= e; si-- { smalls[si] -= (c - u32(nbps)) }
}
rilmt = nri // new rough size limit!
nbps++ // for one partial sieve base prime pass!
}
// combine results so far; adjust for over subtraction of base prime count...
mut ans := larges[0] + i64(((rilmt + 2 * (nbps - 1)) * (rilmt - 1) / 2))
for ri in 1 .. rilmt { ans -= larges[ri] } // combine results so far!
// add quotient for product of pairs of primes, quad root to cube root...
for ri := 1; ; ri++ { // break controlled below
p := u64(roughs[ri])
m := n / p
ei := int(smalls[half(int(m / p))]) - nbps
if ei <= ri { break } // break out when only multiple equal to `p`
ans -= i64((ei - ri) * (nbps + ri - 1)) // adjust for over addition below!
for sri in ri + 1 .. ei + 1 { // add for products of 1st and 2nd primes...
ans += i64(smalls[half(divide(m, u64(roughs[sri])))]) }
}
return ans + 1 // add one for only even prime of two!
}
start_time := time.now()
mut pow := u64(1)
for i in 0 .. 10 {
println("π(10**$i) = ${count_primes_to(pow)}")
pow *= 10 }
duration := (time.now() - start_time).milliseconds()
println("This took $duration milliseconds.")</syntaxhighlight>
{{out}}
<pre>π(10**0) = 0
π(10**1) = 4
π(10**2) = 25
π(10**3) = 168
π(10**4) = 1229
π(10**5) = 9592
π(10**6) = 78498
π(10**7) = 664579
π(10**8) = 5761455
π(10**9) = 50847534
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.
=={{header|Wren}}==
| |||