///Bilgisayar programlama II ve Sayısal Çözümleme
// 
//  
// Brent yöntemi ile kök bulmak 
// ve kökün sınırlarını saptamak
// Dr. M. Turhan Çoban

import java.util.*;
import java.awt.*;
import java.applet.Applet;
import java.awt.event.*;
import javax.swing.*;
//============= Tanımı gereken fonksiyon ================
abstract class fi_xi
{
  // multifunction multi independent variable
  // vector of dependent variables are returned
  // example f[0]=x[0]+sin(x[1])
  //         f[1]=x[0]*x[0]-x[1]
  // func(x) returns the value of f[0] and f[1]
  // as a two dimensional vector

  abstract double[] func(double x[]);
}

abstract class fij_xi
{
  // multifunction multi independent variable
  // for n independent variable n*n matrix of
  // dependent variables are returned
  // example f[0][0]=x[0]+Math.sin(x[1])   f[0][1]=x[0]-x[1]
  //         f[1][0]=x[0]*x[0]-x[1]        f[1][1]=Math.exp(x[0]+x[1]*x[1]
  // func(x) returns the value of f[0][0], f[0][1]
  //                              f[1][0], f[1][1]
  // as a two dimensional matrix
  abstract double[][] func(double x[]);
}

class fa extends fi_xi
{
public double[] func(double x[])
{
//çözümü istenen fonksiyonun türev fonksiyonları	
//çözümü istenen fonksiyon	
double ff[]=new double[2];
ff[0]=x[0]*x[0]+x[0]*x[1]-10.0;
ff[1]=x[1]+3.0*x[0]*x[1]*x[1]-57.0;
return ff;
} 
}

class fb extends fij_xi
{
public double[][] func(double x[])
{
//çözümü istenen fonksiyonun türev fonksiyonları	
double a[][]=new double[2][2];
a[0][0]=4.0*x[0]+x[1];
a[0][1]=x[0];
a[1][0]=3.0*x[1]*x[1];
a[1][1]=1.0+6.0*x[0]*x[1];
return a;
} 
}

