Skip to main content
Update multiplication algorithm as suggested by Keith Randall
Source Link
justhalf
  • 2.3k
  • 1
  • 19
  • 32

C++ (GMP) - (11287,300000,000 / 422,000) = 26680.7809

Shamelessly combine Kummer's Theorem by xnor and GMP by qwr. Still not even closeStill not even close to the Go solution, not sure why.

Edit: Thanks Keith Randall for the reminder that multiplication is faster if the number is similar in size. I implemented multi-level multiplication, similar to memory coalescing concept on memory management. And the Go solutionresult is impressive. What used to take 51s, not sure whynow takes only 0.5s (i.e., 100-fold improvement!!)

#include <gmpxx.h>
#include <iostream>
#include <time.h>
#include <cstdio>
#include <cstring>

const int MAX=12500000;

int primes[819000], factors[819000], count;
bool sieve[MAX];
int max_idx=0;

void run_sieve(){
    sieve[2] = true;
    primes[0] = 2;
    count = 1;
    for(int i=3; i<MAX; i+=2){
        sieve[i] = true;
    }
    for(int i=3; i<4000; i+=2){
        if(!sieve[i]) continue;
        for(int j = i*i; j<MAX; j+=i){
            sieve[j] = false;
        }
    }
    for(int i=3; i<MAX; i+=2){
        if(sieve[i]) primes[count++] = i;
    }
}

mpz_class sum_digits(mpz_class n){
    clock_t t = clock();
    char* str = mpz_get_str(NULL, 10, n.get_mpz_t());
    int result = 0;
    for(int i=0;str[i]>0;i++){
        result+=str[i]-48;
    }
    printf("Done summing in %.3fs\n", ((float)(clock()-t))/CLOCKS_PER_SEC);
    return result;
}

mpz_class nc2_fast(const mpz_class &x){
    clock_t t = clock();
    int prime;
    const unsigned int n = mpz_get_ui(x.get_mpz_t());
    const unsigned int n2 = n/2;
    unsigned int m;
    unsigned int digit;
    unsigned int carry=0;
    unsigned int carries=0;
    mpz_class result = 1;
    mpz_class tmp;
    for(int i=0; i< count; i++){
        prime = primes[i];
        carry=0;
        carries=0;
        if(prime > n) break;
        if(prime > n2){
            mpz_mul_ui(result.get_mpz_t(), result.get_mpz_t(), prime);
            continue;
        }
        m=n2;
        while(m>0){
            digit = m%prime;
            carry = (2*digit + carry >= prime) ? 1 : 0;
            carries += carry;
            m/=prime;
        }
        if(carries>0){
            mpz_ui_pow_ui(tmp.get_mpz_t(), prime, carries);
            mpz_mul(result.get_mpz_t(), result.get_mpz_t(), tmp.get_mpz_t());
        }
    }
    printf("Done calculating binom in %.3fs\n", ((float)(clock()-t))/CLOCKS_PER_SEC);
    return result;
}

int main(){
    const mpz_class n = 11300000;
    clock_t t = clock();
    run_sieve();
    printf("Done sieving in %.3fs\n", ((float)(clock()-t))/CLOCKS_PER_SEC);
    std::cout << sum_digits(nc2_fast(n)) << std::endl;
    return 0;
}

OLD CODE (n=14,000,000)
Done sieving in 0.343s
Done calculating binom in 51.929s
Done summing in 0.901s
14000000: 18954729

real    0m53.194s
user    0m53.116s
sys 0m0.060s

NEW CODE (n=14,000,000)
Done sieving in 0.343s
Done calculating binom in 0.552s
Done summing in 0.902s
14000000: 18954729

real    0m1.804s
user    0m1.776s
sys 0m0.023s

The run for n=287,000,000


Done sieving in 4.211s
Done calculating binom in 17.934s
Done summing in 37.677s
287000000: 388788354

real    0m59.928s
user    0m58.759s
sys 0m1.116s

The code. Compile with -lgmp -lgmpxx -O3

#include <gmpxx.h>
#include <iostream>
#include <time.h>
#include <cstdio>

const int MAX=287000000;
const int PRIME_COUNT=15700000;

int primes[PRIME_COUNT], factors[PRIME_COUNT], count;
bool sieve[MAX];
int max_idx=0;

void run_sieve(){
    sieve[2] = true;
    primes[0] = 2;
    count = 1;
    for(int i=3; i<MAX; i+=2){
        sieve[i] = true;
    }
    for(int i=3; i<17000; i+=2){
        if(!sieve[i]) continue;
        for(int j = i*i; j<MAX; j+=i){
            sieve[j] = false;
        }
    }
    for(int i=3; i<MAX; i+=2){
        if(sieve[i]) primes[count++] = i;
    }
}

