# Sum of smallest prime factors

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

• 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 – The Coder Jul 13 '15 at 6:54
• @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. – Nicolás Siplis Jul 13 '15 at 16:18
• Yes, as long as it runs on my own machine and outputs the correct answer in 60 seconds or less, it is acceptable. – Nicolás Siplis Jul 13 '15 at 16:22
• It has 4 threads. – Nicolás Siplis Jul 13 '15 at 16:24
• Does third party Libraries are allowed? It is ok if the program is creating threads? – The Coder Jul 13 '15 at 17:26

## 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)

• 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). – Nicolás Siplis Jul 13 '15 at 20:42
• You also don't need the return in f's definition. Single-expression procs automatically return. – kirbyfan64sos Jul 13 '15 at 20:54
• This isn't the first fastest-code that Nim has won by a noticeable margin. Might be worth investigating. – primo Jul 14 '15 at 0:41
• I'm curious to see how it performs when using GMP, but couldn't implement it correctly despite my efforts. – Nicolás Siplis Jul 14 '15 at 0:43
• Nim is definitely going on my to-learn list! – Sp3000 Jul 14 '15 at 12:36

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; }  • If memory is a concern, one bit per number should be enough. You can use a bitmask to store the flags. – Reto Koradi Jul 13 '15 at 15:48 • @RetoKoradi Unfortunately, that would probably slow the program down enough to put it over the 1 minute mark. – isaacg Jul 13 '15 at 18:03 • What do you need assert.h for? – bot47 Jul 14 '15 at 9:57 • @MaxRied It was left over from an earlie version. – isaacg Jul 14 '15 at 9:57 # Perl, brute force factoring use ntheory ":all"; sub T { my$sum=0;
for (1..$_[0]) {$sum += !($_%2) ? 2 : !($_%3) ? 3 : !($_%5) ? 5 : (factor($_))[0];
}
$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; }  # 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  • 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. – The Coder Jul 16 '15 at 7:16 • @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! – sadakatsu Jul 16 '15 at 12:26 • 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.. – The Coder Jul 16 '15 at 12:30 • You could achieve x2 result if you implement MultiThreading. But it would need to precalculate the whole array before running the Prime calculation. – The Coder Jul 16 '15 at 12:36 • @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. – sadakatsu Jul 16 '15 at 12:42 # 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  • 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) – mawir Jul 13 '15 at 15:04 ## 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


# 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.

# 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.)

• 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)))) – Nicolás Siplis Jul 13 '15 at 17:00
• 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>". – Candles Jul 13 '15 at 18:32
• I tried but still no luck, just in case I tried compiling online with Ideone but that didn't work either. – Nicolás Siplis Jul 13 '15 at 19:08
• Oh, sorry, I forgot the "do" keyword on line 6. It should run now though, give it another shot. – Candles Jul 13 '15 at 19:15
• 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. – Nicolás Siplis Jul 13 '15 at 19:24

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);
}


## 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());
}

}