Largest Prime Factor redux

Question

There have been a number of questions asking for comments on code for the Project Euler problem to compute the largest prime factor of a number. A primary intent of this task is to find the factor as quickly as possible.

I decided to implement the code in C++ using at least some approximation of the best techniques/algorithms of which I was aware¹.

#include <iostream>
#include <vector>
#include <time.h>

typedef unsigned long long integer;

const integer input = 60085147514377;

class prime_table {
    std::vector<integer> primes;
public:
    prime_table(integer num) {
        integer max = ceil(sqrt(num));
        std::vector<bool> table;
        table.resize(max, false);

        for (int i = 2; i < max; i++) {
            if (!table[i]) {
                primes.push_back(i);
                for (integer temp = 2 * i; temp < max; temp += i)
                    table[temp] = true;
            }
        }
    }

    unsigned size() {   return primes.size();   }   

    integer back() { return primes.back(); }

    integer operator[](unsigned index) {    return primes[index];   }

    integer prime_factor(integer input) {
        for (int i = 0; primes[i] * primes[i] <= input; i++)
            if (input % primes[i] == 0)
                return prime_factor(input / primes[i]);
        return input;
    }
};

int main() { 
    prime_table s(input);

    integer current = 1;

    clock_t start = clock();
    for (int i = 0; i < s.size(); i++) {
        if ((input % s[i]) == 0) {
            if (s[i]>current)
                current = s[i];
            integer candidate = input / s[i];
            candidate = s.prime_factor(candidate);
            if (candidate > current) {
                current = candidate;
                if (current > s.back())
                    break;
            }
        }
    }
    clock_t stop = clock();
    std::cout << "\n" << current << "\n";
    std::cout << "Time: " << double(stop - start) / CLOCKS_PER_SEC << " seconds\n";
}

Unfortunately, this has the exact opposite problem most attempts do: rather than taking far too long to run, it runs so fast that much of the time it simply shows the minimum interval (1 millisecond) the clock implementation I'm using can measure. Worse, a substantial part of the time, it shows up as running in 0 milliseconds.

To try to get a better idea of the speed, I chose a somewhat larger number to factor -- specifically I tacked 77 on the end to get of the input suggested by Project Euler to get 60085147514377. The code above gives the largest prime factor as 345937 (which Wolfram Alpha agrees is the largest prime factor of this number) in about 9 or 10 milliseconds.

Specific question: implicit in almost any question like this is one of whether further improvements in speed can be achieved easily. I'm reasonably certain there is room for improvement in the sieving part.

I'm also somewhat uncertain about the number theory part--I think this is correct, and it produces correct results for every test I've tried, but I feel slightly less than 100% certain that there couldn't be some input for which it could produce an incorrect result.

---

^{1. I'm aware, for example, that the sieve of Atkins is at least theoretically faster than the sieve of Eratosthenes, but it apparently requires fairly non-trivial code to actually achieve that theoretical speedup so I haven't attempted to implement it here.}

Answer 1

Slight improvement in the sieve.
I doubt it is significant for max.

prime_table(integer num) {
    integer max = ceil(sqrt(num));
    std::vector<char> table(max, false);   // Making this char rather than bool
                                           // probably makes a tiny difference in speed.
                                           // Don't do an extra re-size do it all up
                                           // front in one allocation.

    // Initialize with the first 2 primes.
    // We don't need to mark multiples of 2 and 3 in the table as true.
    // because of the way we increment in the loop. We basically skip
    // all cells in table {for table[x] where (x%2 ==  0 || x%3 == 0)}
    // so we don't even consider multiples of 2 and 3.
    primes.push_back(2);
    primes.push_back(3);
    int inc = 4;

    // Increment by 2 then 4 (repeat)  (represented by `inc`)
    // This effectively ignores all multiples of 2 and 3.
    for (int i = 5; i < max; i += inc) {
        inc = 6 - inc;
        if (!table[i]) {
            primes.push_back(i);
            for (integer temp = 2 * i; temp < max; temp += i)
                table[temp] = true;
        }
    }
}

