19
\$\begingroup\$

SF(n) is a function which computes the smallest prime factor for a given number n.

We'll call T(N) the sum of every SF(n) with 2 <= n <= N.

T(1) = 0 (the sum is over 0 summands)

T(2) = 2 (2 is the first prime)

T(3) = 5 = 2 + 3

T(4) = 7 = 2 + 3 + 2

T(5) = 12 = 2 + 3 + 2 + 5

...

T(10000) = 5786451

The winner will be the one who manages to compute the largest T(N) in 60 seconds on my own laptop (Toshiba Satellite L845, Intel Core i5, 8GB RAM).


Current top score: Nicolás Siplis - 3.6e13 points - Nim
\$\endgroup\$
  • \$\begingroup\$ Pf(2) = 2, Pf(3) = 3, So, T(3) = 2 + 3 = 5. Am I right? I do programmed to find prime factors, but could you elaborate about the current requirement. Thank you \$\endgroup\$ – The Coder Jul 13 '15 at 6:54
  • 1
    \$\begingroup\$ @ToddLehman I'm running each code in my own laptop (Sony Vaio SVF14A16CLB), so if it takes less than 60 seconds I'll increase the number and decrease it when it takes longer. \$\endgroup\$ – Nicolás Siplis Jul 13 '15 at 16:18
  • 1
    \$\begingroup\$ Yes, as long as it runs on my own machine and outputs the correct answer in 60 seconds or less, it is acceptable. \$\endgroup\$ – Nicolás Siplis Jul 13 '15 at 16:22
  • 1
    \$\begingroup\$ It has 4 threads. \$\endgroup\$ – Nicolás Siplis Jul 13 '15 at 16:24
  • 1
    \$\begingroup\$ Does third party Libraries are allowed? It is ok if the program is creating threads? \$\endgroup\$ – The Coder Jul 13 '15 at 17:26

12 Answers 12

12
\$\begingroup\$

Nim, 3.6e13

Simply sieving is not the best answer when trying to calculate the highest N possible since the memory requirements become too high. Here's a different approach (started with Nim a couple days ago and fell in love with the speed and syntax, any suggestions to make it faster or more readable are welcome!).

import math
import sequtils
import nimlongint # https://bitbucket.org/behrends/nimlongint/

proc s(n : int) : int128 =
    var x = toInt128(n)
    (x * x + x) div 2 - 1

proc sum_pfactor(N : int) : int128 =    
    var
        root = int(sqrt(float(N)))
        u = newSeqWith(root+1,false)
        cntA,cntB,sumA,sumB = newSeq[int128](root+1)
        pcnt,psum,ret : int128
        interval,finish,d,q,t : int

    for i in 0..root:
        cntA[i] = i-1
        sumA[i] = s(i)

    for i in 1..root:
        cntB[i] = N div i - 1
        sumB[i] = s(N div i)

    for p in 2..root:
        if cntA[p] == cntA[p-1]:
            continue

        pcnt = cntA[p - 1]
        psum = sumA[p - 1]
        q = p * p
        ret = ret + p * (cntB[p] - pcnt)
        cntB[1] = cntB[1] - cntB[p] + pcnt
        sumB[1] = sumB[1] - (sumB[p] - psum) * p
        interval = (p and 1) + 1
        finish = min(root,N div q)

        for i in countup(p+interval,finish,interval):

            if u[i]:
                continue

            d = i * p

            if d <= root:
                cntB[i] = cntB[i] - cntB[d] + pcnt
                sumB[i] = sumB[i] - (sumB[d] - psum) * p
            else:
                t = N div d
                cntB[i] = cntB[i] - cntA[t] + pcnt
                sumB[i] = sumB[i] - (sumA[t] - psum) * p

        if q <= root:
            for i in countup(q,finish-1,p*interval):
                u[i] = true

        for i in countdown(root,q-1):
            t = i div p
            cntA[i] = cntA[i] - cntA[t] + pcnt
            sumA[i] = sumA[i] - (sumA[t] - psum) * p

    sumB[1] + ret