mpz_class sum_digits(mpz_class n){
    clock_t t = clock();
    char* str = mpz_get_str(NULL, 10, n.get_mpz_t());
    int result = 0;
    for(int i=0;str[i]>0;i++){
        result+=str[i]-48;
    }
    printf("Done summing in %.3fs\n", ((float)(clock()-t))/CLOCKS_PER_SEC);
    return result;
}

mpz_class nc2_fast(const mpz_class &x){
    clock_t t = clock();
    int prime;
    const unsigned int n = mpz_get_ui(x.get_mpz_t());
    const unsigned int n2 = n/2;
    unsigned int m;
    unsigned int digit;
    unsigned int carry=0;
    unsigned int carries=0;
    mpz_class result = 1;
    mpz_class prime_prods = 1;
    mpz_class tmp;
    mpz_class tmp_prods[32], tmp_prime_prods[32];
    for(int i=0; i<32; i++){
        tmp_prods[i] = (mpz_class)NULL;
        tmp_prime_prods[i] = (mpz_class)NULL;
    }
    for(int i=0; i< count; i++){
        prime = primes[i];
        carry=0;
        carries=0;
        if(prime > n) break;
        if(prime > n2){
            tmp = prime;
            for(int j=0; j<32; j++){
                if(tmp_prime_prods[j] == NULL){
                    tmp_prime_prods[j] = tmp;
                    break;
                } else {
                    mpz_mul(tmp.get_mpz_t(), tmp.get_mpz_t(), tmp_prime_prods[j].get_mpz_t());
                    tmp_prime_prods[j] = (mpz_class)NULL;
                }
            }
            continue;
        }
        m=n2;
        while(m>0){
            digit = m%prime;
            carry = (2*digit + carry >= prime) ? 1 : 0;
            carries += carry;
            m/=prime;
        }
        if(carries>0){
            tmp = 0;
            mpz_ui_pow_ui(tmp.get_mpz_t(), prime, carries);
            for(int j=0; j<32; j++){
                if(tmp_prods[j] == NULL){
                    tmp_prods[j] = tmp;
                    break;
                } else {
                    mpz_mul(tmp.get_mpz_t(), tmp.get_mpz_t(), tmp_prods[j].get_mpz_t());
                    tmp_prods[j] = (mpz_class)NULL;
                }
            }
        }
    }
    result = 1;
    prime_prods = 1;
    for(int j=0; j<32; j++){
        if(tmp_prods[j] != NULL){
            mpz_mul(result.get_mpz_t(), result.get_mpz_t(), tmp_prods[j].get_mpz_t());
        }
        if(tmp_prime_prods[j] != NULL){
            mpz_mul(prime_prods.get_mpz_t(), prime_prods.get_mpz_t(), tmp_prime_prods[j].get_mpz_t());
        }
    }
    mpz_mul(result.get_mpz_t(), result.get_mpz_t(), prime_prods.get_mpz_t());
    printf("Done calculating binom in %.3fs\n", ((float)(clock()-t))/CLOCKS_PER_SEC);
    return result;
}

int main(int argc, char* argv[]){
    const mpz_class n = atoi(argv[1]);
    clock_t t = clock();
    run_sieve();
    printf("Done sieving in %.3fs\n", ((float)(clock()-t))/CLOCKS_PER_SEC);
    std::cout << n << ": " << sum_digits(nc2_fast(n)) << std::endl;
    return 0;
}

C++ (GMP) - (11,300,000 / 422,000) = 26.78

Shamelessly combine Kummer's Theorem by xnor and GMP by qwr. Still not even close to the Go solution, not sure why.

#include <gmpxx.h>
#include <iostream>
#include <time.h>
#include <cstdio>
#include <cstring>

const int MAX=12500000;

int primes[819000], factors[819000], count;
bool sieve[MAX];
int max_idx=0;

void run_sieve(){
    sieve[2] = true;
    primes[0] = 2;
    count = 1;
    for(int i=3; i<MAX; i+=2){
        sieve[i] = true;
    }
    for(int i=3; i<4000; i+=2){
        if(!sieve[i]) continue;
        for(int j = i*i; j<MAX; j+=i){
            sieve[j] = false;
        }
    }
    for(int i=3; i<MAX; i+=2){
        if(sieve[i]) primes[count++] = i;
    }
}

mpz_class sum_digits(mpz_class n){
    clock_t t = clock();
    char* str = mpz_get_str(NULL, 10, n.get_mpz_t());
    int result = 0;
    for(int i=0;str[i]>0;i++){
        result+=str[i]-48;
    }
    printf("Done summing in %.3fs\n", ((float)(clock()-t))/CLOCKS_PER_SEC);
    return result;
}

