Problem 95 - Amicable Chains

#include <vector>

std::vector<int> generatePrimes(int n) {
	bool *primeSieve = new bool[n]();
	for (int i = 2; i < n; i++) primeSieve[i] = true;
	for (int i = 2; i*i < n; i++)
		if (primeSieve[i])
			for (int j = i*i; j < n; j += i)
				primeSieve[j] = false;
	std::vector<int> primes;
	for (int i = 0; i < n; i++)
		if (primeSieve[i]) primes.push_back(i);
	delete[] primeSieve;
	return primes;
}

int num = n;	// make a copy of n to reduce 
int sigma = 1;
// Find the distinct prime factors of n up to sqrt(n)
for (auto p : primes) {
	if (p*p > n) break;
	if (n % p != 0) continue;
	int m = p;			// m = p^(a+1)
	while (n % m == 0)	// get the multiplicity of p
		m *= p;
	num /= (m/p);		// reduce n at each step
	sigma *= (m-1)/(p-1);
}
// If n is prime or has a prime factor > sqrt(n) it is a distinct prime factor with multiplicity of 1
if (num > 1) sigma *= (num+1);    
aliquot[n] = sigma - n;

enum class color {
	white,	// unvisited
	grey,	// visiting
	black	// visited
};

// initialize all vertices as white (unvisited)
std::vector<color> visited = std::vector<color>(N+1, color::white);

bool dfs(currVertex, vertices, visited, N, &cycle)
// return true if currVertex is part of a cycle, false if not

(1) If the current vertex has already been visited or is out of range
	(i) Return false

(2) If the current vertex is a dead-end (leads to "1")
	(i) Mark the current vertex as visited
	(ii) Return false

(2) If the current vertex has been seen before
	the current vertex is part of a cycle
	(i) Add the current vertex to cycle
	(ii) Mark the current vertex as visited
	(iii) Return true

(3) The current vertex is unvisited, mark it as visiting 

(4) Recursively call dfs on the current vertex's single neighbor
	(i) If this call returns true and the current vertex is still marked as visiting
		the current vertex is part of a cycle
		(a) Add the current vertex to cycle
		(b) Mark the current vertex as visited
		(c) Return true
	(i) Else
		the current vertex is not part of a cycle
		(a) Mark the current vertex as visited
		(b) Return false

int *aliquot;	// graph adjacency list
std::vector<int> cycle, longestCycle;
for (int n = 2; n < N+1; n++) {
	if (aliquot[n] == 1 || visited[n] == color::black)
		continue;
	dfs(n, aliquot, N+1, visited, cycle);
	if (cycle.size() > longestCycle.size())
		longestCycle = cycle;
	cycle.clear();
}

#include <vector>
#include <cmath>		// std::sqrt
#include <algorithm>	// std::min_element
#include <iostream>		// std::cout

// unvisited, visiting, visited
enum class color { white, grey, black };

std::vector<int> generatePrimes(int N) {
	bool *primeSieve = new bool[N]();
	for (int i = 2; i < N; i++) primeSieve[i] = true;
	for (int i = 2; i*i < N; i++)
		if (primeSieve[i])
			for (int j = i*i; j < N; j += i) primeSieve[j] = false;
	std::vector<int> primes;
	for (int i = 0; i < N; i++)
		if (primeSieve[i]) primes.push_back(i);
	delete[] primeSieve;
	return primes;
}

void getAliquotSums(int *aliquot, int N, const std::vector<int> &primes) {
	for (int n = 2; n < N; n++) {
		int num = n;	// make a copy of n to reduce 
		int sigma = 1;
		// Find the distinct prime factors of n up to sqrt(n)
		for (auto p : primes) {
			if (p*p > n) break;
			if (n % p != 0) continue;
			int m = p;			// m = p^(a+1)
			while (n % m == 0)	// get the multiplicity of p
				m *= p;
			num /= (m/p);		// reduce n with each prime factor
			sigma *= (m-1)/(p-1);
		}
		// If n is prime or has a prime factor > sqrt(n),
		// it is a distinct prime factor with multiplicity of 1
		if (num > 1) sigma *= (num+1);    
		aliquot[n] = sigma - n;
	}
}

bool dfs(int curr, int *vertices, int N, std::vector<color> &visited, std::vector<int> &cycle) {
	// Return true if the curr vertex is part of a cycle, false if not
    
	// vertex is out of range, already visited, or a dead-end
	if (curr >= N || visited[curr] == color::black) return false;
	if (vertices[curr] == 1) {  // dead-end
		visited[curr] = color::black;
		return false;
	}

	// vertex was seen before, cycle detected
	if (visited[curr] == color::grey) {
		cycle.push_back(curr);
		visited[curr] = color::black;
		return true;
	}

	// vertex is unvisited, mark as visiting
	visited[curr] = color::grey;

	// check the single neighbor of curr
	if (dfs(vertices[curr], vertices, N, visited, cycle) && visited[curr] == color::grey) {
		// the vertex that curr points to is part of a cycle
		// and curr belongs to that cycle
		cycle.push_back(curr);
		visited[curr] = color::black;
		return true;
	} else {
		// vertex is not part of a cycle, mark as visited
		visited[curr] = color::black;
		return false;
	}
}

int longestAmicableChain(int N) {
	// Generate all of the primes under sqrt(N)
	std::vector<int> primes = generatePrimes(static_cast<int>(std::sqrt(N)));

	// Get the aliquot sum for every number up to N
	int *aliquot = new int[N+1]();
	getAliquotSums(aliquot, N, primes);

	// Find every cycle and record the longest one
	std::vector<color> visited = std::vector<color>(N+1, color::white);
	std::vector<int> cycle, longestCycle;
	for (int n = 2; n < N+1; n++) {
		if (aliquot[n] == 1 || visited[n] == color::black)
			continue;
		dfs(n, aliquot, N+1, visited, cycle);
		if (cycle.size() > longestCycle.size())
			longestCycle = cycle;
		cycle.clear();
	}
	delete[] aliquot;
	if (longestCycle.empty()) return 0;	// no cycle found

	// Print chain
	std::cout << "Chain of length " << longestCycle.size() << " found:\n";
	for (auto it = longestCycle.rbegin(); it != longestCycle.rend(); ++it)
		std::cout << *it << "->";
	std::cout << longestCycle.back() << "\n";
	return *std::min_element(longestCycle.begin(), longestCycle.end());
}