function daenhoubutuB(aa,b,c,P)
  bb=P\b;
  cc=c-bb(1)^2/aa(1)-bb(2)^2/aa(2);
  a1=-2*bb(3)/aa(1);
  a2=-2*bb(3)/aa(2);
  xx0=[bb(1)/aa(1);bb(2)/aa(2);cc/2/bb(3)];
  n=15;
  zzzz=0:1/n:1;
  zzz=zzzz.*zzzz;
  tate=ones(2*n+1,1);
  zzz=tate*zzz;
  if a1 < 0
    zzz=-zzz;
  end
  tt=(0:pi/n:2*pi)';
  ct=cos(tt);
  st=sin(tt);
  xxxx=zzzz*sqrt(abs(a1));
  yyyy=zzzz*sqrt(abs(a2));
  xxx=ct*xxxx;
  yyy=st*yyyy;
  Q=[P zeros(3,1);zeros(1,3) 1];
  T=[eye(3) -xx0;zeros(1,3) 1];
  yoko=ones(1,n+1);
  for i=1:2*n+1
    xxxx=Q*T*[xxx(i,:);yyy(i,:);zzz(i,:);yoko];
    x(i,:)=xxxx(1,:);y(i,:)=xxxx(2,:);z(i,:)=xxxx(3,:);
  end
  mesh(x,y,z)
  axis equal
  grid on
end

ax=3;ay=3;az=3;hxy=0;hyz=-1;hzx=2*sqrt(2);bx=-6;by=-3*sqrt(2);bz=3*sqrt(2);c=-1; % ensyuu P.206 mondai 188.1
A=[ax hxy hzx;hxy ay hyz;hzx hyz az];
b=[bx;by;bz];
[P AA]=eig(A);
aa=diag(AA);
if abs(aa(1))<10^-14
  R=[0 0 1;0 1 0;1 0 0];
  aa=R*aa;
  P=P*R;
elseif abs(aa(2))<10^-14
  R=[1 0 0;0 0 1;0 1 0];
  aa=R*aa;
  P=P*R;
end
daenhoubutuB(aa,b,c,P)