function [] = Problem23()
%ProjectEuler/Octave/Problem23.m
%Matthew Ellison
% Created: 03-22-19
%Modified: 03-28-19
%Find the sum of all the positive integers which cannot be written as the sum of two abundant numbers
%{
	Copyright (C) 2019  Matthew Ellison

	This program is free software: you can redistribute it and/or modify
	it under the terms of the GNU Lesser General Public License as published by
	the Free Software Foundation, either version 3 of the License, or
	(at your option) any later version.

	This program is distributed in the hope that it will be useful,
	but WITHOUT ANY WARRANTY; without even the implied warranty of
	MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
	GNU Lesser General Public License for more details.

	You should have received a copy of the GNU Lesser General Public License
	along with this program.  If not, see <https://www.gnu.org/licenses/>.
%}


MAX_NUM = 28123;	%The highest number that will be evaluated


%Setup the variables
divisorSums = [];
%Make sure every element has a 0 in it's location
for cnt = 1 : MAX_NUM
	divisorSums(end + 1) = 0;
end

%Start the timer
startTime = clock();

%Get the sum of the divisors of all numbers < MAX_NUM
for cnt = 1 : MAX_NUM
	div = getDivisors(cnt);
	if(size(div)(2) > 1)
		div(end) = [];
	end
	divisorSums(cnt) = sum(div);
end

%Get the abundant numbers
abund = [];
for cnt = 1 : size(divisorSums)(2)
	if(divisorSums(cnt) > cnt)
		abund(end + 1) = cnt;
	end
end

%Check if each number can be the sum of 2 abundant numbers and add to the sum if no
sumOfNums = 0;
for cnt = 1 : MAX_NUM
	if(~isSum(abund, cnt))
		printf("Added %d to sum\n", cnt)
		sumOfNums += cnt;
	end
end

%Stop the timer
endTime = clock();

%Print the results
printf("The answer is %d\n", sumOfNums)
printf("It took %f seconds to run this algorithm\n", etime(endTime, startTime))

end	%Problem23

function [divisors] = getDivisors(goalNumber)
	divisors = [];
	%Start by checking that the number is positive
	if(goalNumber <= 0)
		return;
	%If the number is 1 return just itself
	elseif(goalNumber == 1)
		divisors(end + 1) = 1;
		return;
	end

	%Start at 3 and loop through all numbers < sqrt(goalNumber) looking for a number that divides it evenly
	topPossibleDivisor = ceil(sqrt(goalNumber));
	possibleDivisor = 1;
	while(possibleDivisor <= topPossibleDivisor)
		%If you find one add it and the number it creates to the list
		if(mod(goalNumber, possibleDivisor) == 0)
			divisors(end + 1) = possibleDivisor;
			%Account for the possibility sqrt(goalNumber) being a divisor
			if(possibleDivisor != topPossibleDivisor)
				divisors(end + 1) = goalNumber / possibleDivisor;
			end
			%Take care of a few occations where a number was added twice
			if(divisors(end) == (possibleDivisor + 1))
				possibleDivisor += 1;
			end
		end
		possibleDivisor += 1;
	end

	%Sort the list before returning for neatness
	divisors = sort(divisors);
end

function [answer] = isSum(abund, num)
	sumOfNums = 0;
	%Pick a number for the first part of the sum
	for firstNum = 1 : size(abund)(2)
		%Pick a number for the second part of the sum
		for secondNum = firstNum : size(abund)(2)
			sumOfNums = abund(firstNum) + abund(secondNum);
			if(sumOfNums == num)
				answer = true;
				return;
			elseif(sumOfNums > num)
				break;
			end
		end
	end
	answer = false;
end

%{
Results: I let this run for 11 hours and did not come up with a result. I added some tracking and tested again and it seems like it is doing what it is supposed to and will come to the correct answer... eventually
However, for now, this remains unproven.
%}