60 lines
1.8 KiB
Chapel
60 lines
1.8 KiB
Chapel
// compile with --fast for maximum speed...
|
|
|
|
use Time;
|
|
|
|
proc countPrimes(lmt: uint(64)): int(64) {
|
|
if lmt < 9 { // when there are no odd primes less than square root...
|
|
if lmt < 3 { if lmt < 2 { return 0; } else { return 1; } }
|
|
return (lmt - (lmt >> 1)): int(64);
|
|
}
|
|
|
|
// Chapel doesn't have closures, so emulate them with a class...
|
|
class LegendrePi {
|
|
var n: uint(64);
|
|
var dom: domain(1);
|
|
var oprms: [dom] uint(32);
|
|
proc init(n: uint(64)) {
|
|
// first, an array of odd primes to the square root of n is generated...
|
|
this.n = n;
|
|
const sqrtn = sqrt(n: real(64)): int(64);
|
|
const rtlmt = (sqrtn - 3) / 2; this.dom = {0 .. rtlmt};
|
|
this.oprms = 0;
|
|
for i in 0 .. rtlmt do this.oprms[i] = (i + i + 3): uint(32);
|
|
var i = 0;
|
|
for i in (0 ..) { // cull the array
|
|
var ci = (i + i) * (i + 3) + 3; if ci > rtlmt { break; }
|
|
const bp = i + i + 3;
|
|
while (ci <= rtlmt) { this.oprms[ci] = 0; ci += bp; }
|
|
}
|
|
var psz = 0;
|
|
for ti in 0 .. rtlmt { // compress the odd primes array...
|
|
const tv = this.oprms[ti];
|
|
if tv != 0 { this.oprms[psz] = tv; psz += 1; }
|
|
}
|
|
this.dom = { 0 ..< psz };
|
|
}
|
|
proc phi(x: uint(64), a: int): int(64) {
|
|
if a <= 0 { return (x - (x >> 1)): int(64); } // take care of prime of 2
|
|
const na = a - 1; const p = this.oprms[na]: uint(64);
|
|
if x <= p { return 1: int(64); }
|
|
return phi(x, na) - phi(x / p, na);
|
|
}
|
|
proc this(): int(64) {
|
|
return phi(n, this.oprms.size) + this.oprms.size: int(64);
|
|
}
|
|
}
|
|
return (new LegendrePi(lmt))();
|
|
}
|
|
|
|
proc main() {
|
|
var timer: Timer;
|
|
timer.start();
|
|
|
|
for i in 0 .. 9 {
|
|
writeln("π(10**", i, ") = ", countPrimesx(10: uint(64) ** i));
|
|
}
|
|
|
|
timer.stop();
|
|
|
|
writeln("This took ", timer.elapsed(TimeUnits.milliseconds), " milliseconds.");
|
|
}
|