mpz_class nc2_fast(const mpz_class &x){
    clock_t t = clock();
    int prime;
    const unsigned int n = mpz_get_ui(x.get_mpz_t());
    const unsigned int n2 = n/2;
    unsigned int m;
    unsigned int digit;
    unsigned int carry=0;
    unsigned int carries=0;
    mpz_class result = 1;
    mpz_class tmp;
    for(int i=0; i< count; i++){
        prime = primes[i];
        carry=0;
        carries=0;
        if(prime > n) break;
        if(prime > n2){
            mpz_mul_ui(result.get_mpz_t(), result.get_mpz_t(), prime);
            continue;
        }
        m=n2;
        while(m>0){
            digit = m%prime;
            carry = (2*digit + carry >= prime) ? 1 : 0;
            carries += carry;
            m/=prime;
        }
        if(carries>0){
            mpz_ui_pow_ui(tmp.get_mpz_t(), prime, carries);
            mpz_mul(result.get_mpz_t(), result.get_mpz_t(), tmp.get_mpz_t());
        }
    }
    printf("Done calculating binom in %.3fs\n", ((float)(clock()-t))/CLOCKS_PER_SEC);
    return result;
}

int main(){
    const mpz_class n = 11300000;
    clock_t t = clock();
    run_sieve();
    printf("Done sieving in %.3fs\n", ((float)(clock()-t))/CLOCKS_PER_SEC);
    std::cout << sum_digits(nc2_fast(n)) << std::endl;
    return 0;
}

C++ (GMP) - (287,000,000 / 422,000) = 680.09

Shamelessly combine Kummer's Theorem by xnor and GMP by qwr. Still not even close to the Go solution, not sure why.

Edit: Thanks Keith Randall for the reminder that multiplication is faster if the number is similar in size. I implemented multi-level multiplication, similar to memory coalescing concept on memory management. And the result is impressive. What used to take 51s, now takes only 0.5s (i.e., 100-fold improvement!!)


OLD CODE (n=14,000,000)
Done sieving in 0.343s
Done calculating binom in 51.929s
Done summing in 0.901s
14000000: 18954729

real    0m53.194s
user    0m53.116s
sys 0m0.060s

NEW CODE (n=14,000,000)
Done sieving in 0.343s
Done calculating binom in 0.552s
Done summing in 0.902s
14000000: 18954729

real    0m1.804s
user    0m1.776s
sys 0m0.023s

The run for n=287,000,000


Done sieving in 4.211s
Done calculating binom in 17.934s
Done summing in 37.677s
287000000: 388788354

real    0m59.928s
user    0m58.759s
sys 0m1.116s

The code. Compile with -lgmp -lgmpxx -O3

#include <gmpxx.h>
#include <iostream>
#include <time.h>
#include <cstdio>

const int MAX=287000000;
const int PRIME_COUNT=15700000;

int primes[PRIME_COUNT], factors[PRIME_COUNT], count;
bool sieve[MAX];
int max_idx=0;

void run_sieve(){
    sieve[2] = true;
    primes[0] = 2;
    count = 1;
    for(int i=3; i<MAX; i+=2){
        sieve[i] = true;
    }
    for(int i=3; i<17000; i+=2){
        if(!sieve[i]) continue;
        for(int j = i*i; j<MAX; j+=i){
            sieve[j] = false;
        }
    }
    for(int i=3; i<MAX; i+=2){
        if(sieve[i]) primes[count++] = i;
    }
}

mpz_class sum_digits(mpz_class n){
    clock_t t = clock();
    char* str = mpz_get_str(NULL, 10, n.get_mpz_t());
    int result = 0;
    for(int i=0;str[i]>0;i++){
        result+=str[i]-48;
    }
    printf("Done summing in %.3fs\n", ((float)(clock()-t))/CLOCKS_PER_SEC);
    return result;
}

