// (c) Arndt Brünner, alle Rechte vorbehalten!
//
// Berechnet die Schnittpunkte zweier Quadriken (Kegelschnitte)
// a1 x² + b1 y² + c1 x y + d1 y + e1 y + f1 = 0
// a2 x² + b2 y² + c2 x y + d2 y + e2 y + f2 = 0
// Benötigt polysolv.js
// In X und Y werden Arrays mit den Koordinaten zurückgegeben

Erl=false;
function kegelschnittschnitte(a1,b1,c1,d1,e1,f1,a2,b2,c2,d2,e2,f2,X,Y)
{
	if((a1==0)&&(a2==0)&&(b1==0)&&(b2==0))return kegelschnittschnitte4(a1,b1,c1,d1,e1,f1,a2,b2,c2,d2,e2,f2,X,Y);
	if((a2==0)&&(b2==0)&&(a1!=0))return kegelschnittschnitte6(a1,b1,c1,d1,e1,f1,a2,b2,c2,d2,e2,f2,X,Y);
	if((a1==0)&&(b1==0)&&(a2!=0))return kegelschnittschnitte6(a2,b2,c2,d2,e2,f2,a1,b1,c1,d1,e1,f1,X,Y);
	if(((a1==0)&&(b2==0))||((a2==0)&&(b1==0)))return kegelschnittschnitte3(a1,b1,c1,d1,e1,f1,a2,b2,c2,d2,e2,f2,X,Y);
	if(a1*a2==0)return kegelschnittschnitte2(a1,b1,c1,d1,e1,f1,a2,b2,c2,d2,e2,f2,X,Y);

	// für a1!=0 und a2!=0:

	//var a12=a1*a1,a13=a12*a1,a22=a2=a2,c22=c2*c2,d22=d2*d2,e22=e2*e2,e12=e1*e1,c12=c1*c1,d12=d1*d1,f12=f1*f1,i=0;
	var A,B,C,D,E,i=0;
	A=(b2*b2)*(a1*a1*a1)+b1*(c2*c2)*(a1*a1)-2*a2*b1*b2*(a1*a1)-b2*c1*c2*(a1*a1)+(a2*a2)*(b1*b1)*a1+a2*b2*(c1*c1)*a1-a2*b1*c1*c2*a1;
	B=2*b2*e2*(a1*a1*a1)-b2*c2*d1*(a1*a1)-b2*c1*d2*(a1*a1)+2*b1*c2*d2*(a1*a1)+(c2*c2)*e1*(a1*a1)-2*a2*b2*e1*(a1*a1)-2*a2*b1*e2*(a1*a1)-c1*c2*e2*(a1*a1)+2*a2*b2*c1*d1*a1-a2*b1*c2*d1*a1-a2*b1*c1*d2*a1+2*(a2*a2)*b1*e1*a1-a2*c1*c2*e1*a1+a2*(c1*c1)*e2*a1;
	C=(e2*e2)*(a1*a1*a1)+2*b2*f2*(a1*a1*a1)+b1*(d2*d2)*(a1*a1)-b2*d1*d2*(a1*a1)+2*c2*d2*e1*(a1*a1)-c2*d1*e2*(a1*a1)-c1*d2*e2*(a1*a1)-2*a2*e1*e2*(a1*a1)+(c2*c2)*f1*(a1*a1)-2*a2*b2*f1*(a1*a1)-2*a2*b1*f2*(a1*a1)-c1*c2*f2*(a1*a1)+a2*b2*(d1*d1)*a1+(a2*a2)*(e1*e1)*a1-a2*b1*d1*d2*a1-a2*c2*d1*e1*a1-a2*c1*d2*e1*a1+2*a2*c1*d1*e2*a1+2*(a2*a2)*b1*f1*a1-a2*c1*c2*f1*a1+a2*(c1*c1)*f2*a1;
	D=2*e2*f2*(a1*a1*a1)+(d2*d2)*e1*(a1*a1)-d1*d2*e2*(a1*a1)+2*c2*d2*f1*(a1*a1)-2*a2*e2*f1*(a1*a1)-c2*d1*f2*(a1*a1)-c1*d2*f2*(a1*a1)-2*a2*e1*f2*(a1*a1)-a2*d1*d2*e1*a1+a2*(d1*d1)*e2*a1-a2*c2*d1*f1*a1-a2*c1*d2*f1*a1+2*(a2*a2)*e1*f1*a1+2*a2*c1*d1*f2*a1;
	E=(f2*f2)*(a1*a1*a1)+(d2*d2)*f1*(a1*a1)-d1*d2*f2*(a1*a1)-2*a2*f1*f2*(a1*a1)+(a2*a2)*(f1*f1)*a1-a2*d1*d2*f1*a1+a2*(d1*d1)*f2*a1;
	//alert(A+"\n"+B+"\n"+C+"\n"+C+"\n"+D+"\n"+E);


	var n=QQuad(1,B/A,C/A,D/A,E/A);
	//alert(x1r+" "+x1i+"i\n"+x2r+" "+x2i+"i\n"+x3r+" "+x3i+"i\n"+x4r+" "+x4i+"i\n");
	if(Math.abs(x1i)<1e-9)Y[i++]=x1r;
	if(Math.abs(x2i)<1e-9)Y[i++]=x2r;
	if((n>2)&&(Math.abs(x3i)<1e-9))Y[i++]=x3r;
	if((n>3)&&(Math.abs(x4i)<1e-9))Y[i++]=x4r;
	var t="";
	for(i=0;i<Y.length;i++)
	{
		if((Y[i]==0)&&(a1!=0))
		{
			X[i++]=(Math.sqrt((d1*d1)-4*a1*f1)-d1)/(2*a1)
			X[i]=-(Math.sqrt((d1*d1)-4*a1*f1)+d1)/(2*a1);
		}
		else if((Y[i]*(a1*c2-a2*c1)+a1*d2-a2*d1)==0)
		{
			X[i++]=(Math.sqrt(-a*(4*(a1*a1)*(b1*(d2*d2)+c2*(c2*f1-d2*e1))-a1*(4*a2*(2*b1*d1*d2+c1*(2*c2*f1-d2*e1)-c2*d1*e1)+((c1*d2-c2*d1)*(c1*d2-c2*d1)))+4*(a2*a2)*(b1*(d1*d1)+c1*(c1*f1-d1*e1))))+a1*(c2*d1-c1*d2))/(2*a1*(a2*c1-a1*c2));
			X[i]=(Math.sqrt(-a*(4*(a1*a1)*(b1*(d2*d2)+c2*(c2*f1-d2*e1))-a1*(4*a2*(2*b1*d1*d2+c1*(2*c2*f1-d2*e1)-c2*d1*e1)+((c1*d2-c2*d1)*(c1*d2-c2*d1)))+4*(a2*a2)*(b1*(d1*d1)+c1*(c1*f1-d1*e1))))+a1*(c1*d2-c2*d1))/(2*a1*(a1*c2-a2*c1));
		}
		else  if((a1!=0)||(c1!=0))
			X[i]=-((Y[i]*Y[i])*(a1*b2-a2*b1)+Y[i]*(a1*e2-a2*e1)+a1*f2-a2*f1)/(Y[i]*(a1*c2-a2*c1)+a1*d2-a2*d1);
		else X[i]=-(b1*(Y[i]*Y[i])+e1*Y[i]+f1)/d1;
		t+=X[i]+"   "+Y[i]+"\n";
	}
	//alert("(1) Schnittpunkte:\n"+t);
}
function kegelschnittschnitte2(a1,b1,c1,d1,e1,f1,a2,b2,c2,d2,e2,f2,X,Y)
{
	if(b1*b2==0)return kegelschnittschnitte5(a1,b1,c1,d1,e1,f1,a2,b2,c2,d2,e2,f2,X,Y);

	// für b1!=0 und b2!=0:

	var A=b1*((a1*a1)*(b2*b2)+a1*(c2*(b1*c2-b2*c1)-2*a2*b1*b2)+a2*(a2*(b1*b1)-c1*(b1*c2-b2*c1)));
	var B=-b1*(a1*(2*b1*(b2*d2-c2*e2)-b2*(2*b2*d1-c1*e2-c2*e1))-a2*(2*(b1*b1)*d2-b1*(2*b2*d1+c1*e2+c2*e1)+2*b2*c1*e1)+(b1*c2-b2*c1)*(c1*d2-c2*d1));
	var C=-b1*(a1*(b1*(2*b2*f2-(e2*e2))-2*(b2*b2)*f1+b2*e1*e2)-a2*(2*(b1*b1)*f2-b1*(2*b2*f1+e1*e2)+b2*(e1*e1))-(b1*b1)*(d2*d2)+b1*(2*b2*d1*d2+c1*(c2*f2+d2*e2)-c2*(c2*f1+2*d1*e2-d2*e1))-b2*(b2*(d1*d1)+(c1*c1)*f2-c1*(c2*f1+d1*e2-2*d2*e1)-c2*d1*e1));
	var D=(2*(b1*b1*b1)*d2*f2-(b1*b1)*(2*b2*(d1*f2+d2*f1)+c1*e2*f2+c2*(e1*f2-2*e2*f1)-e2*(d1*e2-d2*e1))+b1*b2*(2*b2*d1*f1+c1*(2*e1*f2-e2*f1)-c2*e1*f1-d1*e1*e2+d2*(e1*e1)));
	var E=(b1*b1*b1)*(f2*f2)-(b1*b1)*(2*b2*f1*f2+e2*(e1*f2-e2*f1))+b1*b2*(b2*(f1*f1)+(e1*e1)*f2-e1*e2*f1);
/*
	var A=b1*(a1*c2-a2*c1)*(b1*c2-b2*c1);
	var B=(a1*b1*(2*b1*c2*e2-b2*(c1*e2+c2*e1))+a2*b1*(2*b2*c1*e1-b1*(c1*e2+c2*e1))+c1*(b1*c2-b2*c1)*(b2*d1-b1*d2));
	var C=(a1*b1*(b1*(e2*e2)+b2*(1-e1*e2))+a2*b1*(b2*(e1*e1)-b1*(e1*e2+1))+(b1*b1)*(c1*(c2*(e2-f2)-d2*e2)+c2*e1*(c2*f1-d2))+b1*(c2*d1-b2*((c1*c1)*(e2-f2)+c1*(c2*(e1*(2*f1+1)-f1)-d1*e2-2*d2*e1)-c2*d1*e1))+b2*c1*(b2*(c1*(e1*(f1+1)-f1)-2*d1*e1)-d1));
	var D=((b1*b1)*(c1*e2*(e2-f2)+c2*e1*(2*e2*f1-f2)-d2*(e1*e2+1))+b1*(d1*e2-b2*(c1*(2*e1*(e2*(f1+1)-f2)-e2*f1)+c2*e1*f1*(2*e1-1)-d1*(e1*e2+1)-d2*(e1*e1)))+b2*e1*(b2*(c1*(e1*(2*f1+1)-2*f1)-d1*e1)-d1));
	var E=(b1*b1)*(e1*e2*(e2*f1-f2)-f2)-b1*b2*((e1*e1)*(2*e2*f1-f2)-e1*e2*f1-f1)+(b2*b2)*(e1*e1)*f1*(e1-1);
*/
	var n=QQuad(A,B,C,D,E);
	//alert(x1r+" "+x1i+"i\n"+x2r+" "+x2i+"i\n"+x3r+" "+x3i+"i\n"+x4r+" "+x4i+"i\n");
	var i=0;
	if(Math.abs(x1i)<1e-9)X[i++]=x1r;
	if(Math.abs(x2i)<1e-9)X[i++]=x2r;
	if((n>2)&&(Math.abs(x3i)<1e-9))X[i++]=x3r;
	if((n>3)&&(Math.abs(x4i)<1e-9))X[i++]=x4r;
	var t="";
	for(i=0;i<X.length;i++)
	{
		if((X[i]==0)&&(b1!=0))
		{
			Y[i++]=(Math.sqrt((e1*e1)-4*b1*f1)-e1)/(2*b1)
			Y[i]=-(Math.sqrt((e1*e1)-4*b1*f1)+e1)/(2*b1);
		}
		else if((X[i]*(b1*c2-b2*c1)+b1*e2-b2*e1)==0)
		{
			Y[i++]=(Math.sqrt(-b*(4*(b1*b1)*(a1*(e2*e2)+c2*(c2*f1-e2*d1))-b1*(4*b2*(2*a1*e1*e2+c1*(2*c2*f1-e2*d1)-c2*e1*d1)+((c1*e2-c2*e1)*(c1*e2-c2*e1)))+4*(b2*b2)*(a1*(e1*e1)+c1*(c1*f1-e1*d1))))+b1*(c2*e1-c1*e2))/(2*b1*(b2*c1-b1*c2));
			Y[i]=(Math.sqrt(-b*(4*(b1*b1)*(a1*(e2*e2)+c2*(c2*f1-e2*d1))-b1*(4*b2*(2*a1*e1*e2+c1*(2*c2*f1-e2*d1)-c2*e1*d1)+((c1*e2-c2*e1)*(c1*e2-c2*e1)))+4*(b2*b2)*(a1*(e1*e1)+c1*(c1*f1-e1*d1))))+b1*(c1*e2-c2*e1))/(2*b1*(b1*c2-b2*c1));
			//alert(-b*(4*(b1*b1)*(a1*(e2*e2)+c2*(c2*f1-e2*d1))-b1*(4*b2*(2*a1*e1*e2+c1*(2*c2*f1-e2*d1)-c2*e1*d1)+((c1*e2-c2*e1)*(c1*e2-c2*e1)))+4*(b2*b2)*(a1*(e1*e1)+c1*(c1*f1-e1*d1))));
		}
		else if((b2!=0)||(c2!=0)) 
			{Y[i]=-((X[i]*X[i])*(b1*a2-b2*a1)+X[i]*(b1*d2-b2*d1)+b1*f2-b2*f1)/(X[i]*(b1*c2-b2*c1)+b1*e2-b2*e1);
			}
		else if(e2!=0)Y[i]=-(a2*(X[i]*X[i])+d2*X[i]+f2)/e2;
		else alert("?");
		t+=X[i]+"   "+Y[i]+"\n";
	}
	//alert("(2) Schnittpunkte:\n"+t);
}
function kegelschnittschnitte3(a1,b1,c1,d1,e1,f1,a2,b2,c2,d2,e2,f2,X,Y)
{
	if((a2==0)&&(b1==0))return kegelschnittschnitte3(a2,b2,c2,d2,e2,f2,a1,b1,c1,d1,e1,f1,X,Y);
	var c14=c1*c1*c1*c1,A,B,C,D,E;
	if(c1!=0)
	{
		A=b1*c14*(a2*b1-c1*c2);
		B=c14*(2*a2*b1*e1-b1*(c1*d2+c2*d1)+c1*(c1*e2-c2*e1));
		C=c14*(a2*(2*b1*f1+(e1*e1))-b1*d1*d2+(c1*c1)*f2-c1*(c2*f1-2*d1*e2+d2*e1)-c2*d1*e1);
		D=c14*(2*a2*e1*f1+c1*(2*d1*f2-d2*f1)-d1*(c2*f1-d1*e2+d2*e1));
		E=c14*(a2*(f1*f1)+d1*(d1*f2-d2*f1));
	}
	else if(d1!=0)
	{
		A=a2*(b1*b1)/(d1*d1);
		B=b1*(2*a2*e1-c2*d1)/(d1*d1);
		C=(a2*(2*b1*f1+(e1*e1))-b1*d1*d2-c2*d1*e1)/(d1*d1);
		D=(2*a2*e1*f1-c2*d1*f1+(d1*d1)*e2-d1*d2*e1)/(d1*d1);
		E=(a2*(f1*f1)+(d1*d1)*f2-d1*d2*f1)/(d1*d1);
	}
	else 
	{
		Y[0]=Y[1]=(Math.sqrt((e1*e1)-4*b1*f1)-e1)/(2*b1);
		Y[2]=Y[3]=-(Math.sqrt((e1*e1)-4*b1*f1)+e1)/(2*b1);
		X[0]=(Math.sqrt(2)*Math.sqrt(Math.sqrt((e1*e1)-4*b1*f1)*(c2*(2*b1*d2-c2*e1)-4*a2*b1*e2)+4*a2*b1*(e1*e2-2*b1*f2)+2*(b1*b1)*(d2*d2)-2*b1*c2*(c2*f1+d2*e1)+(c2*c2)*(e1*e1))-c2*Math.sqrt((e1*e1)-4*b1*f1)-2*b1*d2+c2*e1)/(4*a2*b1);
		X[1]=-(Math.sqrt(2)*Math.sqrt(Math.sqrt((e1*e1)-4*b1*f1)*(c2*(2*b1*d2-c2*e1)-4*a2*b1*e2)+4*a2*b1*(e1*e2-2*b1*f2)+2*(b1*b1)*(d2*d2)-2*b1*c2*(c2*f1+d2*e1)+(c2*c2)*(e1*e1))+c2*Math.sqrt((e1*e1)-4*b1*f1)+2*b1*d2-c2*e1)/(4*a2*b1);
		X[2]=-(Math.sqrt(2)*Math.sqrt(Math.sqrt((e1*e1)-4*b1*f1)*(4*a2*b1*e2-c2*(2*b1*d2-c2*e1))+4*a2*b1*(e1*e2-2*b1*f2)+2*(b1*b1)*(d2*d2)-2*b1*c2*(c2*f1+d2*e1)+(c2*c2)*(e1*e1))-c2*Math.sqrt((e1*e1)-4*b1*f1)+2*b1*d2-c2*e1)/(4*a2*b1);
		X[3]=(Math.sqrt(2)*Math.sqrt(Math.sqrt((e1*e1)-4*b1*f1)*(4*a2*b1*e2-c2*(2*b1*d2-c2*e1))+4*a2*b1*(e1*e2-2*b1*f2)+2*(b1*b1)*(d2*d2)-2*b1*c2*(c2*f1+d2*e1)+(c2*c2)*(e1*e1))+c2*Math.sqrt((e1*e1)-4*b1*f1)-2*b1*d2+c2*e1)/(4*a2*b1);
		return;
	}
	var n=QQuad(A,B,C,D,E);
	//alert("(3)"+"\n"+A+"\n"+B+"\n"+C+"\n"+D+"\n"+E+"\n\n"+x1r+" "+x1i+"i\n"+x2r+" "+x2i+"i\n"+x3r+" "+x3i+"i\n"+x4r+" "+x4i+"i\n");
	var i=0;
	if(Math.abs(x1i)<1e-9)Y[i++]=x1r;
	if(Math.abs(x2i)<1e-9)Y[i++]=x2r;
	if((Math.abs(x3i)<1e-9)&&(n>2))Y[i++]=x3r;
	if((Math.abs(x4i)<1e-9)&&(n>3))Y[i++]=x4r;
	var t="";
	for(i=0;i<Y.length;i++)
	{
		if(c1!=0)X[i]=-(b1*(Y[i]*Y[i])+e1*Y[i]+f1)/(c1*Y[i]+d1);
		else {X[i]=-(b1*(Y[i]*Y[i])+e1*Y[i]+f1)/d1;}
		t+=X[i]+"   "+Y[i]+"\n";
	}
	//alert("(3) Schnittpunkte:\n"+t);
}
function kegelschnittschnitte4(a1,b1,c1,d1,e1,f1,a2,b2,c2,d2,e2,f2,X,Y)
{
	//alert("(4) Schnittpunkte:");
	if((c1==0)&&(c2==0))
	{
		//zwei Geraden
		X[0]=X[1]=(e1*f2-e2*f1)/(d1*e2-d2*e1);Y[0]=Y[1]=(d2*f1-d1*f2)/(d1*e2-d2*e1);
		return;
	}
	var r=(c1*c1)*(f2*f2)-2*c1*(c2*f1*f2+d1*e2*f2+d2*(e1*f2-2*e2*f1))+(c2*c2)*(f1*f1)+2*c2*(d1*(2*e1*f2-e2*f1)-d2*e1*f1)+((d1*e2-d2*e1)*(d1*e2-d2*e1));
	if(r<0)return;r=Math.sqrt(r);
	X[0]=(r+c1*f2-c2*f1-d1*e2+d2*e1)/(2*(c2*d1-c1*d2));
	Y[0]=(r-c1*f2+c2*f1-d1*e2+d2*e1)/(2*(c1*e2-c2*e1));
	if(r==0)return;
	X[1]=(r-c1*f2+c2*f1+d1*e2-d2*e1)/(2*(c1*d2-c2*d1));
	Y[1]=(r+c1*f2-c2*f1+d1*e2-d2*e1)/(2*(c2*e1-c1*e2));
}
function kegelschnittschnitte5(a1,b1,c1,d1,e1,f1,a2,b2,c2,d2,e2,f2,X,Y)
{
	if((a1!=0)||(a2!=0))alert("Fehler");
	if(b2!=0)return kegelschnittschnitte5(a2,b2,c2,d2,e2,f2,a1,b1,c1,d1,e1,f1,X,Y);
	var n,i=0,A,B,C,D;
	if(c2!=0)
	{
		A=b1*(c2*c2*c2);
		B=(c2*c2)*(b1*d2-c1*e2+c2*e1);
		C=-(c2*c2)*(c1*f2-c2*f1+d1*e2-d2*e1);	
		D=(c2*c2)*(d2*f1-d1*f2);
		n=Kub(A,B,C,D);
	}
	else if(d2!=0)
		n=Quad((b1*d2-c1*e2)/d2,-(c1*f2+d1*e2-d2*e1)/d2,-(d1*f2-d2*f1)/d2);
	if(Math.abs(x1i)<1e-9)Y[i++]=x1r;
	if(Math.abs(x2i)<1e-9)Y[i++]=x2r;
	if((Math.abs(x3i)<1e-9)&&(n>2))Y[i++]=x3r;
	if((Math.abs(x4i)<1e-9)&&(n>3))Y[i++]=x4r;
	var t="";
	for(i=0;i<Y.length;i++)
	{
		X[i]=-(e2*Y[i]+f2)/(c2*Y[i]+d2);
		t+=X[i]+"   "+Y[i]+"\n";
	}
	//alert("(5) Schnittpunkte:\n"+t);
}
function kegelschnittschnitte6(a1,b1,c1,d1,e1,f1,a2,b2,c2,d2,e2,f2,X,Y)
{
	//für a1!=0 und a2=b2=0:
	var A,B,C,D,E,n,i=0;
	if(c2!=0)
	{
		A=b1*c2*c2*c2*c2;
		B=(c2*c2*c2)*(2*b1*d2-c1*e2+c2*e1);
		C=(c2*c2)*(a1*(e2*e2)+b1*(d2*d2)-c1*(c2*f2+d2*e2)+c2*(c2*f1-d1*e2+2*d2*e1));
		D=(c2*c2)*(2*a1*e2*f2-c1*d2*f2+c2*(2*d2*f1-d1*f2)-d2*(d1*e2-d2*e1));
		E=(c2*c2)*(a1*(f2*f2)-d2*(d1*f2-d2*f1));
		n=QQuad(A,B,C,D,E);
	}
	else if(d2!=0)
	{
		A=(a1*(e2*e2)+b1*(d2*d2)-c1*d2*e2);
		B=(2*a1*e2*f2-c1*d2*f2-d1*d2*e2+(d2*d2)*e1);
		C=a1*(f2*f2)-d1*d2*f2+(d2*d2)*f1;
		n=Quad(A,B,C);
	}
	else{ alert("Fehler: In Quadrik kommt eine Variable nicht vor");return;}
	if(Math.abs(x1i)<1e-9)Y[i++]=x1r;
	if(Math.abs(x2i)<1e-9)Y[i++]=x2r;
	if((Math.abs(x3i)<1e-9)&&(n>2))Y[i++]=x3r;
	if((Math.abs(x4i)<1e-9)&&(n>3))Y[i++]=x4r;
	var t="";
	for(i=0;i<Y.length;i++)
	{
		if(c2==0)X[i]=-(e2*Y[i]+f2)/d2;
		else X[i]=-(e2*Y[i]+f2)/(c2*Y[i]+d2);
		t+=X[i]+"   "+Y[i]+"\n";
	}
	//alert("(6) Schnittpunkte:\n"+t);
}