Compiler Explorer

Source code

#include <cstdint>
#include <vector>
#include <iostream>
#include <iomanip>
#include <fstream>
#include <stdexcept>
#include <optional>
#include <future>
#include <string>
#include <chrono>
#include <ctime>
#include <sstream>
#include <algorithm>
#include <mutex>
#include <atomic>
#include <thread>

#define ASSERT_MSG(cond, msg) { if (!(cond)) throw std::runtime_error("Assertion (" #cond ") failed at line " + std::to_string(__LINE__) + "! Msg: '" + std::string(msg) + "'."); }
#define ASSERT(cond) ASSERT_MSG(cond, "")
#define COUT(code) { std::lock_guard<std::recursive_mutex> lock(g_cout_mux); std::cout << TimeUtcStr() << ": " code; std::cout << std::flush; }
#define LN { std::cout << "LN " << std::setw(5) << __LINE__ << "  " << __FUNCTION__ << "  " << __FILE__ << std::endl << std::flush; }

static std::recursive_mutex g_cout_mux;

using u8 = uint8_t;
using u32 = uint32_t;
using u64 = uint64_t;

static std::string const out_fname = "output.txt";
u64 constexpr
    out_file_buf_size = 16ULL << 20, // This much bytes is used for output file buffering
    base_range_bits = 63, // Start sieving from 2^base_range_bits + 2^small_range_bits * small_irange_start
    small_range_bits = 24, // Small range bits sieved in single thread, in this case [2^63 + 2^24 * 5, 2^63 + 2^24 * (5 + 1)).
    small_irange_start = 0; // Which small range index should start from, all previous ranges should be sieved already

static double TimeSec() {
    static auto const gtb = std::chrono::high_resolution_clock::now();
    return std::chrono::duration_cast<std::chrono::duration<double>>(std::chrono::high_resolution_clock::now() - gtb)
        .count();
}

static size_t NThreads() {
    return std::thread::hardware_concurrency();
}

static std::string TimeUtcStr() {
    auto now = std::chrono::system_clock::now();
    auto seconds_since_epoch = std::chrono::duration_cast<std::chrono::seconds>(now.time_since_epoch());
    auto fractional_seconds = std::chrono::duration_cast<std::chrono::microseconds>(
        now.time_since_epoch() - seconds_since_epoch);
    std::time_t now_as_time_t = std::chrono::system_clock::to_time_t(now);
    std::tm utc_tm = *std::gmtime(&now_as_time_t);
    char buffer[80];
    std::strftime(buffer, sizeof(buffer), "%Y.%m.%d %H:%M:%S", &utc_tm);
    std::stringstream ss;
    ss << buffer << '.' << std::setfill('0') << std::setw(3) << fractional_seconds.count() / 1000;
    return ss.str();
}

template <typename T>
std::vector<T> GenPrimesSieve_Eratosthenes(uint64_t const end) {
    // https://en.wikipedia.org/wiki/Sieve_of_Eratosthenes
    if (end <= 2)
        return {};
    uint64_t const cnt = end >> 1;
    std::vector<uint8_t> composites((cnt + 7) / 8);
    auto Get = [&](uint64_t const i){
        return bool((composites[i / 8] >> (i % 8)) & 1);
    };
    auto Set = [&](uint64_t const i){
        composites[i / 8] |= uint8_t(1) << (i % 8);
    };
    std::vector<T> primes = {2};
    uint64_t i = 0, ops = 0;
    for (i = 1; i < cnt; ++i) {
        if (Get(i))
            continue;
        uint64_t const
            p = 2 * i + 1,
            start = (p * p) >> 1;
        primes.push_back(p);
        if (start >= cnt)
            break;
        for (uint64_t j = start; j < cnt; j += p)
            Set(j);
    }
    for (i = i + 1; i < cnt; ++i)
        if (!Get(i))
            primes.push_back(2 * i + 1);
    primes.shrink_to_fit();
    return primes;
}