You needed the smallest prime here:

integer prime_factor(integer input) {
    for (int i = 0; primes[i] * primes[i] <= input; i++)
        if (input % primes[i] == 0)
            return prime_factor(input / primes[i]);
    return input;
}

This may result in a lot of wasted work (if primes[0] works you did not need the other primes). To only do as much work as you need you may want to try and implement a lazy sieve. That way you only calculate as many primes as you actually need. That being said the sieve is already fast so the extra code and complexity may outweigh the advantages of being lazy.

But to facilitate your tests (to see if a lazy sieve is worthwhile) you could update the code to use std::begin() and std::end() on the primes container. That way it can use your current technique. But if you change the container for primes into something that calculates primes on the fly you don't need to modify your code.

integer prime_factor(integer input)
{
    for (auto prime: primes)   // use std::begin() and std::end in the background (I think)
    {
        if ((prime * prime) > input) {
            break;
        }
        if (input % prime == 0) {
            return prime_factor(input / prime);
        }
    }
    return input;
}

Here is an attempt at a Lazy prime Sieve. Maybe that will help.

class LazyPrime
{
    std::vector<integer>    primes;
    integer                 max;
    std::vector<char>       table;
    integer                 inc;
    integer                 lastPrimeFound;

    struct LazyPrimeIterator
    {
        LazyPrime&   parent;
        std::size_t  current;
        LazyPrimeIterator(LazyPrime& p, int v) : parent(p), current(v) {}
        LazyPrimeIterator& operator++()     {current = parent.next(++current);return *this;}
        integer            operator*()      {return parent.primes[current];}
        bool               operator!=(LazyPrimeIterator const& rhs) const
        {return current != rhs.current;}
    };

    public:
        typedef LazyPrimeIterator           iterator;
        LazyPrime(integer num)
            : max(ceil(sqrt(num)))
            , table(max, false)
            , inc(2)
            , lastPrimeFound(5)
        {
            primes.push_back(2);
            primes.push_back(3);
            primes.push_back(5);
        }

        std::size_t next(std::size_t current)
        {
            if (current < primes.size())
            {    return current;
            }

            // Assume primes retrieved in order.
            assert(current == primes.size());

            for (integer temp = 2 * lastPrimeFound; temp < max; temp += lastPrimeFound){
                table[temp] = true;
            }

            // Increment by 2 then 4 (repeat)  (represented by `inc`)
            // This effectively ignores all multiples of 2 and 3.
            for (integer temp = lastPrimeFound + inc; temp < max; temp += inc) {
                inc = 6 - inc;
                if (!table[temp]) {
                    lastPrimeFound = temp;
                    primes.push_back(lastPrimeFound);
                    return current;
                }
            }
            lastPrimeFound = max - inc;
            return -1; // Same as end.
        }

        iterator begin() { return LazyPrimeIterator(*this, 0);}
        iterator end()   { return LazyPrimeIterator(*this, -1);}
};

Answer 2

The sieve can be optimized by concentrating on eliminating only multiples of odd numbers.

So you can step by 2, and in the inner loop, every multiple up to the prime squared is already eliminated, so it can start there. Since that will be an odd multiple and all the even multiples have already been eliminated, you can step, by 2 times the prime, to eliminate the rest:

Update: By estimating the size of the primes vector that you'll need and resizing again after your done, you eliminate all those costly push backs because now you can use assignments instead:

prime_table(integer num)
{
    integer max = ceil(sqrt(num));
    primes.resize((int)(max / (log(max) - 1.08366)) * 2);
    std::vector<bool> table(max, false);
    primes[0] = 2;
    integer index = 1;
    for (integer i = 3; i < max; i += 2) {
        if (!table[i]) {
            primes[index++] = i;
            integer step = i * 2;
            for (integer temp = i * i; temp < max; temp += step)
                table[temp] = true;
        }
    }
    primes.resize(index);
}

A side note. Since you're using this to find prime factors you're only eliminating to the square root of the input number. This can create problems if you do this for other purposes.