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