% Largest triangle given n points on a unit circle
% Writing efficient code

% Part 1
  idx=1;
  for n=100:50:200
      theta = rand(n,1)*2*pi;  % Angle of random pts on the unit circle
      % Determine how long it takes to compute the largest
      % possible triangle obtained by selecting vertices from the
      % points represented by theta
      tic;

      A = 0;
      for i=1:n
          for j=1:n
              for k=1:n
                  % theta --> Cartesian
                  c1 = cos(theta(i));
                  s1 = sin(theta(i));
                  c2 = cos(theta(j));
                  s2 = sin(theta(j));
                  c3 = cos(theta(k));
                  s3 = sin(theta(k));
                  % Area using Heron's Formula
                  a = sqrt((c1-c2)^2 + (s1-s2)^2);
                  b = sqrt((c1-c3)^2 + (s1-s3)^2);
                  c = sqrt((c2-c3)^2 + (s2-s3)^2);
                  s = (a+b+c)/2;
                  Aijk = sqrt((s-a)*(s-b)*(s-c)*s);
                  A = max(A,Aijk);
              end
          end
      end
      t1(idx)=toc;
      idx=idx+1;
  end
  
% Part 2

  idx=1;
  for n=100:50:200
      theta = rand(n,1)*2*pi;  % Angle of random pts on the unit circle
      % Determine how long it takes to compute the largest
      % possible triangle obtained by selecting vertices from the
      % points represented by theta
      tic;

      A = 0;
      for i=1:n-2
          c1 = cos(theta(i));
          s1 = sin(theta(i));
          for j=i+1:n-1
              c2 = cos(theta(j));
              s2 = sin(theta(j));
              for k=j+1:n
                  % theta --> Cartesian
                  c3 = cos(theta(k));
                  s3 = sin(theta(k));
                  % Area using Heron's Formula
                  a = sqrt((c1-c2)^2 + (s1-s2)^2);
                  b = sqrt((c1-c3)^2 + (s1-s3)^2);
                  c = sqrt((c2-c3)^2 + (s2-s3)^2);
                  s = (a+b+c)/2;
                  Aijk = sqrt((s-a)*(s-b)*(s-c)*s);
                  A = max(A,Aijk);
              end
          end
      end
      t2(idx)=toc;
      idx=idx+1;
  end

% Part 3

  idx=1;
  for n=100:50:200
      theta = rand(n,1)*2*pi;  % Angle of random pts on the unit circle
      % Determine how long it takes to compute the largest
      % possible triangle obtained by selecting vertices from the
      % points represented by theta
      tic;

      cx = cos(theta);
      cy = sin(theta);
      D=zeros(n,n);
      for i=1:n
          for j=i:n
              Dij = sqrt((cx(i)-cx(j))^2+(cy(i)-cy(j))^2);
              D(i,j) = Dij;
              D(j,i) = Dij;
          end
      end
      
      A = 0;
      for i=1:n-2
          for j=i+1:n-1
              a = D(i,j);
              for k=j+1:n
                  % theta --> Cartesian
                  % Area using Heron's Formula
                  b = D(i,k);
                  c = D(j,k);
                  s = (a+b+c)/2;
                  Aijk = sqrt((s-a)*(s-b)*(s-c)*s);
                  A = max(A,Aijk);
              end
          end
      end
      t3(idx)=toc;
      idx=idx+1;
  end


% Part 4

n=100:50:200;
plot(n,t1,n,t2,n,t3);
title('Execution time vs. n');
xlabel('n');
ylabel('Running time in seconds');
legend('Method 1', 'Method 2', 'Method 3');


disp('  n      t1/t2     t1/t3')
disp('---------------------------')
idx=1;
for n=100:50:200
    disp(sprintf('%3d %8.1f %10.1f', n, t1(idx)/t2(idx), ...
                                        t1(idx)/t3(idx)))
    idx=idx+1;
end

disp(' ')
disp('The speed-up you get going from method 1 (t1) to method 2 (t2)')
disp('  is the ratio t2/t1, which is about 1/8 (8 times faster).')

disp(' ')
disp('What if n=1000?  You can try to extrapolate from the graphs.')
disp('  Another thought is to try to determine the relationship between t and n.')
disp('  Since there are triply nested loops, each with a range of about n,')
disp('  maybe we can try to relate t to n^3, resulting in the table below.')
disp('  The multiplier 10e6 just makes the ratios easier to read. ')
disp('  These numbers are MACHINE DEPENDENT.')

disp(' ')
disp('  n      t1/n^3*10^6     t2/n^3*10^6     t3/n^3*10^6');
disp('------------------------------------------------------------------');
idx=1;
for n=100:50:200
    disp(sprintf(' %d       %f       %f        %f', n,t1(idx)/n^3*10^6, ...
                                                      t2(idx)/n^3*10^6, ...
                                                      t3(idx)/n^3*10^6));
    idx=idx+1;
end

disp(' ')
disp('You should see that t1/n^3 is similar for all n.  The same can ')
disp('  be said for t2/n^3 and t3/n^3.  You can use these ratios to ')
disp('  estimate t for n=1000.')
disp('E.g., On a particular machine t1/n^3, t2/n^3, and t3/n^3 are ')
disp('  about 2*10e-6, 0.3*10e-6, and 0.07*10e-6, respectively.  Then ')
disp('  I would expect that the required time for n=1000 for the 3 ')
disp('  methods to be about 2000s, 300s, and 70s.  Just as a sanity check, ')
disp('  I ran method 3 for n=1000 ten times--t was between 65 and 67 ')
disp('  seconds for all 10 trials.')