var time = cpuTime()
echo(sum_pfactor(int(3.6e13))," - ",cpuTime() - time)
\$\endgroup\$
  • \$\begingroup\$ I tried implementing the GMP wrapper for Nim into my code but just couldn't make it work (never used GMP before so that certainly didn't help). \$\endgroup\$ – Nicolás Siplis Jul 13 '15 at 20:42
  • \$\begingroup\$ You also don't need the return in f's definition. Single-expression procs automatically return. \$\endgroup\$ – kirbyfan64sos Jul 13 '15 at 20:54
  • 3
    \$\begingroup\$ This isn't the first fastest-code that Nim has won by a noticeable margin. Might be worth investigating. \$\endgroup\$ – primo Jul 14 '15 at 0:41
  • \$\begingroup\$ I'm curious to see how it performs when using GMP, but couldn't implement it correctly despite my efforts. \$\endgroup\$ – Nicolás Siplis Jul 14 '15 at 0:43
  • \$\begingroup\$ Nim is definitely going on my to-learn list! \$\endgroup\$ – Sp3000 Jul 14 '15 at 12:36
5
\$\begingroup\$

C, Prime Sieve: 5e9

Results:

$ time ./sieve 
Finding sum of lowest divisors of n = 2..5000000000
572843021990627911

real    0m57.144s
user    0m56.732s
sys 0m0.456s 

Program:

While it's a rather staightforward program, it took me a while to figure out how to get the memory management right - I only have enough ram for 1 byte per number in the range, so I had to be careful. It's a standard sieve of Erasthones.

#include<stdio.h>
#include<stdlib.h>
#include<math.h>
#include<assert.h>

#define LIMIT ((unsigned long long)5e9 +1)
#define ROOT_LIMIT floor(sqrt(LIMIT))

int main()
{
    printf("Finding sum of lowest divisors of n = 2..%llu\n", LIMIT - 1);
    char * found_divisor;
    found_divisor = malloc(LIMIT * sizeof(char));
    if (found_divisor == NULL) {
        printf("Error on malloc");
        return -1;
    }
    unsigned long long i;
    unsigned long long trial_div;
    unsigned long long multiple;
    unsigned long long sum = 0;

    for (i = 0; i < LIMIT; ++i) {
        found_divisor[i] = 0;
    }

    for (trial_div = 2; trial_div <= ROOT_LIMIT; ++trial_div) {
        if (found_divisor[trial_div] == 0) {
            for (multiple = trial_div * trial_div; multiple < LIMIT; multiple += trial_div) {
                if (found_divisor[multiple] == 0) {
                    found_divisor[multiple] = 1;
                    sum += trial_div;
                }
            }
        }
    }

    for (i = 2; i < LIMIT; ++i) {
        if (found_divisor[i] == 0) {
            sum += i;
        }
    }

    free(found_divisor);
    printf("%lld\n", sum);
    return 0;
}
\$\endgroup\$
  • 1
    \$\begingroup\$ If memory is a concern, one bit per number should be enough. You can use a bitmask to store the flags. \$\endgroup\$ – Reto Koradi Jul 13 '15 at 15:48
  • \$\begingroup\$ @RetoKoradi Unfortunately, that would probably slow the program down enough to put it over the 1 minute mark. \$\endgroup\$ – isaacg Jul 13 '15 at 18:03
  • \$\begingroup\$ What do you need assert.h for? \$\endgroup\$ – Max Ried Jul 14 '15 at 9:57
  • \$\begingroup\$ @MaxRied It was left over from an earlie version. \$\endgroup\$ – isaacg Jul 14 '15 at 9:57
3
\$\begingroup\$

Perl, brute force factoring

use ntheory ":all";
sub T {
  my $sum=0;
  for (1..$_[0]) {
    $sum += !($_%2) ? 2 : !($_%3) ? 3 : !($_%5) ? 5 : (factor($_))[0];
  }
  $sum
}
T(90_000_000);

I can get to about 9e7 in 25 seconds on my Linux machine. It could be faster by digging into the C code, as it is saying after a check for 2/3/5, completely factor the number.

There are much more clever ways to do this using sieving. I thought a simple brute force way would be a start. This is basically Project Euler problem 521, by the way.

