function [price, RMSE]= BS_MC_call(S0,K,T,r,sigma,M)
//    stacksize('max')
    W=rand(1,M,'normal');
    S=S0*exp(sigma*sqrt(T)*W+(r-sigma^2/2)*T);
    C=exp(-r*T)*max(S-K,0);
    price=sum(C)/M;
    varest=sum((price-C).^2)/(M-1);
    RMSE=sqrt(varest)/sqrt(M);
endfunction
