function [x1, dd] = quadratic_prep(a)
[m,dd]=size(a); 
b = mean(a);
b = repmat(b,m,1);
c = a - b;
x1 = c;
d = dd;
for j = 1:dd
  d = d+1; x1(:,d) = x1(:,j).^2;
end
for j = 1:dd-1
  for k = j+1:dd
    d = d + 1; x1(:,d) = x1(:,j).*x1(:,k);
  end
end