Arithmetic Derivative Problem (Project Euler)

arithmetic-derivative.cpp

#include "primes.hpp"
#include <vector>
#include <iostream>
// #include <sstream>
#include <cstdio>
// #include <thread>
// #include <future>
// #include <queue>
#include <algorithm>

using namespace std;

struct Factor {
	int64_t c, n;

	Factor (int64_t c, int64_t n) : c(c), n(n) {}
};


// Equivalent to the python factorize() function
std::vector<Factor> factorize (int64_t x, int64_t & coeff_product, const std::vector<int64_t> & prime_list) {
	x = abs(x);
	coeff_product = 1;
	if (x <= 1)
		return { {x, 1} };

	std::vector<Factor> factors;

	for (auto p : prime_list) {
		if (x % p == 0) {
			x /= p;
			coeff_product *= p;
			int64_t n = 1;
			while (x % p == 0) {
				x /= p;
				++n;
			}
			factors.emplace_back(p, n);
			if (x == 1)
				return factors;
		}
	}

	cerr << "-- Warning: " << x << " is not factorizable (prime list needs to be expanded) --\n";

	factors.emplace_back(x, 1);
	return factors;
}

// equiv to ad2()
int64_t arith_derivative (int64_t x, const std::vector<Factor> & factors) {
	
	int64_t sum = 0;
	for (auto factor : factors) {
		sum += factor.n * x / factor.c;
	}
	return sum;
}

// equiv to ad_gcd3()
int64_t ad_gcd (int64_t x, const std::vector<Factor> & factors, int64_t coeff_product) {

	int64_t term_sum = 0;
	for (auto factor : factors) {
		term_sum += factor.n * coeff_product / factor.c;
	}
	int64_t uncommon_terms = 1;
	for (auto factor : factors) {
		if (term_sum % factor.c != 0) {
			uncommon_terms *= factor.c;
		}
	}
	return x / uncommon_terms;
}

// Unused (nice thing about using signed integers is that it's trivial to detect overflow)
bool overflow (int64_t x) { return x < 0; }

// Overload (actually, *this* is equiv to ad_gcd3(); the prev. fcn deals with already factorized terms)
int64_t ad_gcd (int64_t x, int64_t & derivative, const std::vector<int64_t> & prime_list) {
	int64_t coeff_product;
	auto factors = factorize(x, coeff_product, prime_list);
	derivative = arith_derivative(x, factors);
	return ad_gcd(x, factors, coeff_product);
}

// Calculates the sum of ad_gcd over a given integer range (eg. [ 1 .. 1e5 ]), and prints the results
int64_t print_ad_gcd_sum (int64_t start, int64_t end, const std::vector<int64_t> & prime_list, bool print_all, std::ostream & os) {

	int64_t sum = 0;
	while (start < end) {
		int64_t d, r = ad_gcd(start++, d, prime_list);
		sum += r;
		if (print_all)
			// printf("%lld\t%lld\t%lld\t%lld\n", start, d, r, sum);
			os << start << '\t' << d << '\t' << r << '\t' << sum << '\n';
	}
	os << "sum = " << sum << endl;
	// cout.flush();
	return sum;
}

//
// Multithreading code removed (was massively broken)
//


int main (int argc, const char **argv) {

	// if (argc != 2) {
	// 	cout << "usage: " << argv[0] << " <count>\n";
	// 	return -1;
	// }

	// int64_t count = atoi(argv[1]);
	// if (count < 1) {
	// 	cout << "err: count must be >= 1\n";
	// 	return -1;
	// }

	int64_t count = 5e3;
	bool verbose = true;

	// run_parallel_ad_gcd_sum(2, 1, count, generate_primes(count));
	// parallel_compute_ad_gcd_sum(4, 1, count, generate_primes(count), verbose);
	print_ad_gcd_sum(1, count, generate_primes(count), verbose, cout);
}

arithmetic_derivative.py

import itertools

def rwh_primes(n):
	# http://stackoverflow.com/questions/2068372/fastest-way-to-list-all-primes-below-n-in-python/3035188#3035188
	""" Returns  a list of primes < n """
	sieve = [True] * n
	for i in xrange(3,int(n**0.5)+1,2):
		if sieve[i]:
			sieve[i*i::2*i]=[False]*((n-i*i-1)/(2*i)+1)
	return [2] + [i for i in xrange(3,n,2) if sieve[i]]
print("Calculating primes...")
primes = rwh_primes(5000)
print("done.")

