% Inclusion Exclussion principle in MATLAB
%
% From the set {1,2,3, ..., 315} a number is selected
% at random. 

% Using MATLAB to count favorable outcomes, find the probability
% that selected number is 
% divisible with at least one of 3, 5, or 7. For example, the command
% >> a = 1:N; A =   sum(mod(a, 3) == 0); 
% will count how many numbers form {1,...,N} are divisible by 3.

%
%(b) compare this probability with a naive solution
% 1/3 + 1/5 + 1/7 - 1/15 -1/21 - 1/35 + 1/105, and show 
%that the naive solution is correct!
% 
% (c) Is the naive solution correct for {1,2,...,N} if
%     N = 316?
% (d)  Does the naive solution
% solves the problem for any other N from {289, 340}?
% Plot this probability for N in {289, 340}. Is there any symmetry in the 
% plot?
 
clear all;
close all;
probs=[]; n1=289; n2=340;
for N = n1:n2 
numbers = 1:N;
A =   sum(mod(numbers , 3) == 0);
B =   sum(mod(numbers , 5) == 0);
C =   sum(mod(numbers , 7) == 0) ;
AB = sum(mod(numbers , 3) == 0 & mod(numbers , 5) == 0);
AC = sum(mod(numbers , 3) == 0 & mod(numbers , 7) == 0);
BC = sum(mod(numbers , 5) == 0 & mod(numbers , 7) == 0);
ABC = sum(mod(numbers , 3) == 0 & ...
                      mod(numbers , 5) == 0 & mod(numbers ,7)==0);
favorable = A+B+C-AB-AC-BC+ABC;
prob = favorable/N;
probs = [probs prob];
end
plot((n1:n2), probs, 'ko-')
hold on
limit = 1/3+1/5+1/7 - 1/15-1/21-1/35+1/105; %0.542857
plot((n1:n2), limit * ones(1, n2-n1+1),'r-','linewidth',2)
axis tight