public class FSCO6
{

public static double f(double x)
{ return x*x*x-6.0*x*x+11.0*x-6.0;}
//üç gerçek kök 1 2 3

public static double[] genislet(double x1)
{return genislet(x1,(x1+0.1));}


public static double[] genislet(double x1,double x2)
{
	//kök olan bir bölge bulana kadar verilen sınırı genişletir
	int NTRY=50;
	double a[]=new double[2];
	double FACTOR=1.6;
	int j;
	double f1,f2;
	if (x1 == x2) System.out.println("Bad initial range in ZBRAC");
	f1=f(x1);
	f2=f(x2);
	for (j=1;j<=NTRY;j++) {
		if (f1*f2 < 0.0) {a[0]=x1;a[1]=x2;return a;}
		if (Math.abs(f1) < Math.abs(f2))
			f1=f(x1 += FACTOR*(x1-x2));
		else
			f2=f(x2 += FACTOR*(x2-x1));
	}
	{a[0]=x1;a[1]=x2;return a;}
}

public static double[][] koksinirisapta(double x1,double x2,int n,int nbb)
{
	// koklerin yer aldığı alt bölgeleri saptar
	// n : verilen bölgeyi böldüğümüz alt bölge sayısı
	// x1,x2 : sınır değerleri
	// nbb = aranan bölgedeki köksayısı
	//  
	int nb;
	int i;
	double x,fp,fc,dx;
	double xb[][]=new double[2][nbb];
	nb=0;
	dx=(x2-x1)/n;
	x=x1;
	fp=f(x1);
	for (i=1;i<=n;i++) 
	{
		x+=dx;
		fc=f(x);
		// eğer kök olan bölge bulunduysa.....
		if (fc*fp < 0.0 || fp==0) {
			xb[0][nb]=x-dx;
			xb[1][nb]=x;
			nb++;
		}
		fp=fc;
		if (nbb == nb) return xb;
	}	
if( nb == 0)
	JOptionPane.showMessageDialog(null,"arama tamamlandı kök olan bölge bulunamadı "
    ,"ARAMA SONUÇLARI UYARISI",JOptionPane.WARNING_MESSAGE);	
else if(nb<nbb)
  {	JOptionPane.showMessageDialog(null,"arama tamamlandı sadece "+nb+" adet kök bulundu \n"+
    "siz "+nbb+" adet kök için arama yaptırdınız","ARAMA SONUÇLARI UYARISI",JOptionPane.WARNING_MESSAGE);	
   double xc[][]=new double[2][nb];
   for (i=0;i<nb;i++) {xc[0][i]=xb[0][i];xc[1][i]=xb[1][i];}
   return xc;
  }
return xb;
}

public static double brent(double x1,double x2)
{
    //lineer olmıyan denklemin kökü : brent metodu
	//referans : Numerical Recipes in C, second edition, William H. Press, 
    //Saul A. Teukolsky, William T. Vetterling, Brian P. Flannery
    //cambridge university press	
	double ITMAX=200;
    double EPS=3.0e-10;
	int iter;
	double a=x1,b=x2,c,d,e,min1,min2;
	double tol=1.0e-5;
	c=0;
	d=0;
	e=0;
	double fa=f(a),fb=f(b),fc,p,q,r,s,tol1,xm;
	if (fb*fa > 0.0) JOptionPane.showMessageDialog(null,"Kök sınırları doğru olarak seçilmemiş \n sonuç hatalı olabilir"
	,"KÖK SINIRLARI SEÇİM UYARISI",JOptionPane.WARNING_MESSAGE);
	fc=fb;
	for (iter=1;iter<=ITMAX;iter++) {
		if (fb*fc > 0.0) {
			c=a;
			fc=fa;
			e=d=b-a;
		}
		if (Math.abs(fc) < Math.abs(fb)) {
			a=b;
			b=c;
			c=a;
			fa=fb;
			fb=fc;
			fc=fa;
		}
		tol1=2.0*EPS*Math.abs(b)+0.5*tol;
		xm=0.5*(c-b);
		if (Math.abs(xm) <= tol1 || fb == 0.0) return b;
		if (Math.abs(e) >= tol1 && Math.abs(fa) > Math.abs(fb)) {
			s=fb/fa;
			if (a == c) {
				p=2.0*xm*s;
				q=1.0-s;
			} else {
				q=fa/fc;
				r=fb/fc;
				p=s*(2.0*xm*q*(q-r)-(b-a)*(r-1.0));
				q=(q-1.0)*(r-1.0)*(s-1.0);
			}
			if (p > 0.0)  q = -q;
			p=Math.abs(p);
			min1=3.0*xm*q-Math.abs(tol1*q);
			min2=Math.abs(e*q);
			if (2.0*p < (min1 < min2 ? min1 : min2)) {
				e=d;
				d=p/q;
			} else {
				d=xm;
				e=d;
			}
		} else {
			d=xm;
			e=d;
		}
		a=b;
		fa=fb;
		if (Math.abs(d) > tol1)
			b += d;
		else
			b += (xm > 0.0 ? Math.abs(tol1) : -Math.abs(tol1));
		fb=f(b);
	}
	JOptionPane.showMessageDialog(null,"Uyarı maximum iterasyon sayısı aşıldı \n"+
   " çözüm geçerli olmıyabilir","MAKSİMUM ITERASYON SAYISI UYARISI",JOptionPane.WARNING_MESSAGE);
return b;
}

     public static void main (String args[])  
     { 
     double a[],x0,r;
     double dx=0.1;
     a=genislet(Double.parseDouble(JOptionPane.showInputDialog(" ilk tahmin değeri x0 : ")));
     r= brent(a[0],a[1]);
     JOptionPane.showMessageDialog(null," kök değeri  : "+r+"\nFonksiyon değeri : "+f(r)+"\nsınırlar: a="+a[0]+" b="+a[1],
     "brent yöntemiyle kök hesabı : ",JOptionPane.PLAIN_MESSAGE);
     double b[][]=koksinirisapta(0.0,10.0,10,4);
     System.out.println(Matrix.toString(b));
     b=koksinirisapta(b[0][0],b[1][0],10,4);
     System.out.println(Matrix.toString(b));
     System.exit(0);
     } 
}