clear all
       s=20;
a=13;
b=1;
m=64;
for I=1:m
       n=a*s+b;
u=mod(n,m);
v=u/m;

        H= (256/27)*v*(1-v)^3;
        if v<=H
           s=u;
            n=a*s+b;
              u=mod(n,m);
              v=u/m;
                 Beta=v
             else 
       I=I+1;
       s=u;
   end
end