#include "iostream.h"
#include "math.h"
#include "stdlib.h"
int sgn(double x)
{
if (x>0)
{
return 1;
}
else
{
if(x {
return -1;
}
}
return 0;
}
double bessj0(double x)
{
double p1 = 1.0;
double p2 = -0.001098628627;
double p3 = 0.00002734510407;
double p4 = -0.000002073370639;
double p5 = 2.093887211e-07;
double q1 = -0.01562499995;
double q2 = 0.0001430488765;
double q3 = -0.000006911147651;
double q4 = 7.621095161e-07;
double q5 = -9.34945152e-08;
double r1 = 57568490574.0;
double r2 = -13362590354.0;
double r3 = 651619640.7;
double r4 = -11214424.18;
double r5 = 77392.33017;
double r6 = -184.9052456;
double s1 = 57568490411.0;
double s2 = 1029532985.0;
double s3 = 9494680.718;
double s4 = 59272.64853;
double s5 = 267.8532712;
double s6 = 1.0;
double y,aaa,bbb,ccc,ddd,eee,ax,xx,z;
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)));
return 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));
return sqrt(0.636619772 / ax) * (cos(xx) * aaa - eee);
}
}
double bessj1( double x)
{
double r1,r2,r3,r4,r5,r6,s1,s2,s3,s4,s5,s6,p1,p2,p3,p4,p5,q1,q2,q3,q4,q5;
double y,aaa,bbb,ccc,z,ax,xx;
r1 = 72362614232.0;
r2 = -7895059235.0;
r3 = 242396853.1;
r4 = -2972611.439;
r5 = 15704.4826;
r6 = -30.16036606;
s1 = 144725228442.0;
s2 = 2300535178.0;
s3 = 18583304.74;
s4 = 99447.43394;
s5 = 376.9991397;
s6 = 1.0;
p1 = 1.0;
p2 = 0.00183105;
p3 = -0.00003516396496;
p4 = 0.000002457520174;
p5 = -0.000000240337019;
q1 = 0.04687499995;
q2 = -0.0002002690873;
q3 = 0.000008449199096;
q4 = -0.00000088228987;
q5 = 0.000000105787412;
if (fabs(x) < 8.0)
{
y = x * x;
aaa = r1 + y * (r2 + y * (r3 + y * (r4 + y * (r5 + y * r6))));
bbb = s1 + y * (s2 + y * (s3 + y * (s4 + y * (s5 + y * s6))));
return x * aaa / bbb;
}
else
{
ax = fabs(x);
z = 8.0 / ax;
y = z * z;
xx = ax - 2.356194491;
aaa = p1 + y * (p2 + y * (p3 + y * (p4 + y * p5)));
bbb = q1 + y * (q2 + y * (q3 + y * (q4 + y * q5)));
ccc = sqrt(0.636619772 / ax);
return ccc * (cos(xx) * aaa - z * sin(xx) * bbb * sgn(x));
}
}
double func(double x)
{
return bessj0(x);
}
double deriv(double x)
{
return -bessj1(x);
}
void mnbrak(double& ax, double& bx, double& cx, double& fa, double& fb, double& fc)
{
double r,q,dum,gold = 1.618034;
int glimit = 100;
double u,ulim,fu,tiny = 1e-20;
fa = func(ax);
fb = func(bx);
if (fb > fa)
{
dum = ax;
ax = bx;
bx = dum;
dum = fb;
fb = fa;
fa = dum;
}
cx = bx + gold * (bx - ax);
fc = func(cx);
while (fb >= fc)
{
r = (bx - ax) * (fb - fc);
q = (bx - cx) * (fb - fa);
dum = q - r;
if (fabs(dum) < tiny)
{
dum = tiny;
}
u = bx - ((bx - cx) * q - (bx - ax) * r) / (2 * dum);
ulim = bx + glimit * (cx - bx);
if ((bx - u) * (u - cx) > 0)
{
fu = func(u);
if (fu < fc)
{
ax = bx;
fa = fb;
bx = u;
fb = fu;
return;
}
else
{
if (fu > fb)
{
cx = u;
fc = fu;
return;
}
}
u = cx + gold * (cx - bx);
fu = func(u);
}
else
{
if ((cx - u) * (u - ulim) > 0)
{
fu = func(u);
if (fu < fc)
{
bx = cx;
cx = u;
u = cx + gold * (cx - bx);
fb = fc;
fc = fu;
fu = func(u);
}
}
else
{
if ((u - ulim) * (ulim - cx) >= 0)
{
u = ulim;
fu = func(u);
}
else
{
u = cx + gold * (cx - bx);
fu = func(u);
}
}
}
ax = bx;
bx = cx;
cx = u;
fa = fb;
fb = fc;
fc = fu;
}
}
double dbrent(double ax, double bx, double cx, double tol, double& xmin)
{
int done,iter,itmax = 100;
double d,fu,xm,tol1,tol2,a,b,cgold = 0.381966;
double u,v,w,x,e,fx,fv1,fw,zeps = 0.0000000001;
double dx,dv,dw,d1,d2,u1,u2,olde,du;
bool ok1,ok2;
a = ax;
if (cx < ax)
{
a = cx;
}
b = ax;
if (cx > ax)
{
b = cx;
}
v = bx;
w = v;
x = v;
e = 0.0;
fx = func(x);
fv1 = fx;
fw = fx;
dx = deriv(x);
dv = dx;
dw = dx;
for (iter = 1; iter {
xm = 0.5 * (a + b);
tol1 = tol * fabs(x) + zeps;
tol2 = 2.0 * tol1;
if (fabs(x - xm) {
done = -1;
break;
}
else
{
done = 0;
}
if (fabs(e) > tol1)
{
d1 = 2.0 * (b - a);
d2 = d1;
if (dw != dx)
{
d1 = (w - x) * dx / (dx - dw);
}
if (dv != dx)
{
d2 = (v - x) * dx / (dx - dv);
}
u1 = x + d1;
u2 = x + d2;
ok1 = ((a - u1) * (u1 - b) > 0.0) && (dx * d1 ok2 = ((a - u2) * (u2 - b) > 0.0) && (dx * d2 olde = e;
e = d;
if (ok1 || ok2)
{
if (ok1 && ok2)
{
d = d1;
}
else
{
d = d2;
}
}
else
{
if (ok1)
{
d = d1;
}
else
{
d = d2;
}
}
if (fabs(d) > fabs(0.5 * olde))
{
u = x + d;
if (u - a < tol2 || b - u < tol2)
{
d = fabs(tol1) * sgn(xm - x);
}
}
}
if (dx >= 0.0)
{
e = a - x;
}
else
{
e = b - x;
}
d = 0.5 * e;
if (fabs(d) >= tol1)
{
u = x + d;
fu = func(u);
}
else
{
u = x + fabs(tol1) * sgn(d);
fu = func(u);
if (fu > fx)
{
done = -1;
break;
}
else
{
done = 0;
}
}
du = deriv(u);
if (fu {
if (u >= x)
{
a = x;
}
else
{
b = x;
}
v = w;
fv1 = fw;
dv = dw;
w = x;
fw = fx;
dw = dx;
x = u;
fx = fu;
dx = du;
}
else
{
if (u < x)
{
a = u;
}
else
{
b = u;
}
if (fu {
v = w;
fv1 = fw;
dv = dw;
w = u;
fw = fu;
dw = du;
}
else
{
if (fu {
v = u;
fv1 = fu;
dv = du;
}
}
}
}
if (!done)
{
cout exit(1);
}
else
{
xmin = x;
return fx;
}
}
void main()
{
//program d11r4
//driver for routine dbrent
double dbr,ax,bx,cx,fa,fb,fc,tol = 0.000001;
double xmin;
double eql = 0.0001;
double amin[21];
int iflag,i,nmin = 0;
cout cout cout cout cout.setf(ios::fixed|ios::left);
cout.precision(6);
for (i = 1; i {
ax = i;
bx = i + 1.0;
mnbrak(ax, bx, cx, fa, fb, fc);
dbr = dbrent(ax, bx, cx, tol, xmin);
if (nmin == 0)
{
amin[1] = xmin;
nmin = 1;
cout cout.width(14); cout cout.width(14); cout cout.width(14); cout cout.width(14); cout cout }
else
{
iflag = 0;
for (int j = 1; j {
if (fabs(xmin - amin[j]) {
iflag = 1;
}
}
if (iflag == 0)
{
nmin = nmin + 1;
amin[nmin] = xmin;
cout cout.width(14); cout cout.width(14); cout cout.width(14); cout cout.width(14); cout cout }
}
}
}