\$\endgroup\$
  • \$\begingroup\$ If it's at all useful to know, in Python with a sieve I can only manage T(47000). I'm going to try something similar to what you're doing to see if it's faster. \$\endgroup\$ – Kade Jul 13 '15 at 1:02
  • \$\begingroup\$ Looks like not using a sieve is faster.. I was able to calculate T(493900) with a similar method to yours. \$\endgroup\$ – Kade Jul 13 '15 at 1:15
  • \$\begingroup\$ Never used Perl before but I managed to verify your answer, I'll add you to the list! \$\endgroup\$ – Nicolás Siplis Jul 13 '15 at 1:31
  • \$\begingroup\$ To be fair, this uses my module that does the factoring in C (you can force it to use pure Perl for everything, but of course it isn't as fast). \$\endgroup\$ – DanaJ Jul 13 '15 at 1:53
  • \$\begingroup\$ The answer can be computed using any combination of languages, so that's OK. \$\endgroup\$ – Nicolás Siplis Jul 13 '15 at 2:04
3
\$\begingroup\$

Go, 21e9

Does a sieve to find the minimum factor of each number <= N. Spawns goroutines to count sections of the number space.

Run with "go run prime.go -P 4 -N 21000000000".

package main

import (
    "flag"
    "fmt"
    "runtime"
)

const S = 1 << 16

func main() {
    var N, P int
    flag.IntVar(&N, "N", 10000, "N")
    flag.IntVar(&P, "P", 4, "number of goroutines to use")
    flag.Parse()
    fmt.Printf("N = %d\n", N)
    fmt.Printf("P = %d\n", P)
    runtime.GOMAXPROCS(P)

    // Spawn goroutines to check sections of the number range.
    c := make(chan uint64, P)
    for i := 0; i < P; i++ {
        a := 2 + (N-1)*i/P
        b := 2 + (N-1)*(i+1)/P
        go process(a, b, c)
    }
    var sum uint64
    for i := 0; i < P; i++ {
        sum += <-c
    }
    fmt.Printf("T(%d) = %d\n", N, sum)
}

func process(a, b int, res chan uint64) {
    // Find primes up to sqrt(b).  Compute starting offsets.
    var primes []int
    var offsets []int
    for p := 2; p*p < b; p++ {
        if !prime(p) {
            continue
        }
        primes = append(primes, p)
        off := a % p
        if off != 0 {
            off = p - off
        }
        offsets = append(offsets, off)
    }

    // Allocate sieve array.
    composite := make([]bool, S)

    // Check factors of numbers up to b, a block of S at a time.
    var sum uint64
    for ; a < b; a += S {
        runtime.Gosched()
        // Check divisibility of [a,a+S) by our set of primes.
        for i, p := range primes {
            off := offsets[i]
            for ; off < S; off += p {
                if composite[off] {
                    continue // Divisible by a smaller prime.
                }
                composite[off] = true
                if a+off < b {
                    sum += uint64(p)
                }
            }
            // Remember offset for next block.
            offsets[i] = off - S
        }
        // Any remaining numbers are prime.
        for i := 0; i < S; i++ {
            if composite[i] {
                composite[i] = false // Reset for next block.
                continue
            }
            if a+i < b {
                sum += uint64(a + i)
            }
        }
    }
    res <- sum
}

func prime(n int) bool {
    for i := 2; i*i <= n; i++ {
        if n%i == 0 {
            return false
        }
    }
    return true
}

Note that the answer for N=21e9 is between 2^63 and 2^64, so I had to use unsigned 64-bit ints to count correctly...

\$\endgroup\$
  • \$\begingroup\$ I had to modify it to run on my machine (decreased N to 1e9) but the runtime itself is pretty fast, good work! \$\endgroup\$ – Nicolás Siplis Jul 13 '15 at 16:50
  • \$\begingroup\$ @NicolásSiplis: memory usage has been fixed. \$\endgroup\$ – Keith Randall Jul 13 '15 at 19:20
  • \$\begingroup\$ Runtime was 80 seconds but 1.6e10 was calculated in almost exactly 60! \$\endgroup\$ – Nicolás Siplis Jul 13 '15 at 19:41
2
\$\begingroup\$

C++, 1<<34 ~ 1.7e10

Intel(R) Core(TM) i5-2410M CPU @ 2.30GHz

$ g++ -O2 test3.cpp 
$ time ./a.out 
6400765038917999291

real    0m49.640s
user    0m49.610s
sys 0m0.000s
#include <iostream>
#include <vector>

using namespace std;

const long long root = 1 << 17; // must be a power of two to simplify modulo operation
const long long rootd2 = root >> 1;
const long long rootd2m1 = rootd2 - 1;
const long long mult = root; // must be less than or equal to root
const long long n = root * mult; // unused constant (function argument)

int main() {
  vector < int > sieve(rootd2, 0);
  vector < int > primes;
  vector < long long > nexts;
  primes.reserve(root);
  nexts.reserve(root);
  // initialize sum with result for even numbers
  long long sum = n / 2 * 2;
  // sieve of Erathosthenes for numbers less than root
  // all even numbers are skipped
  for(long long i = 1; i < rootd2; ++i){
    if(sieve[i]){
      sieve[i] = 0;
      continue;
    }
    const long long val = i * 2 + 1;
    primes.push_back(val);
    sum += val;
    long long j;
    for(j = (val + 1) * i; j < rootd2; j += val){
      sum += val * (1 - sieve[j]); // conditionals replaced by multiplication
      sieve[j] = 1;
    }
    nexts.push_back(j);
  }
  int k = primes.size();
  long long last = rootd2;
  // segmented sieve of Erathosthenes
  // all even numbers are skipped
  for(int segment = 2; segment <= mult; ++segment){
    last += rootd2;
    for(int i = 0; i < k; ++i){
      long long next = nexts[i];
      long long prime = primes[i];
      if(next < last){
        long long ptr = next & rootd2m1; // modulo replaced by bitmasking
        while(ptr < rootd2){
          sum += prime * (1 - sieve[ptr]); // conditionals replaced by multiplication
          sieve[ptr] = 1;
          ptr += prime;
        }
        nexts[i] = (next & ~rootd2m1) + ptr;
      }
    }
    for(int i = 0; i < rootd2; ++i){
      sum += ((segment - 1) * root + i * 2 + 1) * (1 - sieve[i]);
      sieve[i] = 0;
    }
  }
  cout << sum << endl;
}
\$\endgroup\$
2
\$\begingroup\$

Java 8: 1.8e8 2.4e8

This entry does not compare to several of the other ones already up, but I wanted to post my answer since I had fun working on this.

The main optimizations of my approach are as follow:

  • Every even number has a smallest factor of 2, so these can be added for free after every odd number is processed. Basically, if you have done the work to calculate T(N) when N % 2 == 1, you know that T(N + 1) == T(N) + 2. This allows me to start my counting at three and to increment by iteration by twos.
  • I store my prime numbers in an array as opposed to a Collection type. This more than doubled the N I can reach.
  • I use the prime numbers to factor a number as opposed to performing the Sieve of Eratosthenes. This means that my memory storage is restricted almost completely to my primes array.
  • I store the square root of the number for which I am trying to find the smallest factor. I tried @user1354678's approach of squaring a prime factor each time, but this cost me about than 1e7 from my score.

That's about all there is to it. My code iterates from 3 on by twos until it detects that it has hit or exceeded the time limit, at which point it spits out the answer.

package sum_of_smallest_factors;

public final class SumOfSmallestFactors {
    private static class Result {
        private final int number;
        int getNumber() {
            return number;
        }

        private final long sum;
        long getSum() {
            return sum;
        }


        Result(int number, long sum) {
            this.number = number;
            this.sum = sum;
        }
    }


    private static final long TIME_LIMIT = 60_000_000_000L; // 60 seconds x 1e9 nanoseconds / second


    public static void main(String[] args) {
        SumOfSmallestFactors main = new SumOfSmallestFactors();
        Result result = main.run();
        int number = result.getNumber();
        long sum = result.getSum();
        System.out.format("T(%,d) = %,d\n", number, sum);
    }


    private int[] primes = new int[16_777_216];
    private int primeCount = 0;
    private long startTime;


    private SumOfSmallestFactors() {}

    private Result run() {
        startClock();
        int number;
        long sumOfSmallestFactors = 2;
        for (number = 3; mayContinue(); number += 2) {
            int smallestFactor = getSmallestFactor(number);
            if (smallestFactor == number) {
                addPrime(number);
            }
            sumOfSmallestFactors += smallestFactor + 2;
        }
        --number;

        Result result = new Result(number, sumOfSmallestFactors);
        return result;
    }

    private void startClock() {
        startTime = System.nanoTime();
    }

    private boolean mayContinue() {
        long currentTime = System.nanoTime();
        long elapsedTime = currentTime - startTime;
        boolean result = (elapsedTime < TIME_LIMIT);
        return result;
    }

    private int getSmallestFactor(int number) {
        int smallestFactor = number;
        int squareRoot = (int) Math.ceil(Math.sqrt(number));

        int index;
        int prime = 3;
        for (index = 0; index < primeCount; ++index) {
            prime = primes[index];

            if (prime > squareRoot) {
                break;
            }

            int remainder = number % prime;
            if (remainder == 0) {
                smallestFactor = prime;
                break;
            }
        }

        return smallestFactor;
    }

    private void addPrime(int prime) {
        primes[primeCount] = prime;
        ++primeCount;
    }
}

Running on a different system (Windows 8.1, Intel core i7 @ 2.5 GHz, 8 GB RAM) with the latest version of Java 8 has markedly better results with no code changes:

T(240,308,208) = 1,537,216,753,010,879
\$\endgroup\$
  • \$\begingroup\$ If you could replace the mayContinue() in for loop condition with just a simple condition, you could achieve higher result. And I like your way of precalculating even sum, then incrementing by two. \$\endgroup\$ – The Coder Jul 16 '15 at 7:16
  • \$\begingroup\$ @user1354678, Thanks for the recommendation. Strangely, it did not work. I tried variations of this code on a different computer and found that the posted version is the fastest one. Removing the clock calls from the code and using a simple threshhold number cost me a little over a second. I even tried switching my startTime to an endTime to eliminate the ~2e7 subtractions, but that cost me 3e7 from my score! \$\endgroup\$ – sadakatsu Jul 16 '15 at 12:26
  • \$\begingroup\$ Did you tried it with System.nanoTime() - startTime < TIME_LIMIT, coz it fastens your code a little for me. It's not blazingly fast, considering the fact, this condition is checked millions of times, it'll be a little fast. One thing I learned from your code is, don't put for inside a for.. After moving for to another method in my code, my code speed gets increased by 40%, Thanks.. One thing I'm still figuring out is, Did arrays are much efficient than ArrayList when considering the fact that it's fetched millions of times.. \$\endgroup\$ – The Coder Jul 16 '15 at 12:30
  • \$\begingroup\$ You could achieve x2 result if you implement MultiThreading. But it would need to precalculate the whole array before running the Prime calculation. \$\endgroup\$ – The Coder Jul 16 '15 at 12:36
  • \$\begingroup\$ @user1354678, moving the check from the mayContinue() method into the for loop costs me 8e6 from my score. This may be an issue of local optimizations. I experimented with several data types for storing the primes when I was developing this solution. I was only able to reach 8.8e7 with ArrayList, but I hit 1.8e8 (now 2.4e8) using an array. There may be some performance boosts involved with the lookup, but there are definite boosts for the memory allocation. I have thought about multi-threading the algorithm, but I ran into problems. \$\endgroup\$ – sadakatsu Jul 16 '15 at 12:42
1
\$\begingroup\$

R, 2.5e7

Simple minded sieve of Eratosthenes, vectorised as much as possible. R's not really designed for this sort of problem and I'm pretty sure it can be made faster.

MAX <- 2.5e7
Tot <- 0
vec <- 2:MAX 
while(TRUE) {
    if (vec[1]*vec[1] > vec[length(vec)]) {
        Tot <- Tot + sum(as.numeric(vec))
        break
    }

    fact <- which(vec %% vec[1] == 0)
    Tot <- Tot + vec[1]*length(vec[fact])
    vec <- vec[-fact]
}
Tot
\$\endgroup\$
  • \$\begingroup\$ Fair point about T. 2:MAX is a vector of integers so for large values of MAX, sum(vec) leads to an integer overflow and returns NA. sum(as.numeric(vec)) is summing a vector of doubles which doesn't overflow (although it might not give the right answer) \$\endgroup\$ – mawir Jul 13 '15 at 15:04
1
\$\begingroup\$

Python, ~7e8

Using an incremental Sieve of Erathostenes. Some care does need to be taken that a marked value is marked with its lowest divisor, but the implementation is otherwise fairly straight forward.

Timing was taken with PyPy 2.6.0, input is accepted as a command line argument.

from sys import argv
from math import sqrt

n = int(argv[1])
sieve = {}
imax = int(sqrt(n))

t = n & -2
i = 3
while i <= n:
  divs = sieve.pop(i, [])
  if divs:
    t += divs[-1]
    for v in divs:
      sieve.setdefault(i+v+v, []).append(v)
  else:
    t += i
    if i <= imax: sieve[i*i] = [i]
  i += 2

print t

Sample Usage

$ pypy sum-lpf.py 10000
5786451

$ pypy sum-lpf.py 100000000
279218813374515
\$\endgroup\$
0
\$\begingroup\$

Julia, 5e7

Surely Julia can do better but this is what I have for now. This does 5e7 in about 60 seconds on JuliaBox but I can't test locally yet. Hopefully by then I'll have thought of a more clever approach.

const PRIMES = primes(2^16)

function lpf(n::Int64)
    isprime(n) && return n
    for p in PRIMES
        n % p == 0 && return p
    end
end

function T(N::Int64)
    local x::Int64
    x = @parallel (+) for i = 2:N
        lpf(i)
    end
    x
end

Here we're creating a function lpf that iterates through sequential primes and checks the input for divisibility by each. The function returns the first divisor encountered, thereby obtaining the least prime factor.

The main function computes lpf on the integers from 2 to the input in parallel and reduces the result by summing.

\$\endgroup\$
0
\$\begingroup\$

Common Lisp, 1e7

(defvar input 10000000)
(defvar numbers (loop for i from 2 to input collect i))
(defvar counter)
(defvar primes)

(setf primes (loop for i from 2 to (floor (sqrt input))
    when (loop for j in primes
        do (if (eq (mod i j) 0) (return nil))
        finally (return t))
    collect i into primes
    finally (return primes)))

(format t "~A~%"    
    (loop for i in primes
        do (setf counter 0)
        summing (progn (setf numbers (remove-if #'(lambda (x) (if (eq (mod x i) 0) (progn (incf counter) t))) numbers))
                (* i counter)) into total
        finally (return (+ total (reduce #'+ numbers)))))

I've opted to first generate a list of prime numbers from 2 to (sqrt input), then test every value with the primes, whereas previously I would test against every number up to (sqrt input), which would be pointless (e.g., if a number is divisible by 4, it's also divisible by 2, so it's already accounted for.)

Thank God for side-effects while I'm at it. The remove-if both lowers the size of the list and counts how many elements were removed, so I just have to multiply that by whatever value the loop is on and add that into the running total.

(Fun fact: delete is the destructive equivalent of remove, but for whatever reason, delete is all sorts of slower than remove in this case.)

\$\endgroup\$
  • \$\begingroup\$ Never used Lisp before, I'm getting a compiler error when trying to run your code: (defvar total 0) (defvar counter 0) (defvar input 10000) (defvar numbers (loop for i from 2 to input collect i)) (loop for i from 2 to (floor (sqrt input)) (setf counter 0) summing (prog2 (nsubstitute-if 0 #'(lambda (x) (if (eq (mod x i) 0) (progn (incf counter) t))) numbers) (* i counter) (setf numbers (remove 0 numbers))) into total finally (return (+ total (reduce #'+ numbers)))) \$\endgroup\$ – Nicolás Siplis Jul 13 '15 at 17:00
  • \$\begingroup\$ I'm using SBCL 1.0.38, but when I get home I'll update to the latest version and see how it goes. If you save it into a file, you can run it with "sbcl --script <filename>". \$\endgroup\$ – Candles Jul 13 '15 at 18:32
  • \$\begingroup\$ I tried but still no luck, just in case I tried compiling online with Ideone but that didn't work either. \$\endgroup\$ – Nicolás Siplis Jul 13 '15 at 19:08
  • \$\begingroup\$ Oh, sorry, I forgot the "do" keyword on line 6. It should run now though, give it another shot. \$\endgroup\$ – Candles Jul 13 '15 at 19:15
  • \$\begingroup\$ Great, it computes 6e6 in 60 seconds on my machine! By the way, if I decide to enter my own code, do you know if I should submit it as an answer? I'm not sure if that would allow new submissions. \$\endgroup\$ – Nicolás Siplis Jul 13 '15 at 19:24
0
\$\begingroup\$

Rust 1.5e9

A very naive Eratosthene sieve, but I felt Rust did not receive any love here !

// Expected (approximate) number of primes
fn hint(n:usize) -> usize {
    if n < 2 { 
        1
    } else {
        n / ((n as f64).ln() as usize) + 1
    }
}

fn main() {
    let n:usize = match std::env::args().nth(1) {
        Some(s) => s.parse().ok().expect("Please enter a number !"),
        None => 10000,
    };
    let mut primes = Vec::with_capacity(hint(n));
    let mut sqrt = 2;
    let s = (2..).map(|n:u32| -> u32 {
        if (sqrt * sqrt) < n {
            sqrt += 1;
        }
        let (div, unseen) = match primes.iter().take_while(|&p| *p <= sqrt).filter(|&p| n % p == 0).next() {
            Some(p) => (*p, false),
            None => (n, true),
        };
        if unseen {
            primes.push(div);
        }
        div
    }).take(n-1).fold(0, |acc, p| acc + p);
    println!("{}", s);
}
\$\endgroup\$
0
\$\begingroup\$

Java 2.14e9

Pure Sieve of Eratosthenes with advantage of BitSet

I have calculated the Smallest Prime factor sum upto Integer.MAX_VALUE - 1 just in 33.89 s. But I can't proceed any larger because any further will lead to Integer Overflow on the Size of Bitset. So I'm working on to create another Bitset for the next set of Ranges. Until then, this is the fastest I'm able to generate.


T(214,74,83,646) = 109931450137817286 in 33.89 s
aka
T(2,147,483,646) = 109931450137817286 in 33.89 s

import java.util.BitSet;

public class SmallPrimeFactorSum {

    static int    limit     = Integer.MAX_VALUE - 1;

    // BitSet is highly efficient against boolean[] when Billion numbers were involved
    // BitSet uses only 1 bit for each number
    // boolean[] uses 8 bits aka 1 byte for each number which will produce memory issues for large numbers
    static BitSet primes    = new BitSet(limit + 1);
    static int    limitSqrt = (int) Math.ceil(Math.sqrt(limit));

    static long   start     = System.nanoTime();

    static long   sum       = 0;

    public static void main(String[] args) {
        genPrimes();
    }

    // Generate Primes by Sieve of Eratosthenes
    // Sieve of Eratosthenes is much efficient than Sieve of Atkins as
    // Sieve of Atkins involes Division, Modulus, Multiplication, Subtraction, Addition but
    // Sieve of Eratosthenes involves only addition
    static void genPrimes() {

        // Inverse the Bit values
        primes.flip(0, limit + 1);

        // Now all Values in primes will now be true,
        // True  represents     prime number 
        // False represents not prime number

        // Set 0 and 1 as not Prime number
        primes.clear(0, 2);

        // Set all multiples of each Prime as not Prime;
        for ( int prime = 2; prime > 0 && prime <= limit && prime > 0; prime = primes.nextSetBit(prime + 1) ) {
            // Add Current Prime as its Prime factor
            sum += prime;
            // Skip marking if Current Prime > SQRT(limit)
            if ( prime > limitSqrt ) {
                continue;
            }
            // Mark Every multiple of current Prime as not Prime
            for ( int multiple = prime + prime; multiple <= limit && multiple > 0; multiple += prime ) {
                // Mark as not Prime only if it's true already
                if ( primes.get(multiple) ) {
                    // Add Current Prime as multiple's Prime factor
                    sum += prime;
                    primes.clear(multiple);
                }
            }
        }

        System.out.printf("T(%d) = %d in %.2f s", limit, sum, (System.nanoTime() - start) / 1000000000.0);
        //System.out.printf("Total Primes upto %d : %d\n", limit, primes.cardinality());
    }

}
\$\endgroup\$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.