y=1/sqrt(2)/2;
A=[y -y ;y y];
B=[y y;-y y];
a=[y; y];
b=[-y ; y];
N=10000;
x=grand(2,1,'def');
p=0.5;

for j=1:N
        x=bool2s(grand(1,1,'def')>p)*(A*x+a-B*x-b)+B*x+b;
        plot2d(x(1),x(2),0,style=5,rect=[-0.3,0.4,0.3,0.7]);
end