#define _USE_MATH_DEFINES
#include <float.h>
#include <math.h>
#include <stdio.h>
#include <algorithm>
#include <assert.h>
#define DIFFSIGN(a,b) (a <= 0 && b >= 0 || a >= 0 && b <= 0)
double fun(double x)
{
return sin(x);
}
bool between(double x, double a, double b)
{
if (b > a)
return (x >= a && x <= b);
else
return (x >= b && x <= a);
}
double hfun(double b, double c)
{
assert(b != c);
if (c > b)
return b + fabs(b*DBL_EPSILON);
else
return b - fabs(b*DBL_EPSILON);
}
double mfun(double b, double c)
{
return 0.5*(b + c);
}
double lfun(double b, double a, double fb, double fa)
{
if (fb != fa)
return b - fb*(b - a) / (fb - fa); //secant method;
else if (fa != 0)
return INFINITY;
else return b;
}
double ffun(double a, double b, double fa, double fb)
{
assert(a != b);
return (fa - fb) / (a - b);
}
double rfun(double b, double a, double d, double fb, double fa, double fd)
{
double alpha = ffun(b, d, fb, fd)*fa;
double beta = ffun(a, d, fa, fd)*fb;
if (beta != alpha)
{
return b - beta*(b - a) / (beta - alpha);
}
else if (alpha != 0)
return INFINITY;
else
return 0; //beta == alpha == 0
}
double wfun(double l, double b, double c)
{
double h = hfun(b, c);
double m = mfun(b, c);
if (between(l, h, m))
{
return l;
}
else if (fabs(l - b) <= fabs(b*DBL_EPSILON) && !between(l, b, m))
{
return h;
}
else
{
return m;
}
}
double Dekker(double x0, double x1, double delta)
{
double x, xp, fx, fxp, a, b, c, fa, fb;
double lambda, rho, d, fd;
int itercnt = 0;
//initialization
fxp = fun(x0);
fx = fun(x1);
if (x0 == x1) return fx;
if (!DIFFSIGN(fx, fxp)) return NAN;
if (fabs(fx) <= fabs(fxp))
{
b = x1;
a = c = x0;
fa = fxp;
fb = fx;
}
else
{
b = x0;
a = c = x1;
fa = fx;
fb = fxp;
}
double xk = xp = x0;
double fxk = fxp;
x = x1;
d = NAN;
fd = NAN;
itercnt = 1;
int age = 0;
double bp = b;
double cp = c;
double ap = a;
double fap, fbp;
while (fabs(b - c) > 2 * delta)
{
itercnt++;
age++;
if (fabs(b - c) <= (0.5 + 2 * DBL_EPSILON)*(fabs(bp - cp) + fabs(b)*DBL_EPSILON))
age = 1;
xp = x;
if (itercnt == 2)
{
lambda = lfun(b, a, fb, fa);
if (fabs(lambda - b) < fabs(b*DBL_EPSILON)) break;
x = wfun(lambda, b, c);
}
else if (itercnt >= 3 && age <= 3)
{
rho = rfun(b, a, d, fb, fa, fd);
if (fabs(rho - b) < fabs(b*DBL_EPSILON)) break;
x = wfun(rho, b, c);
}
else if (itercnt >= 3 && age == 4)
{
rho = rfun(b, a, d, fb, fa, fd);
if (fabs(rho - b) < fabs(b*DBL_EPSILON)) break;
x = wfun(2 * rho - b, b, c);
}
else
{
x = mfun(b, c);
}
fxp = fx;
fx = fun(x);
if (DIFFSIGN(fxp, fx))
{
xk = xp;
fxk = fxp;
}
bp = b;
fbp = fb;
ap = a;
fap = fa;
cp = c;
if (fabs(fx) <= fabs(fxk))
{
a = b;
fa = fb;
b = x;
fb = fx;
c = xk;
}
else
{
b = xk;
fb = fxk;
a = c = x;
fa = fx;
}
if (b == x || b == bp)
{
d = ap;
fd = fap;
}
else
{
d = bp;
fd = fbp;
}
}
return b;
}
int main()
{
double res = Dekker(1, 10, 1e-12);
printf("%g\n", res);
return 0;
}
#include <stdio.h>
#include <stdlib.h>
#include <mpfr.h>
#define BITS_PER_DIGIT 3.32192809488736234787
#define DOUBLE_PREC 53
#define QUADRUPLE_PREC 113
#define OCTUPLE_PREC 237
#define BIT_PREC (DOUBLE_PREC)
void fun(mpfr_t *result, mpfr_t x)
{
mpfr_sin(*result, x, MPFR_RNDN);
}
int between(mpfr_t x, mpfr_t a, mpfr_t b)
{
if (mpfr_greater_p(b,a))
return (mpfr_greaterequal_p(x,a) && mpfr_lessequal_p(x,b));
else
return (mpfr_greaterequal_p(x,b) && mpfr_lessequal_p(x,a));
}
void hfun(mpfr_t *result, mpfr_t b, mpfr_t c, mpfr_t machineEps)
{
mpfr_t eps;
mpfr_init2(eps, BIT_PREC);
mpfr_mul(eps, b, machineEps, MPFR_RNDN);
mpfr_abs(eps, eps, MPFR_RNDN);
if (mpfr_greater_p(c,b))
mpfr_add(*result, b, eps, MPFR_RNDN);
else
mpfr_sub(*result, b, eps, MPFR_RNDN);
mpfr_clear(eps);
}
void mfun(mpfr_t *result, mpfr_t b, mpfr_t c)
{
mpfr_add(*result, b, c, MPFR_RNDN);
mpfr_div_2ui(*result, *result, 1, MPFR_RNDN);
}
void lfun(mpfr_t *result, mpfr_t b, mpfr_t a, mpfr_t fb, mpfr_t fa)
{
if (mpfr_lessgreater_p(fb,fa))
{
mpfr_t tmp;
mpfr_init2(tmp, BIT_PREC);
mpfr_sub(tmp, fb, fa, MPFR_RNDN);
mpfr_sub(*result, b, a, MPFR_RNDN);
mpfr_mul(*result, fb, *result, MPFR_RNDN);
mpfr_div(*result, *result, tmp, MPFR_RNDN);
mpfr_sub(*result, b, *result, MPFR_RNDN);
mpfr_clear(tmp);
}
else
{
mpfr_t zero;
mpfr_init2(zero, BIT_PREC);
mpfr_set_d(zero, 0, MPFR_RNDN);
if (mpfr_lessgreater_p(fa,zero))
mpfr_set_d(*result, 1.0/0.0, MPFR_RNDN);
else
mpfr_set(*result, b, MPFR_RNDN);
mpfr_clear(zero);
}
}
void ffun(mpfr_t *result, mpfr_t a, mpfr_t b, mpfr_t fa, mpfr_t fb)
{
mpfr_t tmp;
mpfr_init2(tmp, BIT_PREC);
mpfr_sub(tmp, a, b, MPFR_RNDN);
mpfr_sub(*result, fa, fb, MPFR_RNDN);
mpfr_div(*result, *result, tmp, MPFR_RNDN);
mpfr_clear(tmp);
}
void rfun(mpfr_t *result, mpfr_t b, mpfr_t a, mpfr_t d, mpfr_t fb, mpfr_t fa, mpfr_t fd)
{
mpfr_t alpha, beta;
mpfr_init2(alpha, BIT_PREC);
mpfr_init2(beta, BIT_PREC);
ffun(&alpha, b, d, fb, fd);
mpfr_mul(alpha, alpha, fa, MPFR_RNDN);
ffun(&beta, a, d, fa, fd);
mpfr_mul(beta, beta, fb, MPFR_RNDN);
if (mpfr_lessgreater_p(beta,alpha))
{
mpfr_t tmp;
mpfr_init2(tmp, BIT_PREC);
mpfr_sub(tmp, beta, alpha, MPFR_RNDN);
mpfr_sub(*result, b, a, MPFR_RNDN);
mpfr_mul(*result, beta, *result, MPFR_RNDN);
mpfr_div(*result, *result, tmp, MPFR_RNDN);
mpfr_sub(*result, b, *result, MPFR_RNDN);
mpfr_clear(tmp);
}
else
{
mpfr_t zero;
mpfr_init2(zero, BIT_PREC);
mpfr_set_d(zero, 0, MPFR_RNDN);
if (mpfr_lessgreater_p(alpha,zero))
mpfr_set_d(*result, 1.0/0.0, MPFR_RNDN);
else
mpfr_set_d(*result, 0, MPFR_RNDN);
mpfr_clear(zero);
}
mpfr_clear(alpha);
mpfr_clear(beta);
}
void wfun(mpfr_t *result, mpfr_t l, mpfr_t b, mpfr_t c, mpfr_t machineEps)
{
mpfr_t h,m;
mpfr_init2(h, BIT_PREC);
mpfr_init2(m, BIT_PREC);
hfun(&h, b, c, machineEps);
mfun(&m, b, c);
if (between(l, h, m))
{
mpfr_set(*result, l, MPFR_RNDN);
}
else
{
mpfr_t tmp1,tmp2;
mpfr_init2(tmp1, BIT_PREC);
mpfr_init2(tmp2, BIT_PREC);
mpfr_sub(tmp1, l, b, MPFR_RNDN);
mpfr_abs(tmp1, tmp1, MPFR_RNDN);
mpfr_mul(tmp2, b, machineEps, MPFR_RNDN);
mpfr_abs(tmp2, tmp2, MPFR_RNDN);
if (mpfr_lessequal_p(tmp1, tmp2) && !between(l, b, m))
{
mpfr_set(*result, h, MPFR_RNDN);
}
else
{
mpfr_set(*result, m, MPFR_RNDN);
}
mpfr_clear(tmp1);
mpfr_clear(tmp2);
}
mpfr_clear(h);
mpfr_clear(m);
}
int DiffSign(mpfr_t a,mpfr_t b)
{
mpfr_t zero;
mpfr_init2(zero, BIT_PREC);
mpfr_set_d(zero, 0, MPFR_RNDN);
int res = (mpfr_lessequal_p(a,zero) && mpfr_greaterequal_p(b,zero))
|| (mpfr_greaterequal_p(a,zero) && mpfr_lessequal_p(b,zero));
mpfr_clear(zero);
return res;
}
void Dekker(mpfr_t *result, mpfr_t x0, mpfr_t x1, mpfr_t delta)
{
mpfr_t x, xp, fx, fxp, a, b, c, fa, fb;
mpfr_t tmp1,tmp2,tmp3, tmp4;
int itercnt = 0;
mpfr_init2(fx, BIT_PREC);
fun(&fx, x1);
if (mpfr_equal_p(x0,x1))
{
mpfr_set(*result, fx, MPFR_RNDN);
mpfr_clear(fx);
return;
}
mpfr_init2(fxp, BIT_PREC);
fun(&fxp, x0);
if (!DiffSign(fx, fxp))
{
mpfr_set_d(*result, 0.0/0.0, MPFR_RNDN);
mpfr_clear(fx);
mpfr_clear(fxp);
return;
}
mpfr_init2(x, BIT_PREC);
mpfr_init2(xp, BIT_PREC);
mpfr_init2(a, BIT_PREC);
mpfr_init2(b, BIT_PREC);
mpfr_init2(c, BIT_PREC);
mpfr_init2(fa, BIT_PREC);
mpfr_init2(fb, BIT_PREC);
mpfr_t lambda, rho, d, fd;
mpfr_init2(lambda, BIT_PREC);
mpfr_init2(rho, BIT_PREC);
mpfr_init2(d, BIT_PREC);
mpfr_init2(fd, BIT_PREC);
mpfr_init2(tmp1, BIT_PREC);
mpfr_init2(tmp2, BIT_PREC);
mpfr_init2(tmp3, BIT_PREC);
mpfr_init2(tmp4, BIT_PREC);
mpfr_abs(tmp1, fx, MPFR_RNDN);
mpfr_abs(tmp2, fxp, MPFR_RNDN);
if (mpfr_lessequal_p(tmp1, tmp2))
{
mpfr_set(b, x1, MPFR_RNDN);
mpfr_set(c, x0, MPFR_RNDN);
mpfr_set(a, c, MPFR_RNDN);
mpfr_set(fa, fxp, MPFR_RNDN);
mpfr_set(fb, fx, MPFR_RNDN);
}
else
{
mpfr_set(b, x0, MPFR_RNDN);
mpfr_set(c, x1, MPFR_RNDN);
mpfr_set(a, c, MPFR_RNDN);
mpfr_set(fa, fx, MPFR_RNDN);
mpfr_set(fb, fxp, MPFR_RNDN);
}
mpfr_t xk,fxk;
mpfr_init2(xk, BIT_PREC);
mpfr_set(xp, x0, MPFR_RNDN);
mpfr_set(xk, xp, MPFR_RNDN);
mpfr_init2(fxk, BIT_PREC);
mpfr_set(fxk, fxp, MPFR_RNDN);
mpfr_set(x, x1, MPFR_RNDN);
mpfr_set_d(d, 0.0/0.0, MPFR_RNDN);
mpfr_set_d(fd, 0.0/0.0, MPFR_RNDN);
itercnt = 1;
int age = 0;
mpfr_t bp, cp, ap, fap, fbp;
mpfr_init2(bp, BIT_PREC);
mpfr_init2(cp, BIT_PREC);
mpfr_init2(ap, BIT_PREC);
mpfr_init2(fap, BIT_PREC);
mpfr_init2(fbp, BIT_PREC);
mpfr_set(bp, b, MPFR_RNDN);
mpfr_set(cp, c, MPFR_RNDN);
mpfr_set(ap, a, MPFR_RNDN);
mpfr_t machineEps;
mpfr_init2(machineEps, BIT_PREC);
mpfr_set_d(machineEps, 1, MPFR_RNDN);
mpfr_div_2ui(machineEps, machineEps, BIT_PREC-1, MPFR_RNDN);
while (1)
{
mpfr_sub(tmp1, b, c, MPFR_RNDN);
mpfr_abs(tmp1, tmp1, MPFR_RNDN);
mpfr_mul_2ui(tmp2, delta, 1, MPFR_RNDN);
if (mpfr_lessequal_p(tmp1, tmp2)) break;
itercnt++;
age++;
mpfr_mul_2ui(tmp2, machineEps, 1, MPFR_RNDN);
mpfr_add_d(tmp2, tmp2, 0.5, MPFR_RNDN);
mpfr_sub(tmp3, bp, cp, MPFR_RNDN);
mpfr_abs(tmp3, tmp3, MPFR_RNDN);
mpfr_abs(tmp4, b, MPFR_RNDN);
mpfr_mul(tmp4, tmp4, machineEps, MPFR_RNDN);
mpfr_add(tmp3, tmp3, tmp4, MPFR_RNDN);
mpfr_mul(tmp2, tmp2, tmp3, MPFR_RNDN);
if (mpfr_lessequal_p(tmp1, tmp2))
age = 1;
mpfr_set(xp, x, MPFR_RNDN);
if (itercnt == 2)
{
lfun(&lambda, b, a, fb, fa);
mpfr_sub(tmp1, lambda,b, MPFR_RNDN);
mpfr_abs(tmp1, tmp1, MPFR_RNDN);
mpfr_mul(tmp2, b, machineEps, MPFR_RNDN);
mpfr_abs(tmp2, tmp2, MPFR_RNDN);
if (mpfr_less_p(tmp1, tmp2)) break;
wfun(&x, lambda, b, c, machineEps);
}
else if (itercnt >= 3 && age <= 3)
{
rfun(&rho, b, a, d, fb, fa, fd);
mpfr_sub(tmp1, rho, b, MPFR_RNDN);
mpfr_abs(tmp1, tmp1, MPFR_RNDN);
mpfr_mul(tmp2, b, machineEps, MPFR_RNDN);
mpfr_abs(tmp2, tmp2, MPFR_RNDN);
if (mpfr_less_p(tmp1, tmp2)) break;
wfun(&x, rho, b, c, machineEps);
}
else if (itercnt >= 3 && age == 4)
{
rfun(&rho, b, a, d, fb, fa, fd);
mpfr_sub(tmp1, rho, b, MPFR_RNDN);
mpfr_abs(tmp1, tmp1, MPFR_RNDN);
mpfr_mul(tmp2, b, machineEps, MPFR_RNDN);
mpfr_abs(tmp2, tmp2, MPFR_RNDN);
if (mpfr_less_p(tmp1, tmp2)) break;
mpfr_mul_2ui(tmp1, rho, 1, MPFR_RNDN);
mpfr_sub(tmp1, tmp1, b, MPFR_RNDN);
wfun(&x, tmp1, b, c, machineEps);
}
else
{
mfun(&x, b, c);
}
mpfr_set(fxp, fx, MPFR_RNDN);
fun(&fx, x);
if (DiffSign(fxp, fx))
{
mpfr_set(xk, xp, MPFR_RNDN);
mpfr_set(fxk, fxp, MPFR_RNDN);
}
mpfr_set(bp, b, MPFR_RNDN);
mpfr_set(fbp, fb, MPFR_RNDN);
mpfr_set(ap, a, MPFR_RNDN);
mpfr_set(fap, fa, MPFR_RNDN);
mpfr_set(cp, c, MPFR_RNDN);
mpfr_abs(tmp1, fx, MPFR_RNDN);
mpfr_abs(tmp2, fxk, MPFR_RNDN);
if (mpfr_lessequal_p(tmp2, tmp2))
{
mpfr_set(a, b, MPFR_RNDN);
mpfr_set(fa, fb, MPFR_RNDN);
mpfr_set(b, x, MPFR_RNDN);
mpfr_set(fb, fx, MPFR_RNDN);
mpfr_set(c, xk, MPFR_RNDN);
}
else
{
mpfr_set(b, xk, MPFR_RNDN);
mpfr_set(fb, fxk, MPFR_RNDN);
mpfr_set(c, x, MPFR_RNDN);
mpfr_set(a, c, MPFR_RNDN);
mpfr_set(fa, fx, MPFR_RNDN);
}
if (mpfr_equal_p(b, x) || mpfr_equal_p(b, bp))
{
mpfr_set(d, ap, MPFR_RNDN);
mpfr_set(fd, fap, MPFR_RNDN);
}
else
{
mpfr_set(d, bp, MPFR_RNDN);
mpfr_set(fd, fbp, MPFR_RNDN);
}
}
mpfr_set(*result, b, MPFR_RNDN);
mpfr_clear(machineEps);
mpfr_clear(xk);
mpfr_clear(fxk);
mpfr_clear(tmp1);
mpfr_clear(tmp2);
mpfr_clear(tmp3);
mpfr_clear(tmp4);
mpfr_clear(x);
mpfr_clear(xp);
mpfr_clear(fx);
mpfr_clear(fxp);
mpfr_clear(a);
mpfr_clear(b);
mpfr_clear(c);
mpfr_clear(fa);
mpfr_clear(fb);
mpfr_clear(lambda);
mpfr_clear(rho);
mpfr_clear(d);
mpfr_clear(fd);
}
int main()
{
mpfr_t a,b,t, res;
mpfr_init2(a, BIT_PREC);
mpfr_set_d(a, 1, MPFR_RNDN);
mpfr_init2(b, BIT_PREC);
mpfr_set_d(b, 10, MPFR_RNDN);
mpfr_init2(t, BIT_PREC);
mpfr_set_d(t, 1e-16, MPFR_RNDN);
mpfr_init2(res, BIT_PREC);
Dekker(&res, a, b, t);
printf("result= ");
mpfr_out_str (stdout, 10, 0, res, MPFR_RNDD);
printf("\n");
return 0;
}