mpz_class nc2_fast(const mpz_class &x){
    clock_t t = clock();
    int prime;
    const unsigned int n = mpz_get_ui(x.get_mpz_t());
    const unsigned int n2 = n/2;
    unsigned int m;
    unsigned int digit;
    unsigned int carry=0;
    unsigned int carries=0;
    mpz_class result = 1;
    mpz_class prime_prods = 1;
    mpz_class tmp;
    mpz_class tmp_prods[32], tmp_prime_prods[32];
    for(int i=0; i<32; i++){
        tmp_prods[i] = (mpz_class)NULL;
        tmp_prime_prods[i] = (mpz_class)NULL;
    }
    for(int i=0; i< count; i++){
        prime = primes[i];
        carry=0;
        carries=0;
        if(prime > n) break;
        if(prime > n2){
            tmp = prime;
            for(int j=0; j<32; j++){
                if(tmp_prime_prods[j] == NULL){
                    tmp_prime_prods[j] = tmp;
                    break;
                } else {
                    mpz_mul(tmp.get_mpz_t(), tmp.get_mpz_t(), tmp_prime_prods[j].get_mpz_t());
                    tmp_prime_prods[j] = (mpz_class)NULL;
                }
            }
            continue;
        }
        m=n2;
        while(m>0){
            digit = m%prime;
            carry = (2*digit + carry >= prime) ? 1 : 0;
            carries += carry;
            m/=prime;
        }
        if(carries>0){
            tmp = 0;
            mpz_ui_pow_ui(tmp.get_mpz_t(), prime, carries);
            for(int j=0; j<32; j++){
                if(tmp_prods[j] == NULL){
                    tmp_prods[j] = tmp;
                    break;
                } else {
                    mpz_mul(tmp.get_mpz_t(), tmp.get_mpz_t(), tmp_prods[j].get_mpz_t());
                    tmp_prods[j] = (mpz_class)NULL;
                }
            }
        }
    }
    result = 1;
    prime_prods = 1;
    for(int j=0; j<32; j++){
        if(tmp_prods[j] != NULL){
            mpz_mul(result.get_mpz_t(), result.get_mpz_t(), tmp_prods[j].get_mpz_t());
        }
        if(tmp_prime_prods[j] != NULL){
            mpz_mul(prime_prods.get_mpz_t(), prime_prods.get_mpz_t(), tmp_prime_prods[j].get_mpz_t());
        }
    }
    mpz_mul(result.get_mpz_t(), result.get_mpz_t(), prime_prods.get_mpz_t());
    printf("Done calculating binom in %.3fs\n", ((float)(clock()-t))/CLOCKS_PER_SEC);
    return result;
}

int main(int argc, char* argv[]){
    const mpz_class n = atoi(argv[1]);
    clock_t t = clock();
    run_sieve();
    printf("Done sieving in %.3fs\n", ((float)(clock()-t))/CLOCKS_PER_SEC);
    std::cout << n << ": " << sum_digits(nc2_fast(n)) << std::endl;
    return 0;
}
Source Link
justhalf
  • 2.3k
  • 1
  • 19
  • 32

C++ (GMP) - (11,300,000 / 422,000) = 26.78

Shamelessly combine Kummer's Theorem by xnor and GMP by qwr. Still not even close to the Go solution, not sure why.

#include <gmpxx.h>
#include <iostream>
#include <time.h>
#include <cstdio>
#include <cstring>

const int MAX=12500000;

int primes[819000], factors[819000], count;
bool sieve[MAX];
int max_idx=0;

void run_sieve(){
    sieve[2] = true;
    primes[0] = 2;
    count = 1;
    for(int i=3; i<MAX; i+=2){
        sieve[i] = true;
    }
    for(int i=3; i<4000; i+=2){
        if(!sieve[i]) continue;
        for(int j = i*i; j<MAX; j+=i){
            sieve[j] = false;
        }
    }
    for(int i=3; i<MAX; i+=2){
        if(sieve[i]) primes[count++] = i;
    }
}

mpz_class sum_digits(mpz_class n){
    clock_t t = clock();
    char* str = mpz_get_str(NULL, 10, n.get_mpz_t());
    int result = 0;
    for(int i=0;str[i]>0;i++){
        result+=str[i]-48;
    }
    printf("Done summing in %.3fs\n", ((float)(clock()-t))/CLOCKS_PER_SEC);
    return result;
}

mpz_class nc2_fast(const mpz_class &x){
    clock_t t = clock();
    int prime;
    const unsigned int n = mpz_get_ui(x.get_mpz_t());
    const unsigned int n2 = n/2;
    unsigned int m;
    unsigned int digit;
    unsigned int carry=0;
    unsigned int carries=0;
    mpz_class result = 1;
    mpz_class tmp;
    for(int i=0; i< count; i++){
        prime = primes[i];
        carry=0;
        carries=0;
        if(prime > n) break;
        if(prime > n2){
            mpz_mul_ui(result.get_mpz_t(), result.get_mpz_t(), prime);
            continue;
        }
        m=n2;
        while(m>0){
            digit = m%prime;
            carry = (2*digit + carry >= prime) ? 1 : 0;
            carries += carry;
            m/=prime;
        }
        if(carries>0){
            mpz_ui_pow_ui(tmp.get_mpz_t(), prime, carries);
            mpz_mul(result.get_mpz_t(), result.get_mpz_t(), tmp.get_mpz_t());
        }
    }
    printf("Done calculating binom in %.3fs\n", ((float)(clock()-t))/CLOCKS_PER_SEC);
    return result;
}

int main(){
    const mpz_class n = 11300000;
    clock_t t = clock();
    run_sieve();
    printf("Done sieving in %.3fs\n", ((float)(clock()-t))/CLOCKS_PER_SEC);
    std::cout << sum_digits(nc2_fast(n)) << std::endl;
    return 0;
}