def factorize (x):
	''' Computes the prime factors of x (in terms of x = product([ c ** n for (c, n) in factors) ]))
		using trial division. '''
	x = abs(x)
	if x <= 1:
		yield (x, 1)
		return
	for p in primes:
		if x % p == 0:
			n = 1
			x /= p
			while x % p == 0:
				n += 1
				x /= p
			yield (p, n)
			if x == 1:
				return
	# If we reach here, then we have exhausted the list of primes, and x has at least one *really* large
	# prime factor. Not quite sure how to handle this, so currently just prints a warning to the user
	# (problem will go away if the prime list is sufficiently large), and return the remainder (x) as
	# a factor. Since x may or may not be prime, this will probably cause errors in the ad_gcd routine,
	# but possibly won't. Either way, just fix it by resizing the list >.<
	print("\n-- Warning: %d is not factorizable (prime list needs to be expanded) --\n"%x)
	yield (x, 1)

def ad_gcd3 (x):
	''' Computes the gcd of x, x' (where x' is the arithmetic derivative of x). Third (current) impl '''
	terms = []
	coeff_product = 1
	for (c, k) in factorize(x):
		coeff_product *= c
		terms += [(c, k)]

	# terms only in x'
	ad_term_sum = 0	# term only in x' (sum of k[i] * product(c[j] for all j != i) for all i)
	for (c, k) in terms:
		ad_term_sum += k * coeff_product / c
	uncommon_terms = 1
	for (c, _) in terms:
		if ad_term_sum % c != 0: # not a common factor of x or x'
			uncommon_terms *= c
	return x / uncommon_terms

def ad_gcd_ (terms):
	''' Used by ad_gcd (original impl) '''
	not_identity = lambda (x, n): n != 0
	minus_1 = lambda (x, n): (x, n-1)
	common_terms1 = filter(not_identity, map(minus_1, terms))

	term_product = 1
	for (x, n) in terms:
		term_product *= x

	term_sum = 0
	for i in range(len(terms)):
		_, p = terms[i]
		for j in range(len(terms)):
			if j != i:
				x, _ = terms[j]
				p *= x
		term_sum += p

	#common_terms2 = factorize(gcd(term_sum, term_product))
	common_terms2 = intersect(factorize(term_sum), factorize(term_product))
	return combine(common_terms1, common_terms2)


def ad_gcd (x):
	''' Computes the gcd of x, x' (original implementation). '''
	terms = ad_gcd_(factorize(x))
	terms = map(lambda (x, n): x ** n, terms)
	return reduce(lambda a, b: a * b, terms)

def ad (x):
	''' Computes the arithmetic deriviative of x (first, inefficient impl) '''
	common_term_product = 1
	terms_, coeff_product = [], 1
	for (c, n) in factorize(x):
		common_term_product *= c ** (n - 1)
		terms_ += [(c, n)]
		coeff_product *= c

	term_sum = 0
	for (c, n) in terms_:
		term_sum += coeff_product / c * n

	return common_term_product * term_sum

def ad2 (x):
	''' Computes the arithmetic derivative of x (second impl) '''
	return sum(map(lambda (c, n): n * x / c, factorize(x)))


if __name__ == '__main__':
	sum_ = 0
	for x in range(1, 5000):
		c = ad_gcd3(x)
		sum_ += c
		print(x, ad(x), ad2(x), c, sum_)
	print(sum_)

primes.hpp

#ifndef __PRIMES_HPP__
#define __PRIMES_HPP__
#include <vector>
#include <cmath>

// Generates a list of prime numbers from 1 to n using the sieve of Eratosthenes (sieve implemented
// as a bit array to save space).
std::vector<int64_t> generate_primes (int64_t n) {
	// Note: not the most efficient impl of a prime sieve (eg. sieve could just store odds, and skip 1 and 2, 
	// *but* would require fancy indexing arithmetic that could cause bugs and add computational overhead).
	
	std::vector<bool> sieve (n);
	std::vector<int64_t> primes { 2 };
	std::fill(sieve.begin(), sieve.end(), true);

	// int64_t max = sqrt(n);
	for (int64_t i = 3; i < n; i += 2) {
		if (sieve[i] == true) {
			// std::cout << i << std::endl;
			primes.push_back(i);
			for (auto j = i*2; j < n; j += i) {
				// std::cout << "j = " << j << "\n";
				sieve[j] = false;
			}
			// std::cout << join(", ", primes) << std::endl;
		}
	}
	return primes;
}

#endif