Sieve of Eratosthenes - prime factor

Question

I just wrote the following to utilize the sieve to find the largest prime factor of some natural number greater than 2.

The program builds and runs and works for small test values, but simply crashes for values larger than, say, 1000000.

I wrote this myself--and believe it will probably be extremely inefficient--after finding out what the sieve is about.

Can you please suggest improvements?

Thank you.

//LARGEST PRIME FACTOR w/SIEVE OF ERATHOSTHENES
#include <iostream>
#include <math.h>

using namespace std;

unsigned long long int numberToCheck=0;

void sieve(unsigned long long int checkNumber)
{
    checkNumber=numberToCheck;
    unsigned long long int root=(int)(sqrt(checkNumber));
    unsigned long long int primeFlagger[checkNumber+1];

    for(unsigned long long int i=0;i<=checkNumber;i++)
    {
        primeFlagger[i]=1;
    }

    primeFlagger[0]=0;
    primeFlagger[1]=0;

    for(unsigned long long int j=2;j<=checkNumber;j++)
    {
        if(primeFlagger[j]==1)
        {
            for(unsigned long long int k=2;(k*j)<=checkNumber;k++)
            {
                primeFlagger[j*k]=0;
            }
        }
    }

    for(unsigned long long int l=checkNumber;l>=0;l--)
    {
        if(primeFlagger[l]==1)
        {
            if(checkNumber%l==0)
            {
                cout<<l<<" is the largest prime factor"<<endl;
                break;
            }
        }
    }
}

int main()
{
    cout<<"Largest prime factor less then or equal to? "<<endl;
    cin>>numberToCheck;
    cout<<endl;
    cout<<"Retrieving largest prime factor..."<<endl;

    sieve(numberToCheck);



    return 0;
}

Answer 1

Other answer describes most of problems, but from my side I'll tell one most obvious problem regarding speed - you're checking primes up to number itself, but it will be much much faster if you check only till square root of number.

Instead of fixing your bugs, as they are described in other answer, I decided to implement my own code, that I think solves your task in best and most beautiful way. If you wanted only fixing bugs then please look at my answer as for educational purpose, maybe my code implementation would be of interest for you.

I implemented two functions, one finds largest factor using Sieve of Eratosthenes, another finds it using Trial Division algorithm. Trial Division iterates through all odd integers till square root of a number and checks if they divide given number, without using any sieve.

It appeared that on average regular Trial Division is faster about 10-20x compared to Sieve of Eratosthenes method. But Sieve of Eratosthenes is much faster if it is reused very many times. After comparing both timings I think that Sieve of Eratosthenes can be more profitable if re-used for at least 20 different numbers, if re-used for millions of numbers then it will be much much more profitable speedwise.

You can notice that in Sieve of Eratosthenes I did several optimizations. First to store bit vector signifying if a number is composite or prime, I used bit packing, used single bit per flag, for that I packed bits inside each uint8_t byte. Besides using less memory which is very necessary if you need to compute very many primes, it is also necessary to increase speed of a program because smaller memory has bigger chance to fit into L1-L3 caches of CPU.

Second optimization of Sieve of Eratosthenes that I did is keeping inside bit vector only bits of odd numbers, excluding even numbers. This further reduced occupied memory 2x times more, and increased speed of a program also 2x times more.

Inside code function LargestFactor_Sieve() solves your task using Sieve of Eratosthenes, and LargestFactor_Loop() solves your task using Trial Division.

Code run may take around 10 minutes only because it does timings of all bit sizes till 64 bits, and this 64 bit size takes most of time. See timings output after my code below.

Try it online!

#include <cstdint>
#include <vector>
#include <cmath>

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

