Advertisement
Not a member of Pastebin yet?
Sign Up,
it unlocks many cool features!
- Input Output
- 90 24
- 3000 430
- 9000 1117
- 4000000 283146 <--- input = 4*10^6
- 800000000 41146179 <--- input = 9*10^8
- 1100000000 55662470 <--- input = 1.1*10^9
- Input Output
- 1907000000 93875448 <--- input = 1.907*10^9
- 1337000000 66990613 <--- input = 1.337*10^9
- 1240000000 62366021 <--- input = 1.24*10^9
- 660000000 34286170 <--- input = 6.6*10^8
- 99820000 5751639 <--- input = 9.982*10^7
- 40550000 2465109 <--- input = 4.055*10^7
- 24850000 1557132 <--- input = 2.485*10^7
- 41500 4339
- #!/bin/bash
- a=(1907000000 1337000000 1240000000 660000000 99820000 40550000 24850000 41500)
- echo DennisC
- exec 2>> times/dennisc.txt
- time for j in ${a[@]}; do ./dennisc $j; done >> /dev/null;
- echo DennisPy
- exec 2>> times/dennispy.txt
- time for j in ${a[@]}; do pypy dennispy.py <<< $j; done >> /dev/null;
- echo arjandelumens
- exec 2>> times/arjandelumens.txt
- time for j in ${a[@]}; do ./arjandelumens $j; done >> /dev/null;
- echo orlp
- exec 2>> times/orlp.txt
- time for j in ${a[@]}; do ./orlp $j; done >> /dev/null;
- # echo mwr247
- # time node-v4.3.1-linux-x64/bin/node mwr247.js
- # mwr247 using js seems a bit longer, so I am going to run the fastest
- # and then come back to his.
- # mwr247 provided a function, so I appended
- # console.log( F( <argument> ) )
- # to his code, for each argument.
- $ ./timepi '-march=native -O3' pi 1000
- pi(x) = 93875448 real time: 2774958 ns
- pi(x) = 66990613 real time: 2158491 ns
- pi(x) = 62366021 real time: 2023441 ns
- pi(x) = 34286170 real time: 1233158 ns
- pi(x) = 5751639 real time: 384284 ns
- pi(x) = 2465109 real time: 239783 ns
- pi(x) = 1557132 real time: 196248 ns
- pi(x) = 4339 real time: 60597 ns
- 0.00572879 s
- real 0m28.006s
- user 0m15.703s
- sys 0m14.319s
- real 0m8.934s
- user 0m8.795s
- sys 0m0.150s
- real 0m8.956s
- user 0m8.818s
- sys 0m0.150s
- real 0m8.907s
- user 0m8.775s
- sys 0m0.145s
- real 0m8.904s
- user 0m8.775s
- sys 0m0.141s
- real 0m8.902s
- user 0m8.783s
- sys 0m0.132s
- real 0m9.087s
- user 0m8.923s
- sys 0m0.176s
- real 0m8.905s
- user 0m8.778s
- sys 0m0.140s
- real 0m9.005s
- user 0m8.859s
- sys 0m0.158s
- real 0m8.911s
- user 0m8.789s
- sys 0m0.135s
- real 0m8.907s
- user 0m8.781s
- sys 0m0.138s
- $ time for i in 1.907e9 1.337e9 1.24e9 6.6e8 9.982e7 4.055e7 2.485e7 41500
- > do pypy pi.py <<< $i; done
- 93875448
- 66990613
- 62366021
- 34286170
- 5751639
- 2465109
- 1557132
- 4339
- real 0m1.696s
- user 0m1.360s
- sys 0m0.332s
- real 3m56.961s
- user 3m38.802s
- sys 0m18.512s
- 0000000: cafe babe 0000 0034 0030 0a00 0900 1709 .......4.0......
- 0000010: 0018 0019 0a00 1a00 1b0a 0008 001c 0a00 ................
- 0000020: 1d00 1e0a 0008 001f 0a00 2000 2107 0022 .......... .!.."
- 0000030: 0700 2301 0006 3c69 6e69 743e 0100 0328 ..#...<init>...(
- 0000040: 2956 0100 0443 6f64 6501 000f 4c69 6e65 )V...Code...Line
- 0000050: 4e75 6d62 6572 5461 626c 6501 0004 6d61 NumberTable...ma
- 0000060: 696e 0100 1628 5b4c 6a61 7661 2f6c 616e in...([Ljava/lan
- 0000070: 672f 5374 7269 6e67 3b29 5601 0008 6e75 g/String;)V...nu
- 0000080: 6d50 7269 6d65 0100 0428 4929 4901 000d mPrime...(I)I...
- 0000090: 5374 6163 6b4d 6170 5461 626c 6501 0007 StackMapTable...
- 00000a0: 6973 5072 696d 6501 0004 2849 295a 0100 isPrime...(I)Z..
- 00000b0: 0a53 6f75 7263 6546 696c 6501 0006 452e .SourceFile...E.
- 00000c0: 6a61 7661 0c00 0a00 0b07 0024 0c00 2500 java.......$..%.
- 00000d0: 2607 0027 0c00 2800 290c 0010 0011 0700 &..'..(.).......
- 00000e0: 2a0c 002b 002c 0c00 1300 1407 002d 0c00 *..+.,.......-..
- 00000f0: 2e00 2f01 0001 4501 0010 6a61 7661 2f6c ../...E...java/l
- 0000100: 616e 672f 4f62 6a65 6374 0100 106a 6176 ang/Object...jav
- 0000110: 612f 6c61 6e67 2f53 7973 7465 6d01 0003 a/lang/System...
- 0000120: 6f75 7401 0015 4c6a 6176 612f 696f 2f50 out...Ljava/io/P
- 0000130: 7269 6e74 5374 7265 616d 3b01 0011 6a61 rintStream;...ja
- 0000140: 7661 2f6c 616e 672f 496e 7465 6765 7201 va/lang/Integer.
- 0000150: 0008 7061 7273 6549 6e74 0100 1528 4c6a ..parseInt...(Lj
- 0000160: 6176 612f 6c61 6e67 2f53 7472 696e 673b ava/lang/String;
- 0000170: 2949 0100 136a 6176 612f 696f 2f50 7269 )I...java/io/Pri
- 0000180: 6e74 5374 7265 616d 0100 0770 7269 6e74 ntStream...print
- 0000190: 6c6e 0100 0428 4929 5601 000e 6a61 7661 ln...(I)V...java
- 00001a0: 2f6c 616e 672f 4d61 7468 0100 0473 7172 /lang/Math...sqr
- 00001b0: 7401 0004 2844 2944 0021 0008 0009 0000 t...(D)D.!......
- 00001c0: 0000 0004 0001 000a 000b 0001 000c 0000 ................
- 00001d0: 001d 0001 0001 0000 0005 2ab7 0001 b100 ..........*.....
- 00001e0: 0000 0100 0d00 0000 0600 0100 0000 0100 ................
- 00001f0: 0900 0e00 0f00 0100 0c00 0000 2c00 0300 ............,...
- 0000200: 0100 0000 10b2 0002 2a03 32b8 0003 b800 ........*.2.....
- 0000210: 04b6 0005 b100 0000 0100 0d00 0000 0a00 ................
- 0000220: 0200 0000 0300 0f00 0400 0a00 1000 1100 ................
- 0000230: 0100 0c00 0000 6600 0200 0300 0000 2003 ......f....... .
- 0000240: 3c03 3d1c 1aa2 0018 1b1c b800 0699 0007 <.=.............
- 0000250: 04a7 0004 0360 3c84 0201 a7ff e91b ac00 .....`<.........
- 0000260: 0000 0200 0d00 0000 1600 0500 0000 0600 ................
- 0000270: 0200 0700 0900 0800 1800 0700 1e00 0900 ................
- 0000280: 1200 0000 1800 04fd 0004 0101 5001 ff00 ............P...
- 0000290: 0000 0301 0101 0002 0101 fa00 0700 0a00 ................
- 00002a0: 1300 1400 0100 0c00 0000 9700 0300 0300 ................
- 00002b0: 0000 4c1a 05a2 0005 03ac 1a05 9f00 081a ..L.............
- 00002c0: 06a0 0005 04ac 1a05 7099 0009 1a06 709a ........p.....p.
- 00002d0: 0005 03ac 1a87 b800 078e 0460 3c10 063d ...........`<..=
- 00002e0: 1c1b a300 1b1a 1c04 6470 9900 0b1a 1c04 ........dp......
- 00002f0: 6070 9a00 0503 ac84 0206 a7ff e604 ac00 `p..............
- 0000300: 0000 0200 0d00 0000 2200 0800 0000 0c00 ........".......
- 0000310: 0700 0d00 1300 0e00 2100 0f00 2a00 1000 ........!...*...
- 0000320: 3200 1100 4400 1000 4a00 1200 1200 0000 2...D...J.......
- 0000330: 1100 0907 0901 0b01 fd00 0b01 0114 01fa ................
- 0000340: 0005 0001 0015 0000 0002 0016 ............
- > time java E 41500;time java E 24850000;time java E 40550000;time java E 99820000;time java E 660000000;time java E 1240000000;time java E 1337000000;time java E 1907000000
- 4339
- real 0m0.236s
- user 0m0.112s
- sys 0m0.024s
- 1557132
- real 0m8.842s
- user 0m8.784s
- sys 0m0.060s
- 2465109
- real 0m18.442s
- user 0m18.348s
- sys 0m0.116s
- 5751639
- real 1m15.642s
- user 1m8.772s
- sys 0m0.252s
- 34286170
- real 40m35.810s
- user 16m5.240s
- sys 0m5.820s
- 62366021
- real 104m12.628s
- user 39m32.348s
- sys 0m13.584s
- 66990613
- real 110m22.064s
- user 42m28.092s
- sys 0m11.320s
- 93875448
- real 171m51.650s
- user 68m39.968s
- sys 0m14.916s
- use std::env;
- const EMPTY: usize = std::usize::MAX;
- const MAX_X: usize = 800;
- fn main() {
- let args: Vec<_> = env::args().collect();
- let x: usize = args[1].trim().parse().expect("expected a number");
- let root = (x as f64).sqrt() as usize;
- let y = (x as f64).powf(0.3333333333333) as usize + 1;
- let sieve_size = x / y + 2;
- let mut sieve = vec![true; sieve_size];
- let mut primes = vec![0; sieve_size];
- sieve[0] = false;
- sieve[1] = false;
- let mut a = 0;
- let mut num_primes = 1;
- let mut num_primes_smaller_root = 0;
- // find all primes up to x/y ~ x^2/3 aka sieve_size
- for i in 2..sieve_size {
- if sieve[i] {
- if i <= root {
- if i <= y {
- a += 1;
- }
- num_primes_smaller_root += 1;
- }
- primes[num_primes] = i;
- num_primes += 1;
- let mut multiples = i;
- while multiples < sieve_size {
- sieve[multiples] = false;
- multiples += i;
- }
- }
- }
- let interesting_primes = primes[a + 1..num_primes_smaller_root + 1].iter();
- let p_2 =
- interesting_primes
- .map(|ip| primes.iter().take_while(|&&p| p <= x / ip).count())
- .enumerate()
- .map(|(i, v)| v - 1 - i - a)
- .fold(0, |acc, v| acc + v);
- let mut phi_results = vec![EMPTY; (a + 1) * MAX_X];
- println!("pi({}) = {}", x, phi(x, a, &primes, &mut phi_results) + a - 1 - p_2);
- }
- fn phi(x: usize, b: usize, primes: &[usize], phi_results: &mut [usize]) -> usize {
- if b == 0 {
- return x;
- }
- if x < MAX_X && phi_results[x + b * MAX_X] != EMPTY {
- return phi_results[x + b * MAX_X];
- }
- let value = phi(x, b - 1, primes, phi_results) - phi(x / primes[b], b - 1, primes, phi_results);
- if x < MAX_X {
- phi_results[x + b * MAX_X] = value;
- }
- value
- }
- #!/bin/bash
- a=(1907000000 1337000000 1240000000 660000000 99820000 40550000 24850000 41500)
- for i in {0..100}; do
- for j in ${a[@]}; do
- ./target/release/pi_n $j > /dev/null;
- done;
- done;
- [me@localhost pi_n]$ time ./time.sh
- real 0m37.011s
- user 0m34.752s
- sys 0m2.410s
- #include <cmath>
- #include <cstdint>
- #include <iostream>
- #include <vector>
- #include <boost/dynamic_bitset.hpp>
- uint64_t num_primes(uint64_t n) {
- // http://stackoverflow.com/questions/4643647/fast-prime-factorization-module
- uint64_t pi = (n >= 2) + (n >= 3);
- if (n < 5) return pi;
- n += 1;
- uint64_t correction = n % 6 > 1;
- uint64_t wheels[6] = { n, n - 1, n + 4, n + 3, n + 2, n + 1 };
- uint64_t limit = wheels[n % 6];
- boost::dynamic_bitset<> sieve(limit / 3);
- sieve.set();
- sieve[0] = false;
- for (uint64_t i = 0, upper = uint64_t(std::sqrt(limit))/3; i <= upper; ++i) {
- if (sieve[i]) {
- uint64_t k = (3*i + 1) | 1;
- for (uint64_t j = (k*k) / 3; j < limit/3; j += 2*k) sieve[j] = false;
- for (uint64_t j = (k*k + 4*k - 2*k*(i & 1))/3; j < limit/3; j += 2*k) sieve[j] = false;
- }
- }
- pi += sieve.count();
- for (uint64_t i = limit / 3 - correction; i < limit / 3; ++i) pi -= sieve[i];
- return pi;
- }
- int main(int argc, char** argv) {
- if (argc <= 1) {
- std::cout << "Usage: " << argv[0] << " nn";
- return 0;
- }
- std::cout << num_primes(std::stoi(argv[1])) << "n";
- return 0;
- }
- real 0m22.760s
- user 0m22.704s
- sys 0m0.080s
- real 0m22.854s
- user 0m22.800s
- sys 0m0.077s
- real 0m22.742s
- user 0m22.700s
- sys 0m0.066s
- real 0m22.484s
- user 0m22.450s
- sys 0m0.059s
- real 0m22.653s
- user 0m22.597s
- sys 0m0.080s
- real 0m22.665s
- user 0m22.602s
- sys 0m0.088s
- real 0m22.528s
- user 0m22.489s
- sys 0m0.062s
- real 0m22.510s
- user 0m22.474s
- sys 0m0.060s
- real 0m22.819s
- user 0m22.759s
- sys 0m0.084s
- real 0m22.488s
- user 0m22.459s
- sys 0m0.053s
- 10^4 (10,000) => 0.001 seconds
- 10^5 (100,000) => 0.003 seconds
- 10^6 (1,000,000) => 0.009 seconds
- 10^7 (10,000,000) => 0.074 seconds
- 10^8 (100,000,000) => 1.193 seconds
- 10^9 (1,000,000,000) => 14.415 seconds
- #include <cstdint>
- #include <vector>
- #include <iostream>
- #include <limits>
- #include <cmath>
- #include <array>
- // uses posix ffsll
- #include <string.h>
- #include <algorithm>
- #include <thread>
- constexpr uint64_t wheel_width = 2;
- constexpr uint64_t buf_size = 1<<(10+6);
- constexpr uint64_t dtype_width = 6;
- constexpr uint64_t dtype_mask = 63;
- constexpr uint64_t buf_len = ((buf_size*wheel_width)>>dtype_width);
- constexpr uint64_t seg_len = 6*buf_size;
- constexpr uint64_t limit_i_max = 0xfffffffe00000001ULL;
- typedef std::vector<uint64_t> buf_type;
- void mark_composite(buf_type& buf, uint64_t prime,
- std::array<uint64_t, 2>& poff,
- uint64_t seg_start, uint64_t max_j)
- {
- const auto p = 2*prime;
- for(uint64_t k = 0; k < wheel_width; ++k)
- {
- for(uint64_t j = 2*poff[k]+(k==0); j < max_j; j += p)
- {
- buf[(j-seg_start)>>dtype_width] |= 1ULL << (j & dtype_mask);
- poff[k] += prime;
- }
- }
- }
- struct prime_counter
- {
- buf_type buf;
- uint64_t n;
- uint64_t seg_a, seg_b;
- uint64_t nj;
- uint64_t store_max;
- uint64_t& store_res;
- prime_counter(uint64_t n, uint64_t seg_a, uint64_t seg_b, uint64_t nj, uint64_t store_max,
- uint64_t& store_res) :
- buf(buf_len), n(n), nj(nj), seg_a(seg_a), seg_b(seg_b),
- store_max(store_max), store_res(store_res)
- {}
- prime_counter(const prime_counter&) = default;
- prime_counter(prime_counter&&) = default;
- prime_counter& operator =(const prime_counter&) = default;
- prime_counter& operator =(prime_counter&&) = default;
- void operator()(uint64_t nsmall_segs,
- const std::vector<uint64_t>& primes,
- std::vector<std::array<uint64_t, 2> > poffs)
- {
- uint64_t res = 0;
- // no new prime added portion
- uint64_t seg_start = buf_size*wheel_width*seg_a;
- uint64_t seg_min = seg_len*seg_a+5;
- if(seg_a > nsmall_segs)
- {
- uint64_t max_j = buf_size*wheel_width*nsmall_segs+(seg_a-nsmall_segs)*(buf_len<<dtype_width);
- for(size_t k = 0; k < wheel_width; ++k)
- {
- for(uint64_t i = 0; i < poffs.size() && max_j >= (2*poffs[i][k]+(k==0)); ++i)
- {
- // adjust poffs
- // TODO: might be a more efficient way
- auto w = (max_j-(2*poffs[i][k]+(k==0)));
- poffs[i][k] += primes[i]*(w/(2*primes[i]));
- if(w % (2*primes[i]) != 0)
- {
- poffs[i][k]+=primes[i];// += primes[i]*(w/(2*primes[i])+1);
- }
- /*else
- {
- }*/
- }
- }
- }
- for(uint64_t seg = seg_a; seg < seg_b; ++seg)
- {
- std::fill(buf.begin(), buf.end(), 0);
- const uint64_t limit_i = std::min<uint64_t>((((seg_len+seg_min) >= limit_i_max) ?
- std::numeric_limits<uint32_t>::max() :
- ceil(sqrt(seg_len+seg_min))),
- store_max);
- uint64_t max_j = std::min(seg_start+(buf_len<<dtype_width), nj);
- for(uint64_t i = 0; i < primes.size() && primes[i] <= limit_i; ++i)
- {
- mark_composite(buf, primes[i], poffs[i], seg_start, max_j);
- }
- // sieve
- uint64_t val;
- const uint64_t stop = std::min(seg_min+seg_len, n);
- for(uint64_t i = ffsll(~(buf[0]))-((~buf[0]) != 0)+64*((~buf[0]) == 0);
- (val = 6ULL*(i>>1)+seg_min+2ULL*(i&1ULL)) < stop;)
- {
- if(!(buf[i>>dtype_width] & (1ULL << (i & dtype_mask))))
- {
- ++res;
- ++i;
- }
- else
- {
- uint64_t mask = buf[i>>dtype_width]>>(i&dtype_mask);
- const int64_t inc = ffsll(~mask)-((~mask) != 0)+64*((~mask) == 0);
- i += inc;
- }
- }
- seg_min += seg_len;
- seg_start += buf_size*wheel_width;
- }
- store_res = res;
- }
- };
- uint64_t num_primes(uint64_t n)
- {
- uint64_t res = (n >= 2) + (n >= 3);
- if(n >= 5)
- {
- buf_type buf(buf_len);
- // compute and store primes < sqrt(n)
- const uint64_t store_max = ceil(sqrt(n));
- // only primes >= 5
- std::vector<uint64_t> primes;
- std::vector<std::array<uint64_t, 2> > poffs;
- primes.reserve(ceil(1.25506*store_max/log(store_max)));
- poffs.reserve(ceil(1.25506*store_max/log(store_max)));
- uint64_t seg_start = 0;
- uint64_t seg_min = 5;
- const uint64_t num_segs = 1+(n-seg_min)/seg_len;
- const uint64_t nj = (n-seg_min)/3+1;
- // compute how many small segments there are
- const uint64_t nsmall_segs = 1+(store_max-seg_min)/seg_len;
- for(uint64_t seg = 0; seg < nsmall_segs; ++seg)
- {
- std::fill(buf.begin(), buf.end(), 0);
- // mark off small primes
- const uint64_t limit_i = std::min<uint64_t>((((seg_len+seg_min) >= limit_i_max) ?
- std::numeric_limits<uint32_t>::max() :
- ceil(sqrt(seg_len+seg_min))),
- store_max);
- uint64_t max_j = std::min(seg_start+(buf_len<<dtype_width), nj);
- for(uint64_t i = 0; i < primes.size() && primes[i] <= limit_i; ++i)
- {
- mark_composite(buf, primes[i], poffs[i], seg_start, max_j);
- }
- // sieve
- uint64_t val;
- const uint64_t stop = std::min(seg_min+seg_len, n);
- for(uint64_t i = ffsll(~(buf[0]))-((~buf[0]) != 0)+64*((~buf[0]) == 0);
- (val = 6ULL*(i>>1)+seg_min+2ULL*(i&1ULL)) < stop;)
- {
- if(!(buf[i>>dtype_width] & (1ULL << (i & dtype_mask))))
- {
- if(val <= store_max)
- {
- // add prime and poffs
- primes.push_back(val);
- poffs.emplace_back();
- poffs.back()[0] = (val*val-1)/6-1;
- if(i&1)
- {
- // 6n+1 prime
- poffs.back()[1] = (val*val+4*val-5)/6;
- }
- else
- {
- // 6n+5 prime
- poffs.back()[1] = (val*val+2*val-5)/6;
- }
- // mark-off multiples
- mark_composite(buf, val, poffs.back(), seg_start, max_j);
- }
- ++res;
- ++i;
- }
- else
- {
- uint64_t mask = buf[i>>dtype_width]>>(i&dtype_mask);
- const int64_t inc = ffsll(~mask)-((~mask) != 0)+64*((~mask) == 0);
- i += inc;
- }
- }
- seg_min += seg_len;
- seg_start += buf_size*wheel_width;
- }
- // multi-threaded sieving for remaining segments
- std::vector<std::thread> workers;
- auto num_workers = std::min<uint64_t>(num_segs-nsmall_segs, std::thread::hardware_concurrency());
- std::vector<uint64_t> store_reses(num_workers);
- workers.reserve(num_workers);
- auto num_segs_pw = ceil((num_segs-nsmall_segs)/static_cast<double>(num_workers));
- for(size_t i = 0; i < num_workers; ++i)
- {
- workers.emplace_back(prime_counter(n, nsmall_segs+i*num_segs_pw,
- std::min<uint64_t>(nsmall_segs+(i+1)*num_segs_pw,
- num_segs),
- nj, store_max, store_reses[i]),
- nsmall_segs, primes, poffs);
- }
- for(size_t i = 0; i < num_workers; ++i)
- {
- workers[i].join();
- res += store_reses[i];
- }
- }
- return res;
- }
- int main(int argc, char** argv)
- {
- if(argc <= 1)
- {
- std::cout << "usage: " << argv[0] << " nn";
- return -1;
- }
- std::cout << num_primes(std::stoll(argv[1])) << 'n';
- }
- g++ -std=c++11 -o sieve_mt -O3 -march=native -pthread sieve_mt.cpp
- real 4m7.215s
- user 23m54.086s
- sys 0m1.239s
- #include <stdlib.h>
- #include <stdio.h>
- #include <string.h>
- unsigned char p[2000000001];
- int main(int argc, char **argv)
- {
- unsigned int n, c, i, j;
- n = atoi(argv[1]);
- memset(p, 1, n + 1);
- p[1] = p[0] = 0;
- for (i = 2, c = 0; i <= n; i++)
- {
- if (p[i])
- {
- c++;
- for (j = i + i; j <= n; j += i)
- p[j] = 0;
- }
- }
- printf("%d: %dn", n, c);
- return 0;
- }
- real 2m42.657s
- user 2m42.065s
- sys 0m0.757s
- real 2m42.947s
- user 2m42.400s
- sys 0m0.708s
- real 2m42.827s
- user 2m42.282s
- sys 0m0.703s
- real 2m42.800s
- user 2m42.300s
- sys 0m0.665s
- real 2m42.562s
- user 2m42.050s
- sys 0m0.675s
- real 2m42.788s
- user 2m42.192s
- sys 0m0.756s
- real 2m42.631s
- user 2m42.074s
- sys 0m0.720s
- real 2m42.658s
- user 2m42.115s
- sys 0m0.707s
- real 2m42.710s
- user 2m42.219s
- sys 0m0.657s
- real 2m42.674s
- user 2m42.110s
- sys 0m0.730s
- real 1m21.227s
- user 1m20.755s
- sys 0m0.576s
- real 1m20.944s
- user 1m20.426s
- sys 0m0.640s
- real 1m21.052s
- user 1m20.581s
- sys 0m0.573s
- real 1m21.328s
- user 1m20.862s
- sys 0m0.570s
- real 1m21.253s
- user 1m20.780s
- sys 0m0.588s
- real 1m20.925s
- user 1m20.460s
- sys 0m0.576s
- real 1m21.011s
- user 1m20.512s
- sys 0m0.601s
- real 1m21.011s
- user 1m20.550s
- sys 0m0.564s
- real 1m20.875s
- user 1m20.409s
- sys 0m0.569s
- real 1m21.703s
- user 1m21.088s
- sys 0m0.701s
- public class PrimeCounter
- {
- public static final String START_CODE="=",
- TEST_FORMAT="Input = %d , Output = %d , calculated in %f seconds%n",
- PROMPT="Enter numbers to compute pi(x) for (Type ""+START_CODE+"" to start):%n",
- WAIT="Calculating, please wait...%n",
- WARNING="Probably won't work with values close to or more than 2^31%n",
- TOTAL_OUTPUT_FORMAT="Total time for all inputs is %f seconds%n";
- public static final int NUM_THREADS=16,LOW_LIM=1,HIGH_LIM=1<<28;
- private static final Object LOCK=new Lock();
- private static final class Lock{}
- /**
- * Generates and counts primes using an optimized but naive iterative algorithm.
- * Uses MultiThreading for arguments above LOW_LIM
- * @param MAX : argument x for pi(x), the limit to which to generate numbers.
- */
- public static long primeCount(long MAX){
- long ctr=1;
- if(MAX<1<<7){
- for(long i=3;i<=MAX;i+=2){
- if(isPrime(i))++ctr;
- }
- }else{
- long[] counts=new long[NUM_THREADS];
- for(int i=0;i<NUM_THREADS;++i){
- counts[i]=-1;
- }
- long range=Math.round((double)MAX/NUM_THREADS);
- for(int i=0;i<NUM_THREADS;++i){
- long start=(i==0)?3:i*range+1,end=(i==NUM_THREADS-1)?MAX:(i+1)*range;
- final int idx=i;
- new Thread(new Runnable(){
- public void run(){
- for(long j=start;j<=end;j+=2){
- if(isPrime(j))++counts[idx];
- }
- }
- }).start();
- }
- synchronized(LOCK){
- while(!completed(counts)){
- try{
- LOCK.wait(300);}catch(InterruptedException ie){}
- }
- LOCK.notifyAll();
- }
- for(long count:counts){
- ctr+=count;
- }
- ctr+=NUM_THREADS;
- }
- return ctr;
- }
- /**
- * Checks for completion of threads
- * @param array : The array containing the completion data
- */
- private static boolean completed(long[] array){
- for(long i:array){
- if(i<0)return false;
- }return true;
- }
- /**
- * Checks if the parameter is prime or not.
- * 2,3,5,7 are hardcoded as factors.
- * @param n : the number to check for primality
- */
- private static boolean isPrime(long n){
- if(n==2||n==3||n==5||n==7)return true;
- else if(n%2==0||n%3==0||n%5==0||n%7==0)return false;
- else{
- for(long i=11;i<n;i+=2){
- if(n%i==0)return false;
- }
- return true;
- }
- }
- /**
- * Calculates primes using the atandard Sieve of Eratosthenes.
- * Uses 2,3,5,7 wheel factorization for elimination (hardcoded for performance reasons)
- * @param MAX : argument x for pi(x)
- * Will delegate to <code>primeCount(long)</code> for MAX<LOW_LIM and to <code>bitPrimeSieve(long)</code>
- * for MAX>HIGH_LIM, for performance reasons.
- */
- public static long primeSieve(long MAX){
- if(MAX<=1)return 0;
- else if(LOW_LIM>0&&MAX<LOW_LIM){return primeCount(MAX);}
- else if(HIGH_LIM>0&&MAX>HIGH_LIM){return bitPrimeSieve(MAX);}
- int n=(int)MAX;
- int sn=(int)Math.sqrt(n),ctr=2;
- if(sn%2==0)--sn;
- boolean[]ps=new boolean[n+1];
- for(int i=2;i<=n;++i){
- if(i==2||i==3||i==5||i==7)ps[i]=true;
- else if(i%2!=0&&i%3!=0&&i%5!=0&&i%7!=0)ps[i]=true;
- else ++ctr;
- }
- for(int i=(n>10)?11:3;i<=sn;i+=2){
- if(ps[i]){
- for(int j=i*i;j<=n;j+=i){
- if(ps[j]){ ps[j]=false;++ctr;}
- }
- }
- }
- return (n+1-ctr);
- }
- /**
- * Calculates primes using bitmasked Sieve of Eratosthenes.
- * @param MAX : argument x for pi(x)
- */
- public static long bitPrimeSieve(long MAX) {
- long SQRT_MAX = (long) Math.sqrt(MAX);
- if(SQRT_MAX%2==0)--SQRT_MAX;
- int MEMORY_SIZE = (int) ((MAX+1) >> 4);
- byte[] array = new byte[MEMORY_SIZE];
- for (long i = 3; i <= SQRT_MAX; i += 2) {
- if ((array[(int) (i >> 4)] & (byte) (1 << ((i >> 1) & 7))) == 0) {
- for(long j=i*i;j<=MAX;j+=i<<1) {
- if((array[(int) (j >> 4)] & (byte) (1 << ((j >> 1) & 7))) == 0){
- array[(int) (j >> 4)] |= (byte) (1 << ((j >> 1) & 7));
- }
- }
- }
- }
- long pi = 1;
- for (long i = 3; i <= MAX; i += 2) {
- if ((array[(int) (i >> 4)] & (byte) (1 << ((i >> 1) & 7))) == 0) {
- ++pi;
- }
- }
- return pi;
- }
- /**
- * Private testing and timer function
- * @param MAX : input to be passed on to <code>primeSieve(long)</code>
- */
- private static long sieveTest(long MAX){
- long start=System.nanoTime();
- long ps=primeSieve(MAX);
- long end=System.nanoTime();
- System.out.format(TEST_FORMAT,MAX,ps,((end-start)/1E9));
- return end-start;
- }
- /**
- * Main method: accepts user input and shows total execution time taken
- * @param args : The command-line arguments
- */
- public static void main(String[]args){
- double total_time=0;
- java.util.Scanner sc=new java.util.Scanner(System.in);
- java.util.ArrayList<Long> numbers=new java.util.ArrayList<>();
- System.out.format(PROMPT+WARNING);
- String line=sc.nextLine();
- while(!line.equals(START_CODE)/*sc.hasNextLine()&&Character.isDigit(line.charAt(0))*/){
- numbers.add(Long.valueOf(line));
- line=sc.nextLine();
- }
- System.out.format(WAIT);
- for(long num:numbers){
- total_time+=sieveTest(num);
- }
- System.out.format(TOTAL_OUTPUT_FORMAT,total_time/1e9);
- }
- }
- javac PrimeCounter.java
- java PrimeCounter
- Enter numbers to compute pi(x) for (Type "=" to start):
- Probably won't work with values close to or more than 2^31
- 41500
- 24850000
- 40550000
- 99820000
- 660000000
- 1240000000
- 1337000000
- 1907000000
- =
- Calculating, please wait...
- Input = 41500 , Output = 4339 , calculated in 0.002712 seconds
- Input = 24850000 , Output = 1557132 , calculated in 0.304792 seconds
- Input = 40550000 , Output = 2465109 , calculated in 0.523999 seconds
- Input = 99820000 , Output = 5751639 , calculated in 1.326542 seconds
- Input = 660000000 , Output = 34286170 , calculated in 4.750049 seconds
- Input = 1240000000 , Output = 62366021 , calculated in 9.160406 seconds
- Input = 1337000000 , Output = 66990613 , calculated in 9.989093 seconds
- Input = 1907000000 , Output = 93875448 , calculated in 14.832107 seconds
- Total time for all inputs is 40.889700 seconds
- import sys
- sys.setrecursionlimit(sys.maxsize)
- n = int(sys.argv[-1])
- if n < 4:
- print(0 if n < 2 else n-1)
- exit()
- p = [0, 0] + [True] * n
- i = 0
- while i < pow(n, 0.5):
- if p[i]:
- j = pow(i, 2)
- while j < n:
- p[j] = False
- j += i
- i += 1
- print(sum(p) - 2)
- $ time python3 test.py 90
- 24
- real 0m0.045s
- user 0m0.031s
- sys 0m0.010s
- #include <cstdint>
- #include <vector>
- #include <iostream>
- #include <limits>
- #include <cmath>
- #include <array>
- // uses posix ffsll
- #include <string.h>
- #include <algorithm>
- constexpr uint64_t wheel_width = 2;
- constexpr uint64_t buf_size = 1<<(10+6);
- constexpr uint64_t dtype_width = 6;
- constexpr uint64_t dtype_mask = 63;
- constexpr uint64_t buf_len = ((buf_size*wheel_width)>>dtype_width);
- typedef std::vector<uint64_t> buf_type;
- void mark_composite(buf_type& buf, uint64_t prime,
- std::array<uint64_t, 2>& poff,
- uint64_t seg_start, uint64_t max_j)
- {
- const auto p = 2*prime;
- for(uint64_t k = 0; k < wheel_width; ++k)
- {
- for(uint64_t j = 2*poff[k]+(k==0); j < max_j; j += p)
- {
- buf[(j-seg_start)>>dtype_width] |= 1ULL << (j & dtype_mask);
- poff[k] += prime;
- }
- }
- }
- uint64_t num_primes(uint64_t n)
- {
- uint64_t res = (n >= 2) + (n >= 3);
- if(n >= 5)
- {
- buf_type buf(buf_len);
- // compute and store primes < sqrt(n)
- const uint64_t store_max = ceil(sqrt(n));
- // only primes >= 5
- std::vector<uint64_t> primes; // 5,7,11
- std::vector<std::array<uint64_t, 2> > poffs;// {{3,0},{0,5},{8,1}};
- primes.reserve(ceil(1.25506*store_max/log(store_max)));
- poffs.reserve(ceil(1.25506*store_max/log(store_max)));
- uint64_t seg_start = 0;
- uint64_t seg_min = 5;
- constexpr uint64_t seg_len = 6*buf_size;///wheel_width;
- constexpr uint64_t limit_i_max = 0xfffffffe00000001ULL;
- const uint64_t num_segs = 1+(n-seg_min)/seg_len;
- const uint64_t nj = (n-seg_min)/3+1;
- for(uint64_t seg = 0; seg < num_segs; ++seg)
- {
- std::fill(buf.begin(), buf.end(), 0);
- // mark off small primes
- const uint64_t limit_i = std::min<uint64_t>((((seg_len+seg_min) >= limit_i_max) ?
- std::numeric_limits<uint32_t>::max() :
- ceil(sqrt(seg_len+seg_min))),
- store_max);
- uint64_t max_j = std::min(seg_start+(buf_len<<dtype_width), nj);
- for(uint64_t i = 0; i < primes.size() && primes[i] <= limit_i; ++i)
- {
- mark_composite(buf, primes[i], poffs[i], seg_start, max_j);
- }
- // sieve
- uint64_t val;
- const uint64_t stop = std::min(seg_min+seg_len, n);
- for(uint64_t i = ffsll(~(buf[0]))-((~buf[0]) != 0)+64*((~buf[0]) == 0);
- (val = 6ULL*(i>>1)+seg_min+2ULL*(i&1ULL)) < stop;)
- {
- if(!(buf[i>>dtype_width] & (1ULL << (i & dtype_mask))))
- {
- if(val <= store_max)
- {
- // add prime and poffs
- primes.push_back(val);
- poffs.emplace_back();
- poffs.back()[0] = (val*val-1)/6-1;
- if(i&1)
- {
- // 6n+1 prime
- poffs.back()[1] = (val*val+4*val-5)/6;
- }
- else
- {
- // 6n+5 prime
- poffs.back()[1] = (val*val+2*val-5)/6;
- }
- // mark-off multiples
- mark_composite(buf, val, poffs.back(), seg_start, max_j);
- }
- ++res;
- ++i;
- }
- else
- {
- uint64_t mask = buf[i>>dtype_width]>>(i&dtype_mask);
- const int64_t inc = ffsll(~mask)-((~mask) != 0)+64*((~mask) == 0);
- i += inc;
- }
- }
- seg_min += seg_len;
- seg_start += buf_size*wheel_width;
- }
- }
- return res;
- }
- int main(int argc, char** argv)
- {
- if(argc <= 1)
- {
- std::cout << "usage: " << argv[0] << " nn";
- return -1;
- }
- std::cout << num_primes(std::stoll(argv[1])) << 'n';
- }
- g++ -std=c++11 -o sieve -O3 -march=native sieve.cpp
- 41500
- 4339
- real 0m0.001s
- user 0m0.000s
- sys 0m0.000s
- 24850000
- 1557132
- real 0m0.036s
- user 0m0.032s
- sys 0m0.000s
- 40550000
- 2465109
- real 0m0.056s
- user 0m0.052s
- sys 0m0.000s
- 99820000
- 5751639
- real 0m0.149s
- user 0m0.144s
- sys 0m0.000s
- 660000000
- 34286170
- real 0m0.795s
- user 0m0.788s
- sys 0m0.000s
- 1240000000
- 62366021
- real 0m1.468s
- user 0m1.464s
- sys 0m0.000s
- 1337000000
- 66990613
- real 0m1.583s
- user 0m1.576s
- sys 0m0.004s
- 1907000000
- 93875448
- real 0m2.363s
- user 0m2.356s
- sys 0m0.000s
- real 0m9.415s
- user 0m9.414s
- sys 0m0.014s
- real 0m9.315s
- user 0m9.315s
- sys 0m0.013s
- real 0m9.307s
- user 0m9.309s
- sys 0m0.012s
- real 0m9.333s
- user 0m9.330s
- sys 0m0.017s
- real 0m9.288s
- user 0m9.289s
- sys 0m0.012s
- real 0m9.319s
- user 0m9.318s
- sys 0m0.015s
- real 0m9.285s
- user 0m9.284s
- sys 0m0.015s
- real 0m9.342s
- user 0m9.342s
- sys 0m0.014s
- real 0m9.305s
- user 0m9.305s
- sys 0m0.014s
- real 0m9.312s
- user 0m9.313s
- sys 0m0.012s
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement