problem012.m

function N = problem012(P)
%The sequence of triangle numbers is generated by adding the natural numbers. 
%So the 7th triangle number would be 1 + 2 + 3 + 4 + 5 + 6 + 7 = 28.
%We can see that 28 is the first triangle number to have over five divisors.
%What is the value of the first triangle number to have over five hundred (or P in general)divisors?
% The sub-function no_of_factors_S is defined.The function below_primes
% should be defined in the directory.
n = 1;
p1 = 1;%no_of_factors(n); here n = 1
p2 = 2;%no_of_factors(n+1); 
p = p1*p2/2;% a0 is 1 here. The original formula is p = p1*p2*a0/(a0+1)
Prim = primes_to_n(100000);% Matrix B will be having all primes required for calculation.
while p < P 
    if rem(n,2) == 0
        n_temp = n;
    else
        n_temp = n+1;
    end
    
    a0 = 0;
    while rem(n_temp,2) == 0
        n_temp = n_temp/2;
        a0 = a0+1;
    end
    p = p1*p2*a0/(a0+1);
    %keyboard
    n = n+1;
    p1 = p2;
    p2 = no_of_factors_S(n+1,Prim);
end

%by the time it comes out the the value of n changes to n+1 and hence the
%final ans will be
N = n*(n-1)/2;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%defining subfunction no_of_factors_S%%%%%%%%%%%%%
%The function below_primes should be defined in the domain
function s = no_of_factors_S(n,Prim)
% This function returns the no of all the possible factors of the given no
% including 1 and that no itself.
%The function prime_factors.m should be defined in the directory
B = Prim;
k = n;
i = 1;j=1;
A = [ ];
am = 1;
s = 1; p = 1;
while k > 1      
		if rem(k, B(1,i)) == 0%start of first if look
            if j == 1
                A = B(1,i);
                j = 0;
                k = k / B(1,i);
                continue
            end
            
            k = k / B(1,i);
            if B(1,i) ~= A(1,am)
                s = s*(p+1);
                p = 1;
                A = [ A B(1,i)];
                am = am+1;
            else
                p = p + 1;
            end
        else%else corresponding to first if
            i = i + 1;
        end
end
s = s*(p+1);% this step is for last iteration.