function [] = EllipticCurveFactoring(n,B,K)
% factor = EllipticCurveFactoring(n,B,K)
ECurveCount = 1;

while ECurveCount<=K
    %Generate elliptic curve E and point P:
    a = floor(n*rand); x = floor(n*rand); y = floor(n*rand);  
    b = mod(y^2 - x^3 - a*x,n);
    P = [x y];
    fprintf('\rRandomly generated point = P = (%d, %d).\r', x, y)
    fprintf('Randomly generated elliptic curve:  y^2 = x^3 + %dx + %d.\r', a, b)
    ECurveCount = ECurveCount+1; %update for next round
    %We compute B!P internally
    mult = factorial(B);
    Q = [inf inf];
    D = P;

    while mult > 0
    if  mod(mult,2) == 1
        %Q = EllipticCurvePointAdditionModp(Q, D, a, b, n);
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        if Q(1) == inf, SumPoints = D; return, end
        if D(1) == inf, SumPoints = Q; return, end
        if Q == D
            dx = mod(2*Q(2),n);
            dy = mod(3*Q(1)^2+a,n);
            if gcd(dx,n)==1
                m = mod(dy*ModInv(dx,n),n);
                x3 = mod(m^2 - Q(1) - D(1),n);
                y3 = mod(m*(Q(1) - x3) - Q(2),n);
                SumPoints = [x3 y3];
            else
                if gcd(dx,n)> 1 & gcd(dx,n)< n
                    fprintf('Nontrivial Factor d = %d found!\r', gcd(dx,n))
                    factor = gcd(dx,n);
                    return
                end
                m = inf;
                SumPoints = [inf inf];
            end
        else
            if gcd(D(1)-Q(1),n) == 1
                m = mod((D(2)-Q(2))*ModInv((D(1)-Q(1)),n),n);
                x3 = mod(m^2 - Q(1) - D(1),n);
                y3 = mod(m*(Q(1)-x3)-Q(2),n);
                SumPoints = [x3 y3];
            else
                d =gcd(D(1)-Q(1),n);
                if d> 1 & d < n
                   fprintf('Nontrivial Factor d = %d found!\r', d)
                   factor = d;
                   return
                end
                m = inf;
                SumPoints = [inf inf];
            end
        end
        Q = SumPoints;
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    end   
    mult = floor(mult/2)
    %D = EllipticCurvePointAdditionModp(D, D, a, b, p);
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    if D(1) == inf, SumPoints = D; return, end
    if D(1) == inf, SumPoints = D; return, end
    if D == D
        dx = mod(2*D(2),n);
        dy = mod(3*D(1)^2+a,n);
        if gcd(dx,n)==1
            m = mod(dy*ModInv(dx,n),n);
            x3 = mod(m^2 - D(1) - D(1),n);
            y3 = mod(m*(D(1) - x3) - D(2),n);
            SumPoints = [x3 y3];
        else
            if gcd(dx,n)> 1 & gcd(dx,n)< n
               fprintf('Nontrivial Factor d = %d found!\r', gcd(dx,n))
               factor = gcd(dx,n);
               return
            end
            m = inf;
            SumPoints = [inf inf];
        end
    else
        if gcd(D(1)-D(1),n) == 1
            m = mod((D(2)-D(2))*ModInv((D(1)-D(1)),n),n);
            x3 = mod(m^2 - D(1) - D(1),n);
            y3 = mod(m*(D(1)-x3)-D(2),n);
            SumPoints = [x3 y3];
        else
            d =gcd(D(1)-Q(1),n);
            if d> 1 & d < n
                fprintf('Nontrivial Factor d = %d found!\r', d)
                factor = d;
                return
            end
            m = inf;
            SumPoints = [inf inf];
        end
    end
    D = SumPoints;
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    end
    
end