#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;
}