% �Ï•ªS1‚ÌŒvŽZ

hold on;

set(0,'defaultAxesFontSize',16);
set(0,'defaultTextFontSize',16);

resultall = [];

for seed=0:2

    rng(seed, 'twister');

    N = 1000000000;
    %N = 1000000;
    M = 0;
    k = 1;
    step = 1;

    count = [];
    result = [];

    for i=1:N
        x = rand;
        y = rand;
        if x*x + y*y <= 1
            M = M+1;
        end
        if rem(i,step) == 0
            fprintf('%d\t%f\t%f\n', i, M/i, 4.0*M/i);
            count = [count; i];
            result = [result; M/i];
            k = k+1;
        end
        if k == 10
            step = step * 10;
            k = 1;
        end
    end
    
    resultall = [resultall, result];
    %plot(log10(count), result);
end

set(0,'defaultAxesFontSize',16);
set(0,'defaultTextFontSize',16);

plot(log10(count), resultall(:,1), '-k');
plot(log10(count), resultall(:,2), '--k');
plot(log10(count), resultall(:,3), '-.k');

xlim([0 9]);
ylim([0.475 0.925]);
xticks([0 1 2 3 4 5 6 7 8 9]);
xticklabels({'10^0', '10^1', '10^2', '10^3', '10^4', '10^5', '10^6', '10^7', '10^8', '10^9'});
yticks([0.50 0.55 0.60 0.65 0.70 0.75 0.80 0.85 0.90]);
yticklabels({'0.50', '0.55', '0.60', '0.65', '0.70', '0.75', '0.80', '0.85', '0.90'});

sz = size(count);
analytical = pi/4 * ones(sz);
plot(log10(count), analytical, ':k');

%xlabel('iteration');
%ylabel('integral');
%legend('trial #1', 'trial #2', 'trial #3', 'analytical solution', 'Location', 'southeast');
xlabel('‰ñ�”');
ylabel('�Ï•ª’l');
legend('�”’l‰ð #1', '�”’l‰ð #2', '�”’l‰ð #3', '‰ð�Í‰ð', 'Location', 'southeast');

saveas(gcf, 'S1.eps', 'eps');

hold off;
