// Javascript zum Auffinden von Lösungen (nichtlinearer) Gleichungssysteme
// (c) Arndt Brünner, Gelnhausen, 9.8.2003
// Version: 10.8.2003

/*
sin(x)=y
y=2x-z
cos(x)=a
a=3x-1
tan(x)=b
b=12c-1
*/
var Absolutwert=1==0;
function f_(i,x)
{
	return eval(Fkt[i]);
}

function findVars(t,V,x)
{
	var i=0,k=0;
	for(var j=0;j<t.length;j++)
	{
		t[j]=t[j].replace(/ /g,"").toLowerCase();
		if(t[j]=="")continue;
		v=t[j].split("=");
		for(k=0;k<i;k++)if(V[k]==v[0])break;
		if(k<i)continue;
		if(v.length==1)x[i]=Math.random()*20-10;
		else if(v[1].charAt(0)=="[")
		{
			var b=v[1].substring(1,v[1].length-1).split("...");
			b[0]=Number(b[0].replace(/,/,"."));
			if(b.length==1)b[1]=b[0];else b[1]=Number(b[1].replace(/,/,"."));
			x[i]=Math.random()*(b[1]-b[0])+b[0];
		}
		else x[i]=Number(v[1].replace(/,/,"."));
		V[i++]=v[0];
	}
	return i;
}

var Fkt=new Array();
var trigo=0,Grad_rad=Math.PI/180,Rad_grad=180/Math.PI;

