function [epsilonm]=machine_precision
epsilonm=1;
while 1+epsilonm> 1
    epsilonm=epsilonm/2;
end
epsilonm=epsilonm*2;