static std::vector<u64> PrimeSieve(u64 begin, u64 end, std::vector<u32> const & small_primes) {
    if (begin % 2 == 0) ++begin;
    if (end % 2 == 1) --end;
    if (end <= begin) return {};
    
    u64 const cnt = (end - begin + 1) >> 1;
    std::vector<uint8_t> composites((cnt + 7) / 8); composites.shrink_to_fit();
    auto Get = [&](u64 i){ return bool((composites[i / 8] >> (i % 8)) & 1); };
    auto Set = [&](u64 i){ composites[i / 8] |= u8(1) << (i % 8); };
    
    for (u64 const p: small_primes) {
        u64 const
            t = begin - begin % (p * 2) + p,
            start = t < begin ? t + p * 2 : t;
        for (u64 i = (start - begin) >> 1; i < cnt; i += p)
            Set(i);
        if (u64(p) * p >= end) {
            ASSERT_MSG(p < small_primes.back(),
                "To few small primes! Increase Eratosthenes sieving range.");
            break;
        }
    }
    
    std::vector<u64> primes;
    for (u64 i = 0; i < cnt; ++i)
        if (!Get(i))
            primes.push_back(begin + i * 2);
    primes.shrink_to_fit();
    
    return primes;
}

static auto FactorSieve(u64 const begin, u64 const end, std::vector<u32> const & small_primes,
        u64 const fbegin, std::vector<bool> const & filter) {
    std::vector<std::vector<u32>> lo_factors(end - begin); lo_factors.shrink_to_fit();
    std::vector<std::optional<u64>> hi_factors(end - begin); hi_factors.shrink_to_fit();
    std::vector<u64> prods(end - begin, 1); prods.shrink_to_fit();
    
    for (auto const p: small_primes) {
        u64 mp = p;
        while (true) {
            u64 const
                t = begin - begin % mp,
                start = t < begin ? t + mp : t;
            if (start >= t && end + mp >= end) // Check add-overflow
                for (u64 i = start; i < end; i += mp)
                    if (filter[i * 2 + 1 - fbegin]) {
                        lo_factors[i - begin].push_back(p);
                        prods[i - begin] *= p;
                    }
            if ((mp * p) / p != mp)
                break;
            mp *= p;
        }
        if (u64(p) * p >= end) {
            ASSERT_MSG(p < small_primes.back(),
                "To few small primes! Increase Eratosthenes sieving range.");
            break;
        }
    }
    
    for (u64 i = begin; i < end; ++i)
        if (filter[i * 2 + 1 - fbegin]) {
            std::sort(lo_factors[i - begin].begin(), lo_factors[i - begin].end());
            lo_factors[i - begin].shrink_to_fit();
        }
    
    for (u64 i = begin; i < end; ++i)
        if (filter[i * 2 + 1 - fbegin] && prods[i - begin] < i) {
            u64 const f = i / prods[i - begin];
            ASSERT(f * prods[i - begin] == i);
            hi_factors[i - begin] = f;
        }
    
    return std::make_pair(std::move(lo_factors), std::move(hi_factors));
}

static std::atomic<u64> cnt_outputted_ranges = small_irange_start;
static std::mutex output_mutex;