template <typename T>
std::vector<T> GenPrimes_SieveOfEratosthenes(size_t end) {
    // https://en.wikipedia.org/wiki/Sieve_of_Eratosthenes
    if (end <= 2)
        return {};
    size_t const cnt = end >> 1;
    std::vector<u8> composites((cnt + 7) / 8);
    auto Get = [&](size_t i){ return bool((composites[i / 8] >> (i % 8)) & 1); };
    auto Set = [&](size_t i){ composites[i / 8] |= u8(1) << (i % 8); };
    std::vector<T> primes = {2};
    size_t i = 0;
    for (i = 1; i < cnt; ++i) {
        if (Get(i))
            continue;
        size_t const p = 2 * i + 1, start = (p * p) >> 1;
        primes.push_back(p);
        if (start >= cnt)
            break;
        for (size_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);
    return primes;
}

template <typename T>
T LargestFactor_Sieve(T n) {
    if (n <= 0)
        return 0;
    auto primes = GenPrimes_SieveOfEratosthenes<u32>(
        std::llround(std::sqrt(double(n)) - 0.5 + 1e-12) + 1);
    T last = 1;
    for (auto p: primes) {
        if (T(p) * p > n)
            break;
        while (n % p == 0) {
            n /= p;
            last = p;
        }
    }
    return n == 1 ? last : n;
}

template <typename T>
T LargestFactor_Loop(T n) {
    if (n <= 0)
        return 0;
    T last = 1;
    while ((n & 1) == 0) {
        n >>= 1;
        last = 2;
    }
    for (u32 d = 3; T(d) * d <= n; d += 2)
        while (n % d == 0) {
            n /= d;
            last = d;
        }
    return n == 1 ? last : n;
}

#include <chrono>
#include <iostream>
#include <iomanip>
#include <random>

inline double Time() {
    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();
}

template <typename T>
void Test0(size_t bits, size_t ntests) {
    std::mt19937_64 rng{std::random_device{}()};
    std::vector<u64> nums;
    for (size_t i = 0; i < ntests; ++i)
        nums.push_back((rng() | (bits <= 48 ? 0ULL : 1ULL)) >> (64 - bits));
    u64 sum0 = 0, sum1 = 0;
    double tim0 = Time(), tim1 = 0;
    for (auto n: nums)
        sum0 += LargestFactor_Sieve<T>(n);
    tim0 = Time() - tim0;
    tim1 = Time();
    for (auto n: nums)
        sum1 += LargestFactor_Loop<T>(n);
    tim1 = Time() - tim1;
    std::cout << "Bits " << std::setw(2) << bits << ", Sieve " << std::fixed
        << std::setprecision(6) << std::setw(10) << tim0 / ntests << " sec, Loop "
        << std::setw(10) << tim1 / ntests << " sec, LoopBoost "
        << std::setw(5) << std::setprecision(1) << tim0 / tim1 << "x"
        << (sum0 == sum1 ? "" : ", NON-EQUAL!!!") << std::endl;
}

void Test() {
    for (size_t bits = 4; bits <= 64; bits += 4)
        Test0<u64>(bits, bits <= 40 ? 256 : bits <= 48 ? 32 : bits <= 56 ? 16 : 8);
}

int main() {
    Test();
}

Output:

Bits  4, Sieve   0.000001 sec, Loop   0.000000 sec, LoopBoost  30.0x
Bits  8, Sieve   0.000002 sec, Loop   0.000000 sec, LoopBoost  14.9x
Bits 12, Sieve   0.000003 sec, Loop   0.000000 sec, LoopBoost  13.5x
Bits 16, Sieve   0.000005 sec, Loop   0.000000 sec, LoopBoost  11.0x
Bits 20, Sieve   0.000017 sec, Loop   0.000001 sec, LoopBoost  13.6x
Bits 24, Sieve   0.000027 sec, Loop   0.000005 sec, LoopBoost   5.7x
Bits 28, Sieve   0.000081 sec, Loop   0.000016 sec, LoopBoost   5.2x
Bits 32, Sieve   0.000329 sec, Loop   0.000042 sec, LoopBoost   7.8x
Bits 36, Sieve   0.001204 sec, Loop   0.000137 sec, LoopBoost   8.8x
Bits 40, Sieve   0.004266 sec, Loop   0.000506 sec, LoopBoost   8.4x
Bits 44, Sieve   0.015884 sec, Loop   0.001581 sec, LoopBoost  10.0x
Bits 48, Sieve   0.073927 sec, Loop   0.006467 sec, LoopBoost  11.4x
Bits 52, Sieve   0.624178 sec, Loop   0.031799 sec, LoopBoost  19.6x
Bits 56, Sieve   2.458833 sec, Loop   0.048363 sec, LoopBoost  50.8x
Bits 60, Sieve  10.607011 sec, Loop   0.184828 sec, LoopBoost  57.4x
Bits 64, Sieve  94.498041 sec, Loop   0.367877 sec, LoopBoost 256.9x

Answer 2

Your array unsigned long long int primeFlagger[checkNumber+1]; inside sieve func is too long. Use array in global scope, outside of any function or dynamic memory allocation.

Also, you dont need unsigned long long. Its a largest integer data type and you use only one bit of it. Changing type to bool will help you too.

There is other problems:

• unsigned long long int root=(int)(sqrt(checkNumber)); - If number is really large, sqrt(checkNumber) can overflow int;
• unsigned long long int primeFlagger[checkNumber+1]; - type of checkNumber is probably larger than std::size_t - the type for array index and larger than largest memory area that can be allocated. You just cant make array size of unsigned long long.
• checkNumber=numberToCheck; - You dont need this. numberToCheck where already passed into function as parameter checkNumber. Inside sieve checkNummber will be equal to numberToCheck;
• for(unsigned long long int j=2;j<=checkNumber;j++) - this loop should ends then j<=root. This would be enough to mark all non-primal numbers.

If you really need to work with such large numbers, use segmented sieve.