#include <process.h>
#include <math.h>
#include <iomanip.h>
#include <iostream.h>

double bessj0(double x)
{
	double p1,p2,p3,p4,p5,q1,q2,q3,q4,q5;
	double r1,r2,r3,r4,r5,r6,s1,s2,s3,s4,s5,s6;
	double y,bbb,ccc,aaa,temp,eee,ddd,ax,xx,z;
    p1 = 1.0;                  p2 = -0.001098628627;
    p3 = 0.00002734510407;    p4 = -0.000002073370639;
    p5 = 2.093887211e-07;
    q1 = -0.01562499995;      q2 = 0.0001430488765;
    q3 = -0.000006911147651;  q4 = 7.621095161e-07;
    q5 = -9.34945152e-08;
    r1 = 57568490574.0;        r2 = -13362590354.0;
    r3 = 651619640.7;         r4 = -11214424.18;
    r5 = 77392.33017;         r6 = -184.9052456;
    s1 = 57568490411.0;        s2 = 1029532985.0;
    s3 = 9494680.718;         s4 = 59272.64853;
    s5 = 267.8532712;        s6 = 1.0;
    if (fabs(x) < 8.0)
	{
       y = x * x;
       bbb = y * (r4 + y * (r5 + y * r6));
       aaa = r1 + y * (r2 + y * (r3 + bbb));
       ccc = y * (s3 + y * (s4 + y * (s5 + y * s6)));
       temp = aaa / (s1 + y * (s2 + ccc));
	}
    else
	{
       ax = fabs(x);
       z = 8.0 / ax;
       y = z * z;
       xx = ax - 0.785398164;
       ccc = y * (p3 + y * (p4 + y * p5));
       aaa = p1 + y * (p2 + ccc);
       ddd = y * (q3 + y * (q4 + y * q5));
       eee = z * sin(xx) * (q1 + y * (q2 + ddd));
       temp = sqrt(0.636619772 / ax) * (cos(xx) * aaa - eee);
	}
	return temp;
}

double  func(double x)
{
 double t;
 t=bessj0( x);
 return t;
}

void zbrak(double x1, double x2, int n, double xb1[], double xb2[], int& nb)
{
	int nbb,i;
	double dx,x,fc,fp;
    nbb = nb;
    nb = 0;
    x = x1;
    dx = (x2 - x1)/ n;
    fp = func(x);
    for (i = 1 ;i<=n;i++)
	{
        x = x + dx;
        fc = func(x);
        if ((fc * fp) < 0.0 )
		{
            nb = nb + 1;
            xb1[nb] = x - dx;
            xb2[nb] = x;
        }
        fp = fc;
        if (nbb == nb) exit(1);
    }
}

double zbrent(double x1, double x2, double tol)
{
	double a,b,fa,fb,fc,c,d,e,tol1,xm,q,r,p,eps,s,aaa,tt;
	int t,itmax,iter;
    itmax = 100;
    eps = 0.00000003;
    a = x1;
    b = x2;
    fa = func(a);
    fb = func(b);
    if ((fb * fa) > 0.0)
	{
        cout<< "root must be bracketed for zbrent"<<endl;
    }
    fc = fb;
    for (iter = 1;iter<=itmax;iter++)
	{
        if ((fb * fc) > 0.0)
		{
            c = a;
            fc = fa;
            d = b - a;
            e = d;
        }
        if (fabs(fc) < fabs(fb))
		{
            a = b;
            b = c;
            c = a;
            fa = fb;
            fb = fc;
            fc = fa;
        }
        tol1 = 2.0 * eps * fabs(b) + 0.5 * tol;
        xm = 0.5 * (c - b);
        if ((fabs(xm) <= tol1) || (fb == 0.0))
		{
			tt=b;
            return tt;
			exit(1);
            
        }
        if ((fabs(e) >= tol1) && (fabs(fa) > fabs(fb)))
		{
            s = fb / fa;
            if (a == c)
			{
                p = 2.0 * xm * s;
                q = 1 - 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 = fabs(p);
            if ((3.0 * xm * q - fabs(tol1 * q)) < fabs(e * q))
			{
                aaa = 3.0 * xm * q - fabs(tol1 * q);
			}
            else
			{
                aaa = fabs(e * q);
            }
            if ((2.0 * p) < aaa)
			{
                e = d;
                d = p / q;
			}
            else
			{
                d = xm;
                e = d;
            }
		}
        else
		{
            d = xm;
            e = d;
        }
        a = b;
        fa = fb;
        if (fabs(d) > tol1)
            b = b + d;
        else
		{
			if (xm>0)
				t=1;
			if (xm==0) 
				t=0;
			if (xm<0) 
				t=-1;
            b = b + fabs(tol1) * t;
		}
        
        fb = func(b);
	}
	 tt=b;
	 return tt;
     cout<< "zbrent exceeding maximum iterations."<<endl;
}

void main()
{
   //program d10r7
   //driver for routine zbrent
	int n,nbmax,nb,i;
	double  x1,x2,tol,root,xb1[21], xb2[21];
    n = 100;
    nbmax = 20;
    x1 = 1.0;
    x2 = 50.0;
    nb = nbmax;
    zbrak(x1, x2, n, xb1, xb2, nb);
    cout<<endl;
    cout<<"Roots of Bessj0:"<<endl;
    cout<<endl;
    cout<<"                  x              f(x)"<<endl;
    for (i = 1; i<=nb; i++) 
	{

        tol = (0.0000010) * (xb1[i] + xb2[i]) / 2.0;
        root =zbrent(xb1[i], xb2[i], tol);
        cout<< setw(8)<<"root "<<i;
        cout<< setw(12)<<root;
        cout<< setw(19)<<bessj0(root)<<endl;
    }
}