static void ProcessRange(size_t const base_bits, size_t const range_bits,
        u64 const irange, std::vector<u32> const & small_primes) {
    u64 const
        begin = (1ULL << base_bits) + (irange + 0) * (1ULL << range_bits),
        end   = (1ULL << base_bits) + (irange + 1) * (1ULL << range_bits);
    
    auto time_primes = TimeSec();
    auto const ps = PrimeSieve(begin, end, small_primes);
    time_primes = TimeSec() - time_primes;
    
    std::vector<bool> filter(end - begin);
    for (auto p: ps)
        filter.at(p - begin) = true;
    filter.shrink_to_fit();
    u64 const hbegin = begin / 2,
        hend = ((end - 1) / 2 * 2 + 1 > end - 1 ? (end - 2) / 2 : (end - 1) / 2) + 1;
    
    auto time_factor = TimeSec();
    auto const [flo, fhi] = FactorSieve(hbegin, hend, small_primes, begin, filter);
    time_factor = TimeSec() - time_factor;
    
    u64 mem_used = ps.capacity() * sizeof(ps[0]) + filter.capacity() * 1 +
        flo.capacity() * sizeof(flo[0]) + fhi.capacity() * sizeof(fhi[0]);
    for (u64 i = 0; i < flo.size(); ++i)
        mem_used += flo[i].capacity() * sizeof(flo[i][0]);
    
    std::ostringstream file;
    
    u64 hprimes_cnt = 0;
    for (u64 i = begin; i < end; ++i) {
        if (!filter.at(i - begin))
            continue; 
        u64 const q = (i - 1) / 2;
        file << "Range [2^" << base_bits << " + 2^" << range_bits << " * " << irange
            << ", 2^" << base_bits << " + 2^" << range_bits << " * " << (irange + 1)
            << "), prime p = " << i << ", factors of (p - 1) / 2 = " << q << ": ";
        u64 prod = 1;
        auto const & v = flo.at(q - hbegin);
        auto const & h = fhi.at(q - hbegin);
        size_t c = 0;
        for (size_t j = 0; j < v.size(); ++j) {
            ++c;
            prod *= v[j];
            file << v[j] << (j + 1 >= v.size() && !h ? "" : ", ");
        }
        if (h) {
            ++c;
            prod *= *h;
            file << *h;
        }
        file << std::endl;
        if (c == 1)
            ++hprimes_cnt;
        ASSERT_MSG(prod == q, "For p = " + std::to_string(i) + " error in factoring (p - 1) / 2 = " +
                std::to_string(q) + ", product of factors " + std::to_string(prod) + ", factors: " + [&]{
            std::ostringstream ss;
            for (auto e: v)
                ss << e << ", ";
            if (h)
                ss << *h;
            return ss.str();
        }());
        ASSERT(c > 0);
    }
    
    while (true) {
        std::unique_lock<std::mutex> file_lock(output_mutex);
        
        if (irange > cnt_outputted_ranges) {
            file_lock.unlock();
            std::this_thread::yield();
        } else {
            ASSERT_MSG(irange == cnt_outputted_ranges, "irange " + std::to_string(irange) +
                ", cnt_outputted_ranges " + std::to_string(cnt_outputted_ranges));
            std::ofstream ofile(out_fname, std::ios::app);
            ASSERT_MSG(ofile.is_open(), "Can't open file '" + out_fname + "' for write-append!");
            //ofile.rdbuf()->pubsetbuf(nullptr, out_file_buf_size);
            ofile << file.str();
            ofile.close();
            ASSERT_MSG(ofile.good(), "Error closing file '" + out_fname + "'!");
            cnt_outputted_ranges = irange + 1;
            break;
        }
    }
    
    COUT(<< "Processed range [2^" << base_bits << " + 2^" << range_bits << " * " << std::setw(3) << irange
        << ", 2^" << base_bits << " + 2^" << range_bits << " * " << std::setw(3) << (irange + 1)
        << "), prime sieve time " << std::fixed << std::setprecision(1) << std::setw(4) << time_primes
        << " sec, count primes " << ps.size() << ", half factor sieve time " << std::setw(4) << time_factor
        << " sec, count half primes (Sophie Germain) " << hprimes_cnt << ", range mem usage "
        << 1.0 * mem_used / (1 << 20) << " MiB" << std::endl);
}

static void Main() {
    COUT(<< "Creating table of all 32-bit primes! Wait around 1-2 minutes..." << std::endl);
    auto const small_primes = GenPrimesSieve_Eratosthenes<u32>(u32(-1));
    // http://www.umopit.ru/CompLab/primes32eng.htm
    // In total there are 203280221 primes that are 32-bit in size.
    ASSERT_MSG(small_primes.size() == 203280221,
        "small_primes.size() " + std::to_string(small_primes.size()));
    COUT(<< "Small primes count " << small_primes.size()
        << " and mem usage " << std::fixed << std::setprecision(1)
        << 1.0 * small_primes.capacity() * sizeof(small_primes[0]) / (1 << 20)
        << " MiB" << std::endl);
    // Erase file if starting from 0 range
    if (small_irange_start == 0) {
        COUT(<< "Erasing output file!" << std::endl);
        std::ofstream out_cleanup(out_fname);
        ASSERT(out_cleanup.is_open());
    }
    COUT(<< "Starting ranges processing from range index " << small_irange_start << std::endl);
    
    std::vector<std::future<void>> asyncs;
    for (u64 irange = small_irange_start;; ++irange) {
        while (asyncs.size() >= NThreads()) {
            for (size_t j = 0; j < asyncs.size(); ++j)
                if (asyncs[j].wait_for(std::chrono::microseconds(1)) == std::future_status::ready) {
                    asyncs[j].get();
                    asyncs.erase(asyncs.begin() + j);
                    break;
                }
            std::this_thread::yield();
        }
        asyncs.push_back(std::async(std::launch::async, [&, irange]{
            ProcessRange(base_range_bits, small_range_bits, irange, small_primes);
        }));
    }
    
    COUT(<< ": All done!" << std::endl);
}

int main() {
    try {
        Main();
        return 0;
    } catch (std::exception const & ex) {
        COUT(<< "Exception: " << ex.what() << std::endl);
        return -1;
    }
}