Kody źródłowe/Metoda Dekkera

C

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

C z biblioteką MPFR

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