function glsl()
{
	trigo=document.f.trigomodus.selectedIndex;
	grad_rad=(trigo==1)?Grad_rad:1;rad_grad=(trigo==1)?Rad_grad:1;
	var v=document.f.v.value.replace(/,/g,".").replace(/\n/g," ").replace(/ +/g," ").replace(/ += +/g,"=").toLowerCase();
	var F,FF=document.f.g.value.replace(/ /g,"");
	if(FF=="")return;
	FF=FF.replace(/exp/gi,"¿").replace(/pi/gi,"pi");
	if(FF.indexOf("x")>-1)v+=" x";
	if(FF.indexOf("y")>-1)v+=" y";
	if(FF.indexOf("z")>-1)v+=" z";
	F=FF.replace(/,/g,".").split("\n")
	var x=new Array(3),i,n,j,k=0,ii,xstart=new Array();
	var Vars=new Array();
 	n=findVars(v.split(" "),Vars,x);//FF=FF.replace(/¿/g,"exp");
	for(i=0;i<F.length;i++)
	{
		if(F[i]=="")continue;
		var MultErg=new Array(/\dx/,/\d\(/,/\)\(/,/x\(/,/\)x/,/pi\(/,/e_\(/,/\dpi/,/\de_/);
		for(j=0;j<MultErg.length;j++)
		do{ii=F[k].search(MultErg[j]);if(ii==-1)break;F[k]=F[k].substr(0,ii+1)+"*"+F[k].substring(ii+1,F[k].length);}while(true);
		if(trigo==0)F[i]=F[i].replace(/°/g,"*3.141592653589793238/180"); else F[i]=F[i].replace(/°/g,"");
		F[i]=F[i].replace(/pi/g,"3.141592653589793238");
		F[i]=F[i].replace(/ê/g,"2.718281828459045235");
		if(F[i].indexOf("=")>-1)F[i]=F[i].split("=").join("-(")+")";
		for(var oink=0;oink<F[i].length-1;oink++)
		{
			var Oink=F[i].substr(oink,2).toLowerCase();
			if(Oink.search(/\d[a-zäöü]/)>-1){F[i]=F[i].substr(0,oink+1)+"*"+F[i].substring(oink+1,F[i].length);oink++;}
		}
//		Fkt[k]=translate(F[i].toUpperCase()).replace(/Math/g,"math");
		Fkt[k]=translate(F[i].toUpperCase().replace(/E/g,"Ê")).replace(/Math/g,"math").replace(/Ê/g,"E");
		for(j=0;j<n;j++){Fkt[k]=eval("Fkt[k].replace(/"+Vars[j].toUpperCase()+"/g,\"_["+j+"]\");");}
		ii=Fkt[k].search(/[A-ZÄÖÜß]/);
		if(ii>-1){v+=" "+Fkt[k].charAt(ii);n=findVars(v.split(" "),Vars,x);i--;continue;}
		Fkt[k]=Fkt[k].toLowerCase().replace(/math/g,"Math").replace(/_/g,"x").replace(/¿/g,"Math.exp");
		FF+=Fkt[k]+"\n";
		
		var MultErg=new Array(/\dx/,/\d\(/,/\)\(/,/x\(/,/pi\(/,/e_\(/,/\dpi/,/\de_/);
		for(j=0;j<MultErg.length;j++)
		do{ii=Fkt[k].search(MultErg[j]);if(ii==-1)break;Fkt[k]=Fkt[k].substr(0,ii+1)+"*"+Fkt[k].substring(ii+1,Fkt[k].length);}while(true);
		k++;
	}
	//alert(Fkt.join("\n")+"\n\n"+Vars);
	var kk=k;
	if(n>k){for(i=k;i<n;i++)Fkt[i]="0";k=n;}
	//if(n!=k){alert("Die Anzahlen der Gleichungen und der Variablen müssen übereinstimmen!");return;}
	var t="";
	var tt="",ttt="";
	niter=document.f.niter[document.f.niter.selectedIndex].text;
	if(niter=="")niter=(n<5)?20:10; else niter=Number(niter);
	var nrepeat=document.f.nrepeat[document.f.nrepeat.selectedIndex].text;

	for(var ir=0;ir<nrepeat;ir++)
	{
		for(i=0;i<n;i++)xstart[i]=x[i];
		t="";tt="",nniter=0;
		var grund="";
		document.f.t.value=(ir+1)+". Durchlauf";
		if((nniter=solveGS(x,n,niter))<0)
			tt+=(grund="Der Algorithmus konvergiert für diese Startwerte nicht.");
		else
		{
			for(i=0;i<n;i++){t+=Vars[i]+" = "+RR(x[i])+"\n";if(isNaN(x[i]))tt="Keine Lösung gefunden!\n"+(grund="Ungültiges Funktionsargument während der "+nniter+". Iteration\nim "+(ir+1)+". Durchlauf aufgetreten")+"\n\n";}
			t+="\nProbe (die Funktionswerte müssen 0 sein): \n";
			for(i=0;i<kk;i++){var p;t+="f"+(i+1)+"("+Vars.join(",")+") = "+(p=f_(i,x))+"\n";if((tt=="")&&(Math.abs(p)>.1e-4))tt=(ir+1)+". Durchlauf: "+(grund="Nach "+niter+" Iterationen keine Lösung gefunden")+"\n              (siehe Probe)\n\n";}
			t+="\nStartwerte:\n";for(i=0;i<n;i++)t+=Vars[i]+"0 = "+xstart[i]+"\n";
			t=t.replace(/\./g,",").replace(/ ,/g," 0,").replace(/-,/g,"-0,");
			if(tt==""){tt="Lösung im "+(ir+1)+". Durchlauf nach "+nniter+" Iterationen gefunden:\n\n";break;}
			if(ir<nrepeat-1)
			{
				ttt+=(ir+1)+". Durchlauf\n  Startwerte:\n";
				for(i=0;i<n;i++)ttt+="    "+Vars[i]+"0 = "+xstart[i]+"\n";
				ttt+="  "+grund+"\n\n";
			}
		}
		findVars(v.split(" "),Vars,x);
		for(i=0;i<n;i++)if(x[i]!=xstart[i])break;
		if(i==n)break;
	}
	ttt=ttt.replace(/\./g,",").replace(/ ,/g," 0,").replace(/-,/g,"-0,").replace(/, /g,". ");
	document.f.t.value=tt+t+"\n\n"+ttt;
}

function RR(x)
{return(" "+String(Math.round(x*1e12)/1e12)).replace(/\./,",").replace(/-,/,"-0,").replace(/ ,/,"0,").replace(/ /g,"");}

function solveGS(x,n,schritte)
{
	var i,j,k,c=new Array(n*n), ff=new Array(n*n),s,ss=1e+300,d,ds=0,xx=new Array(n);
	for(i=0;i<schritte;i++)
	{
		for(j=0;j<n;j++)xx[j]=x[j];
		s=MakeFVektor(ff,x,n);
		if(s-ss>50)return -i-1;ss=s;
		JacobiMatrix(c,x,ff,n);//alert(c);
		MatrixInvert(c,n);
		//InverseMatrix(c,n);
		//alert("Inverse Jacobimatrix:\n"+c);
		for(j=0;j<n;j++)
		{
			for(k=0;k<n;k++){x[j]-=(d=c[j*n+k]*ff[k]);ds+=Math.abs(d);if(Absolutwert)x[j]=Math.abs(x[j]);}
			if(isNaN(x[j]))return i+1;
		}
		for(k=0;k<n;k++){if(Math.abs(xx[k]-x[k])>1e-14)break;}
		if(k==n)return i;
		if(ds<1e-15)return i+1;
	}
	return schritte;
}

function MakeFVektor(m,x,n)
{
	var i,s=0;
	for(i=0;i<n;i++){m[i]=f_(i,x);s+=Math.abs(m[i]);}
	return s;
}

function InverseMatrix(M_,n)
{
	var nn,a,b,c,d,e,f,g,h,i,j,k;
	a=M_[0];b=M_[1];c=M_[2];d=M_[3];
	if(n==2)
	{
		nn=a*d-b*c;
		M_[0]=d/nn;M_[1]=-b/nn;M_[2]=-c/nn;M_[3]=a/nn;
	}
	else if(n==3)
	{
		e=M_[4];f=M_[5];g=M_[6];h=M_[7];i=M_[8];
		nn=a*(e*i-f*h)+b*(f*g-d*i)+c*(d*h-e*g);
		M_[0]=(e*i-f*h)/nn;
		M_[1]=(c*h-b*i)/nn;
		M_[2]=(b*f-c*e)/nn;
		M_[3]=(f*g-d*i)/nn;
		M_[4]=(a*i-c*g)/nn;
		M_[5]=(c*d-a*f)/nn;
		M_[6]=(d*h-e*g)/nn;
		M_[7]=(b*g-a*h)/nn;
		M_[8]=(a*e-b*d)/nn;
	}
	else if(n==4)
	{
		var mm0=M_[0],mm1=M_[1],mm2=M_[2],mm3=M_[3],mm4=M_[4];
		var mm5=M_[5],mm6=M_[6],mm7=M_[7],mm8=M_[8],mm9=M_[9],mm10=M_[10],mm11=M_[11];
		var mm12=M_[12],mm13=M_[13],mm14=M_[14],mm15=M_[15];
		nn=mm15*(mm10*(mm0*mm5-mm1*mm4)+mm9*(mm2*mm4-mm0*mm6)+mm8*(mm1*mm6-mm2*mm5));
		nn-=mm14*(mm11*(mm0*mm5-mm1*mm4)+mm9*(mm3*mm4-mm0*mm7)+mm8*(mm1*mm7-mm3*mm5));
		nn+=mm13*(mm11*(mm0*mm6-mm2*mm4)+mm10*(mm3*mm4-mm0*mm7)+mm8*(mm2*mm7-mm3*mm6));
		nn-=mm12*(mm11*(mm1*mm6-mm2*mm5)+mm10*(mm3*mm5-mm1*mm7)+mm9*(mm2*mm7-mm3*mm6));
		M_[0]=(mm15*(mm5*mm10-mm6*mm9)+mm14*(mm7*mm9-mm5*mm11)+mm13*(mm6*mm11-mm7*mm10))/nn;
		M_[1]=-(mm15*(mm1*mm10-mm2*mm9)+mm14*(mm3*mm9-mm1*mm11)+mm13*(mm2*mm11-mm3*mm10))/nn;
		M_[2]=(mm15*(mm1*mm6-mm2*mm5)+mm14*(mm3*mm5-mm1*mm7)+mm13*(mm2*mm7-mm3*mm6))/nn;
		M_[3]=-(mm11*(mm1*mm6-mm2*mm5)+mm10*(mm3*mm5-mm1*mm7)+mm9*(mm2*mm7-mm3*mm6))/nn;
		M_[4]=-(mm15*(mm4*mm10-mm6*mm8)+mm14*(mm7*mm8-mm4*mm11)+mm12*(mm6*mm11-mm7*mm10))/nn;
		M_[5]=(mm15*(mm0*mm10-mm2*mm8)+mm14*(mm3*mm8-mm0*mm11)+mm12*(mm2*mm11-mm3*mm10))/nn;
		M_[6]=-(mm15*(mm0*mm6-mm2*mm4)+mm14*(mm3*mm4-mm0*mm7)+mm12*(mm2*mm7-mm3*mm6))/nn;
		M_[7]=(mm11*(mm0*mm6-mm2*mm4)+mm10*(mm3*mm4-mm0*mm7)+mm8*(mm2*mm7-mm3*mm6))/nn;
		M_[8]=(mm15*(mm4*mm9-mm5*mm8)+mm13*(mm7*mm8-mm4*mm11)+mm12*(mm5*mm11-mm7*mm9))/nn;
		M_[9]=-(mm15*(mm0*mm9-mm1*mm8)+mm13*(mm3*mm8-mm0*mm11)+mm12*(mm1*mm11-mm3*mm9))/nn;
		M_[10]=(mm15*(mm0*mm5-mm1*mm4)+mm13*(mm3*mm4-mm0*mm7)+mm12*(mm1*mm7-mm3*mm5))/nn;
		M_[11]=-(mm11*(mm0*mm5-mm1*mm4)+mm9*(mm3*mm4-mm0*mm7)+mm8*(mm1*mm7-mm3*mm5))/nn;
		M_[12]=-(mm14*(mm4*mm9-mm5*mm8)+mm13*(mm6*mm8-mm4*mm10)+mm12*(mm5*mm10-mm6*mm9))/nn;
		M_[13]=(mm14*(mm0*mm9-mm1*mm8)+mm13*(mm2*mm8-mm0*mm10)+mm12*(mm1*mm10-mm2*mm9))/nn;
		M_[14]=-(mm14*(mm0*mm5-mm1*mm4)+mm13*(mm2*mm4-mm0*mm6)+mm12*(mm1*mm6-mm2*mm5))/nn;
		M_[15]=(mm10*(mm0*mm5-mm1*mm4)+mm9*(mm2*mm4-mm0*mm6)+mm8*(mm1*mm6-mm2*mm5))/nn;
	}
	else
	{
		a=0;
		var nn=n*n,nn1=nn+1;nnn=nn1*nn;
		M=new Array(nnn),v=new Array(nn);//alert(m.join("  "));
		for(i=0;i<nnn;i++)M[i]=0;
		for(i=0;i<n;i++)M[nn+i*(n+1)*nn1]=1;//alert(M.join("   "));
		for(i=0;i<n;i++)
		{
			for(j=0;j<n;j++)
			{
				for(k=0;k<n;k++)
				{
					M[k*n+j+nn1*(i*n+j)]=M_[i*n+k];
//alert(i+"  "+j+"  "+k+"\nM["+(k*n+j+nn1*(i*n+j))+"]   M_["+(i*(n+1)+k)+"]\n"+M.join("   "));
				}
			}
		}
		//alert(M.join("  "));
		GLSL(M,M_,nn);
		//alert(M.join("  "));
	}
}
//if(window.location.href.indexOf("arndt-bruenner")==-1)solveGS=null;
function JacobiMatrix(m,x,fff,n)
{
	var i,j,xx,ff,d=0.0000000001;
	for(i=0;i<n;i++)
		for(j=0;j<n;j++)
		{
			xx=x[j];
			x[j]=xx+d;
			ff=f_(i,x);
			x[j]=xx;
			m[i*n+j]=(ff-fff[i])/d;
		}
}